[d82961] | 1 | /*
|
---|
| 2 | * Project: MoleCuilder
|
---|
| 3 | * Description: creates and alters molecular systems
|
---|
| 4 | * Copyright (C) 2011 University of Bonn. All rights reserved.
|
---|
| 5 | * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
|
---|
| 6 | */
|
---|
| 7 |
|
---|
| 8 | /*
|
---|
| 9 | * LinkedCell_Model.cpp
|
---|
| 10 | *
|
---|
| 11 | * Created on: Nov 15, 2011
|
---|
| 12 | * Author: heber
|
---|
| 13 | */
|
---|
| 14 |
|
---|
| 15 | // include config.h
|
---|
| 16 | #ifdef HAVE_CONFIG_H
|
---|
| 17 | #include <config.h>
|
---|
| 18 | #endif
|
---|
| 19 |
|
---|
| 20 | #include "CodePatterns/MemDebug.hpp"
|
---|
| 21 |
|
---|
| 22 | #include "LinkedCell_Model.hpp"
|
---|
| 23 |
|
---|
| 24 | #include <algorithm>
|
---|
| 25 | #include <boost/multi_array.hpp>
|
---|
| 26 | #include <limits>
|
---|
| 27 |
|
---|
| 28 | #include "Box.hpp"
|
---|
| 29 | #include "CodePatterns/Assert.hpp"
|
---|
| 30 | #include "CodePatterns/Info.hpp"
|
---|
| 31 | #include "CodePatterns/Log.hpp"
|
---|
| 32 | #include "LinearAlgebra/RealSpaceMatrix.hpp"
|
---|
| 33 | #include "LinearAlgebra/Vector.hpp"
|
---|
| 34 | #include "LinkedCell/IPointCloud.hpp"
|
---|
| 35 | #include "LinkedCell/LinkedCell.hpp"
|
---|
| 36 | #include "Atom/TesselPoint.hpp"
|
---|
| 37 |
|
---|
| 38 | #include "LinkedCell_Model_inline.hpp"
|
---|
| 39 |
|
---|
| 40 | // initialize static entities
|
---|
| 41 | LinkedCell::tripleIndex LinkedCell::LinkedCell_Model::NearestNeighbors;
|
---|
| 42 |
|
---|
| 43 | /** Constructor of LinkedCell_Model.
|
---|
| 44 | *
|
---|
| 45 | * @param radius desired maximum neighborhood distance
|
---|
| 46 | * @param _domain Box instance with domain size and boundary conditions
|
---|
| 47 | */
|
---|
| 48 | LinkedCell::LinkedCell_Model::LinkedCell_Model(const double radius, const Box &_domain) :
|
---|
| 49 | internal_Sizes(NULL),
|
---|
| 50 | domain(_domain)
|
---|
| 51 | {
|
---|
| 52 | // set default argument
|
---|
| 53 | NearestNeighbors[0] = NearestNeighbors[1] = NearestNeighbors[2] = 1;
|
---|
| 54 |
|
---|
| 55 | setPartition(radius);
|
---|
| 56 |
|
---|
| 57 | // allocate linked cell structure
|
---|
| 58 | AllocateCells();
|
---|
| 59 | }
|
---|
| 60 |
|
---|
| 61 | /** Constructor of LinkedCell_Model.
|
---|
| 62 | *
|
---|
| 63 | * @oaram set set of points to place into the linked cell structure
|
---|
| 64 | * @param radius desired maximum neighborhood distance
|
---|
| 65 | * @param _domain Box instance with domain size and boundary conditions
|
---|
| 66 | */
|
---|
| 67 | LinkedCell::LinkedCell_Model::LinkedCell_Model(IPointCloud &set, const double radius, const Box &_domain) :
|
---|
| 68 | internal_Sizes(NULL),
|
---|
| 69 | domain(_domain)
|
---|
| 70 | {
|
---|
| 71 | Info info(__func__);
|
---|
| 72 |
|
---|
| 73 | // get the partition of the domain
|
---|
| 74 | setPartition(radius);
|
---|
| 75 |
|
---|
| 76 | // allocate linked cell structure
|
---|
| 77 | AllocateCells();
|
---|
| 78 |
|
---|
| 79 | insertPointCloud(set);
|
---|
| 80 | }
|
---|
| 81 |
|
---|
| 82 | /** Destructor of class LinkedCell_Model.
|
---|
| 83 | *
|
---|
| 84 | */
|
---|
| 85 | LinkedCell::LinkedCell_Model::~LinkedCell_Model()
|
---|
| 86 | {
|
---|
| 87 | Reset();
|
---|
| 88 | }
|
---|
| 89 |
|
---|
| 90 |
|
---|
| 91 | /** Allocates as much cells per axis as required by
|
---|
| 92 | * LinkedCell_Model::BoxPartition.
|
---|
| 93 | *
|
---|
| 94 | */
|
---|
| 95 | void LinkedCell::LinkedCell_Model::AllocateCells()
|
---|
| 96 | {
|
---|
| 97 | // resize array
|
---|
| 98 | tripleIndex index;
|
---|
| 99 | for (int i=0;i<NDIM;i++)
|
---|
| 100 | index[i] = static_cast<LinkedCellArray::index>(Dimensions.at(i,i));
|
---|
| 101 | N.resize(index);
|
---|
| 102 | ASSERT(getSize(0)*getSize(1)*getSize(2) < MAX_LINKEDCELLNODES,
|
---|
| 103 | "LinkedCell_Model::AllocateCells() - Number linked of linked cell nodes exceded hard-coded limit, use greater edge length!");
|
---|
| 104 |
|
---|
| 105 | // allocate LinkedCell instances
|
---|
| 106 | for(index[0] = 0; index[0] != static_cast<LinkedCellArray::index>(Dimensions.at(0,0)); ++index[0]) {
|
---|
| 107 | for(index[1] = 0; index[1] != static_cast<LinkedCellArray::index>(Dimensions.at(1,1)); ++index[1]) {
|
---|
| 108 | for(index[2] = 0; index[2] != static_cast<LinkedCellArray::index>(Dimensions.at(2,2)); ++index[2]) {
|
---|
| 109 | LOG(5, "INFO: Creating cell at " << index[0] << " " << index[1] << " " << index[2] << ".");
|
---|
| 110 | N(index) = new LinkedCell(index);
|
---|
| 111 | }
|
---|
| 112 | }
|
---|
| 113 | }
|
---|
| 114 | }
|
---|
| 115 |
|
---|
| 116 | /** Frees all Linked Cell instances and sets array dimensions to (0,0,0).
|
---|
| 117 | *
|
---|
| 118 | */
|
---|
| 119 | void LinkedCell::LinkedCell_Model::Reset()
|
---|
| 120 | {
|
---|
| 121 | // free all LinkedCell instances
|
---|
| 122 | for(iterator3 iter3 = N.begin(); iter3 != N.end(); ++iter3) {
|
---|
| 123 | for(iterator2 iter2 = (*iter3).begin(); iter2 != (*iter3).end(); ++iter2) {
|
---|
| 124 | for(iterator1 iter1 = (*iter2).begin(); iter1 != (*iter2).end(); ++iter1) {
|
---|
| 125 | delete *iter1;
|
---|
| 126 | }
|
---|
| 127 | }
|
---|
| 128 | }
|
---|
| 129 | // set dimensions to zero
|
---|
| 130 | N.resize(boost::extents[0][0][0]);
|
---|
| 131 | }
|
---|
| 132 |
|
---|
| 133 | /** Inserts all points contained in \a set.
|
---|
| 134 | *
|
---|
| 135 | * @param set set with points to insert into linked cell structure
|
---|
| 136 | */
|
---|
| 137 | void LinkedCell::LinkedCell_Model::insertPointCloud(IPointCloud &set)
|
---|
| 138 | {
|
---|
| 139 | if (set.IsEmpty()) {
|
---|
| 140 | ELOG(1, "set is NULL or contains no linked cell nodes!");
|
---|
| 141 | return;
|
---|
| 142 | }
|
---|
| 143 |
|
---|
| 144 | // put each atom into its respective cell
|
---|
| 145 | set.GoToFirst();
|
---|
| 146 | while (!set.IsEnd()) {
|
---|
| 147 | TesselPoint *Walker = set.GetPoint();
|
---|
| 148 | addNode(Walker);
|
---|
| 149 | set.GoToNext();
|
---|
| 150 | }
|
---|
| 151 | }
|
---|
| 152 |
|
---|
| 153 | /** Calculates the required edge length for the given desired distance.
|
---|
| 154 | *
|
---|
| 155 | * We need to make some matrix transformations in order to obtain the required
|
---|
| 156 | * edge lengths per axis. Goal is guarantee that whatever the shape of the
|
---|
| 157 | * domain that always all points at least up to \a distance away are contained
|
---|
| 158 | * in the nearest neighboring cells.
|
---|
| 159 | *
|
---|
| 160 | * @param distance distance of this linked cell array
|
---|
| 161 | */
|
---|
| 162 | void LinkedCell::LinkedCell_Model::setPartition(double distance)
|
---|
| 163 | {
|
---|
| 164 | // generate box matrix of desired edge length
|
---|
| 165 | RealSpaceMatrix neighborhood;
|
---|
| 166 | neighborhood.setIdentity();
|
---|
| 167 | neighborhood *= distance;
|
---|
| 168 |
|
---|
| 169 | // obtain refs to both domain matrix transformations
|
---|
| 170 | //const RealSpaceMatrix &M = domain.getM();
|
---|
| 171 | const RealSpaceMatrix &Minv = domain.getMinv();
|
---|
| 172 |
|
---|
| 173 | RealSpaceMatrix output = Minv * neighborhood;
|
---|
| 174 |
|
---|
| 175 | std::cout << Minv << " * " << neighborhood << " = " << output << std::endl;
|
---|
| 176 |
|
---|
| 177 | Dimensions = output.invert();
|
---|
| 178 | Partition = Minv*Dimensions; //
|
---|
| 179 |
|
---|
| 180 | std::cout << "Dimensions are then " << Dimensions << std::endl;
|
---|
| 181 | std::cout << "Partition matrix is then " << Partition << std::endl;
|
---|
| 182 | }
|
---|
| 183 |
|
---|
| 184 | /** Returns the number of required neighbor-shells to get all neighboring points
|
---|
| 185 | * in the given \a distance.
|
---|
| 186 | *
|
---|
| 187 | * @param distance radius of neighborhood sphere
|
---|
| 188 | * @return number of LinkedCell's per dimension to get all neighbors
|
---|
| 189 | */
|
---|
| 190 | const LinkedCell::tripleIndex LinkedCell::LinkedCell_Model::getStep(const double distance) const
|
---|
| 191 | {
|
---|
| 192 | tripleIndex index;
|
---|
| 193 | index[0] = index[1] = index[2] = 0;
|
---|
| 194 |
|
---|
| 195 | if (fabs(distance) < std::numeric_limits<double>::min())
|
---|
| 196 | return index;
|
---|
| 197 | // generate box matrix of desired edge length
|
---|
| 198 | RealSpaceMatrix neighborhood;
|
---|
| 199 | neighborhood.setIdentity();
|
---|
| 200 | neighborhood *= distance;
|
---|
| 201 |
|
---|
| 202 | const RealSpaceMatrix output = Partition * neighborhood;
|
---|
| 203 |
|
---|
| 204 | //std::cout << "GetSteps: " << Partition << " * " << neighborhood << " = " << output << std::endl;
|
---|
| 205 |
|
---|
| 206 | const RealSpaceMatrix steps = output;
|
---|
| 207 | for (size_t i =0; i<NDIM; ++i)
|
---|
| 208 | index[i] = ceil(steps.at(i,i));
|
---|
| 209 | LOG(2, "INFO: number of shells are ("+toString(index[0])+","+toString(index[1])+","+toString(index[2])+").");
|
---|
| 210 |
|
---|
| 211 | return index;
|
---|
| 212 | }
|
---|
| 213 |
|
---|
| 214 | /** Calculates the index of the cell \a position would belong to.
|
---|
| 215 | *
|
---|
| 216 | * @param position position whose associated cell to calculate
|
---|
| 217 | * @return index of the cell
|
---|
| 218 | */
|
---|
| 219 | const LinkedCell::tripleIndex LinkedCell::LinkedCell_Model::getIndexToVector(const Vector &position) const
|
---|
| 220 | {
|
---|
| 221 | tripleIndex index;
|
---|
| 222 | Vector x(Partition*position);
|
---|
| 223 | LOG(2, "INFO: Transformed position is " << x << ".");
|
---|
| 224 | for (int i=0;i<NDIM;i++) {
|
---|
| 225 | index[i] = static_cast<LinkedCellArray::index>(floor(x[i]));
|
---|
| 226 | }
|
---|
| 227 | return index;
|
---|
| 228 | }
|
---|
| 229 |
|
---|
| 230 | /** Adds a node to the linked cell structure
|
---|
| 231 | *
|
---|
| 232 | * @param Walker node to add
|
---|
| 233 | */
|
---|
| 234 | void LinkedCell::LinkedCell_Model::addNode(TesselPoint *&Walker)
|
---|
| 235 | {
|
---|
| 236 | tripleIndex index = getIndexToVector(Walker->getPosition());
|
---|
| 237 | LOG(2, "INFO: " << *Walker << " goes into cell " << index[0] << ", " << index[1] << ", " << index[2] << ".");
|
---|
| 238 | LOG(2, "INFO: Cell's indices are "
|
---|
| 239 | << N(index)->getIndex(0) << " "
|
---|
| 240 | << N(index)->getIndex(1) << " "
|
---|
| 241 | << N(index)->getIndex(2) << ".");
|
---|
| 242 | N(index)->addPoint(Walker);
|
---|
| 243 | std::pair<MapPointToCell::iterator, bool> inserter = CellLookup.insert( std::make_pair(Walker, N(index)) );
|
---|
| 244 | ASSERT( inserter.second,
|
---|
| 245 | "LinkedCell_Model::addNode() - Walker "
|
---|
| 246 | +toString(*Walker)+" is already present in cell "
|
---|
| 247 | +toString((inserter.first)->second->getIndex(0))+" "
|
---|
| 248 | +toString((inserter.first)->second->getIndex(1))+" "
|
---|
| 249 | +toString((inserter.first)->second->getIndex(2))+".");
|
---|
| 250 | }
|
---|
| 251 |
|
---|
| 252 | /** Removes a node to the linked cell structure
|
---|
| 253 | *
|
---|
| 254 | * We do nothing of Walker is not found
|
---|
| 255 | *
|
---|
| 256 | * @param Walker node to remove
|
---|
| 257 | */
|
---|
| 258 | void LinkedCell::LinkedCell_Model::deleteNode(const TesselPoint *Walker)
|
---|
| 259 | {
|
---|
| 260 | MapPointToCell::iterator iter = CellLookup.begin();
|
---|
| 261 | for (; iter != CellLookup.end(); ++iter)
|
---|
| 262 | if (iter->first == Walker)
|
---|
| 263 | break;
|
---|
| 264 | ASSERT(iter != CellLookup.end(),
|
---|
| 265 | "LinkedCell_Model::deleteNode() - Walker not present in cell stored under CellLookup.");
|
---|
| 266 | if (iter != CellLookup.end()) {
|
---|
| 267 | CellLookup.erase(iter);
|
---|
| 268 | iter->second->deletePoint(Walker);
|
---|
| 269 | }
|
---|
| 270 | }
|
---|
| 271 |
|
---|
| 272 | /** Move Walker from current cell to another on position change.
|
---|
| 273 | *
|
---|
| 274 | * @param Walker node who has moved.
|
---|
| 275 | */
|
---|
| 276 | void LinkedCell::LinkedCell_Model::moveNode(const TesselPoint *Walker)
|
---|
| 277 | {
|
---|
| 278 | ASSERT(0, "LinkedCell_Model::moveNode() - not implemented yet.");
|
---|
| 279 | }
|
---|
| 280 |
|
---|
| 281 | /** Checks whether cell indicated by \a relative relative to LinkedCell_Model::internal_index
|
---|
| 282 | * is out of bounds.
|
---|
| 283 | *
|
---|
| 284 | * \note We do not check for boundary conditions of LinkedCell_Model::domain,
|
---|
| 285 | * we only look at the array sizes
|
---|
| 286 | *
|
---|
| 287 | * @param relative index relative to LinkedCell_Model::internal_index.
|
---|
| 288 | * @return true - relative index is still inside bounds, false - outside
|
---|
| 289 | */
|
---|
| 290 | bool LinkedCell::LinkedCell_Model::checkArrayBounds(const tripleIndex &index) const
|
---|
| 291 | {
|
---|
| 292 | bool status = true;
|
---|
| 293 | for (size_t i=0;i<NDIM;++i) {
|
---|
| 294 | status = status && (
|
---|
| 295 | (index[i] >= 0) &&
|
---|
| 296 | (index[i] < getSize(i))
|
---|
| 297 | );
|
---|
| 298 | }
|
---|
| 299 | return status;
|
---|
| 300 | }
|
---|
| 301 |
|
---|
| 302 | /** Corrects \a index according to boundary conditions of LinkedCell_Model::domain .
|
---|
| 303 | *
|
---|
| 304 | * @param index index to correct according to boundary conditions
|
---|
| 305 | */
|
---|
| 306 | void LinkedCell::LinkedCell_Model::applyBoundaryConditions(tripleIndex &index) const
|
---|
| 307 | {
|
---|
| 308 | for (size_t i=0;i<NDIM;++i) {
|
---|
| 309 | switch (domain.getConditions()[i]) {
|
---|
| 310 | case Box::Wrap:
|
---|
| 311 | if ((index[i] < 0) || (index[i] >= getSize(i)))
|
---|
| 312 | index[i] = (index[i] % getSize(i));
|
---|
| 313 | break;
|
---|
| 314 | case Box::Bounce:
|
---|
| 315 | break;
|
---|
| 316 | case Box::Ignore:
|
---|
| 317 | break;
|
---|
| 318 | default:
|
---|
| 319 | ASSERT(0, "LinkedCell_Model::checkBounds() - unknown boundary conditions.");
|
---|
| 320 | break;
|
---|
| 321 | }
|
---|
| 322 | }
|
---|
| 323 | }
|
---|
| 324 |
|
---|
| 325 | /** Calculates the interval bounds of the linked cell grid.
|
---|
| 326 | * \param index index to give relative bounds to
|
---|
| 327 | * \param step how deep to check the neighbouring cells (i.e. number of layers to check)
|
---|
| 328 | * \return pair of tripleIndex indicating lower and upper bounds
|
---|
| 329 | */
|
---|
| 330 | const LinkedCell::LinkedCell_Model::LinkedCellNeighborhoodBounds LinkedCell::LinkedCell_Model::getNeighborhoodBounds(
|
---|
| 331 | const tripleIndex &index,
|
---|
| 332 | const tripleIndex &step
|
---|
| 333 | ) const
|
---|
| 334 | {
|
---|
| 335 | LinkedCellNeighborhoodBounds neighbors;
|
---|
| 336 |
|
---|
| 337 | // calculate bounds
|
---|
| 338 | for (size_t i=0;i<NDIM;++i) {
|
---|
| 339 | switch (domain.getConditions()[i]) {
|
---|
| 340 | case Box::Wrap:
|
---|
| 341 | if (index[i] - step[i] >= 0) {
|
---|
| 342 | neighbors.first[i] = index[i] - step[i];
|
---|
| 343 | } else {
|
---|
| 344 | neighbors.first[i] = getSize(i) - (index[i] - step[i]);
|
---|
| 345 | }
|
---|
| 346 | if (index[i] + step[i] < getSize(i)) {
|
---|
| 347 | neighbors.second[i] = index[i] + step[i];
|
---|
| 348 | } else {
|
---|
| 349 | neighbors.second[i] = index[i] + step[i] - getSize(i);
|
---|
| 350 | }
|
---|
| 351 | break;
|
---|
| 352 | case Box::Bounce:
|
---|
| 353 | if (index[i] - step[i] >= 0) {
|
---|
| 354 | neighbors.first[i] = index[i] - step[i];
|
---|
| 355 | } else {
|
---|
| 356 | neighbors.first[i] = - (index[i] - step[i]);
|
---|
| 357 | }
|
---|
| 358 | if (index[i] + step[i] < getSize(i)) {
|
---|
| 359 | neighbors.second[i] = index[i] + step[i];
|
---|
| 360 | } else {
|
---|
| 361 | const size_t shift = index[i] + step[i];
|
---|
| 362 | const size_t mod = shift / getSize(i);
|
---|
| 363 | if ((mod / 2)*2 == mod) // even -> come in from lower bound
|
---|
| 364 | neighbors.second[i] = (shift % getSize(i));
|
---|
| 365 | else // odd -> come in from upper bound
|
---|
| 366 | neighbors.second[i] = getSize(i) - (shift % getSize(i));
|
---|
| 367 | }
|
---|
| 368 | break;
|
---|
| 369 | case Box::Ignore:
|
---|
| 370 | ASSERT(index[i] - step[i] >= 0,
|
---|
| 371 | "LinkedCell_Model::getNeighborhoodBounds() - lower out of bounds on 'Ignore' boundary conditions");
|
---|
| 372 | neighbors.first[i] = index[i] - step[i];
|
---|
| 373 | ASSERT (index[i] + step[i] < getSize(i),
|
---|
| 374 | "LinkedCell_Model::getNeighborhoodBounds() - upper out of bounds on 'Ignore' boundary conditions");
|
---|
| 375 | neighbors.second[i] = index[i] + step[i];
|
---|
| 376 | ASSERT(neighbors.second[i] < neighbors.first[i],
|
---|
| 377 | "LinkedCell_Model::getNeighborhoodBounds() - ["
|
---|
| 378 | +toString(neighbors.first[i])+","+toString(neighbors.second[i])
|
---|
| 379 | +"] does not make sense as array bounds.");
|
---|
| 380 | break;
|
---|
| 381 | default:
|
---|
| 382 | ASSERT(0, "LinkedCell_Model::getNeighborhoodBounds() - unknown boundary conditions.");
|
---|
| 383 | break;
|
---|
| 384 | }
|
---|
| 385 | }
|
---|
| 386 |
|
---|
| 387 | return neighbors;
|
---|
| 388 | }
|
---|