1 | /*
|
---|
2 | * Project: MoleCuilder
|
---|
3 | * Description: creates and alters molecular systems
|
---|
4 | * Copyright (C) 2012 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/bind.hpp>
|
---|
26 | #include <boost/multi_array.hpp>
|
---|
27 | #include <limits>
|
---|
28 |
|
---|
29 | #include "Atom/AtomObserver.hpp"
|
---|
30 | #include "Atom/TesselPoint.hpp"
|
---|
31 | #include "Box.hpp"
|
---|
32 | #include "CodePatterns/Assert.hpp"
|
---|
33 | #include "CodePatterns/Info.hpp"
|
---|
34 | #include "CodePatterns/Log.hpp"
|
---|
35 | #include "CodePatterns/Observer/Observer.hpp"
|
---|
36 | #include "CodePatterns/Observer/Notification.hpp"
|
---|
37 | #include "CodePatterns/toString.hpp"
|
---|
38 | #include "LinearAlgebra/RealSpaceMatrix.hpp"
|
---|
39 | #include "LinearAlgebra/Vector.hpp"
|
---|
40 | #include "LinkedCell/IPointCloud.hpp"
|
---|
41 | #include "LinkedCell/LinkedCell.hpp"
|
---|
42 | #include "LinkedCell/LinkedCell_Model_changeModel.hpp"
|
---|
43 | #include "LinkedCell/LinkedCell_Model_LinkedCellArrayCache.hpp"
|
---|
44 | #include "World.hpp"
|
---|
45 |
|
---|
46 | // initialize static entities
|
---|
47 | LinkedCell::tripleIndex LinkedCell::LinkedCell_Model::NearestNeighbors;
|
---|
48 |
|
---|
49 |
|
---|
50 | /** Constructor of LinkedCell_Model.
|
---|
51 | *
|
---|
52 | * @param radius desired maximum neighborhood distance
|
---|
53 | * @param _domain Box instance with domain size and boundary conditions
|
---|
54 | */
|
---|
55 | LinkedCell::LinkedCell_Model::LinkedCell_Model(const double radius, const Box &_domain) :
|
---|
56 | ::Observer(std::string("LinkedCell_Model")+std::string("_")+toString(radius)),
|
---|
57 | Changes( new changeModel(radius) ),
|
---|
58 | internal_Sizes(NULL),
|
---|
59 | N(new LinkedCellArrayCache(Changes, boost::bind(&changeModel::performUpdates, Changes), std::string("N_cached"))),
|
---|
60 | domain(_domain)
|
---|
61 | {
|
---|
62 | // set default argument
|
---|
63 | NearestNeighbors[0] = NearestNeighbors[1] = NearestNeighbors[2] = 1;
|
---|
64 |
|
---|
65 | // get the partition of the domain
|
---|
66 | setPartition(radius);
|
---|
67 |
|
---|
68 | // allocate linked cell structure
|
---|
69 | AllocateCells();
|
---|
70 |
|
---|
71 | // sign in to AtomObserver
|
---|
72 | startListening();
|
---|
73 | }
|
---|
74 |
|
---|
75 | /** Constructor of LinkedCell_Model.
|
---|
76 | *
|
---|
77 | * @oaram set set of points to place into the linked cell structure
|
---|
78 | * @param radius desired maximum neighborhood distance
|
---|
79 | * @param _domain Box instance with domain size and boundary conditions
|
---|
80 | */
|
---|
81 | LinkedCell::LinkedCell_Model::LinkedCell_Model(IPointCloud &set, const double radius, const Box &_domain) :
|
---|
82 | ::Observer(std::string("LinkedCell_Model")+std::string("_")+toString(radius)),
|
---|
83 | Changes( new changeModel(radius) ),
|
---|
84 | internal_Sizes(NULL),
|
---|
85 | N(new LinkedCellArrayCache(Changes, boost::bind(&changeModel::performUpdates, Changes), std::string("N_cached"))),
|
---|
86 | domain(_domain)
|
---|
87 | {
|
---|
88 | Info info(__func__);
|
---|
89 |
|
---|
90 | // get the partition of the domain
|
---|
91 | setPartition(radius);
|
---|
92 |
|
---|
93 | // allocate linked cell structure
|
---|
94 | AllocateCells();
|
---|
95 |
|
---|
96 | insertPointCloud(set);
|
---|
97 |
|
---|
98 | // sign in to AtomObserver
|
---|
99 | startListening();
|
---|
100 | }
|
---|
101 |
|
---|
102 | /** Destructor of class LinkedCell_Model.
|
---|
103 | *
|
---|
104 | */
|
---|
105 | LinkedCell::LinkedCell_Model::~LinkedCell_Model()
|
---|
106 | {
|
---|
107 | // sign off from observables
|
---|
108 | stopListening();
|
---|
109 |
|
---|
110 | // reset linked cell structure
|
---|
111 | Reset();
|
---|
112 | delete N;
|
---|
113 |
|
---|
114 | // delete change queue
|
---|
115 | delete Changes;
|
---|
116 | }
|
---|
117 |
|
---|
118 | /** Signs in to AtomObserver and World to known about all changes.
|
---|
119 | *
|
---|
120 | */
|
---|
121 | void LinkedCell::LinkedCell_Model::startListening()
|
---|
122 | {
|
---|
123 | World::getInstance().signOn(this, World::AtomInserted);
|
---|
124 | World::getInstance().signOn(this, World::AtomRemoved);
|
---|
125 | AtomObserver::getInstance().signOn(this, AtomObservable::PositionChanged);
|
---|
126 | }
|
---|
127 |
|
---|
128 | /** Signs off from AtomObserver and World.
|
---|
129 | *
|
---|
130 | */
|
---|
131 | void LinkedCell::LinkedCell_Model::stopListening()
|
---|
132 | {
|
---|
133 | World::getInstance().signOff(this, World::AtomInserted);
|
---|
134 | World::getInstance().signOff(this, World::AtomRemoved);
|
---|
135 | AtomObserver::getInstance().signOff(this, AtomObservable::PositionChanged);
|
---|
136 | }
|
---|
137 |
|
---|
138 |
|
---|
139 | /** Allocates as much cells per axis as required by
|
---|
140 | * LinkedCell_Model::BoxPartition.
|
---|
141 | *
|
---|
142 | */
|
---|
143 | void LinkedCell::LinkedCell_Model::AllocateCells()
|
---|
144 | {
|
---|
145 | // resize array
|
---|
146 | tripleIndex index;
|
---|
147 | for (int i=0;i<NDIM;i++)
|
---|
148 | index[i] = static_cast<LinkedCellArray::index>(Dimensions.at(i,i));
|
---|
149 | N->setN().resize(index);
|
---|
150 | ASSERT(getSize(0)*getSize(1)*getSize(2) < MAX_LINKEDCELLNODES,
|
---|
151 | "LinkedCell_Model::AllocateCells() - Number linked of linked cell nodes exceeded hard-coded limit, use greater edge length!");
|
---|
152 | LOG(1, "INFO: Allocating array ("
|
---|
153 | << getSize(0) << ","
|
---|
154 | << getSize(1) << ","
|
---|
155 | << getSize(2) << ") for a new LinkedCell_Model.");
|
---|
156 |
|
---|
157 | // allocate LinkedCell instances
|
---|
158 | for(index[0] = 0; index[0] != static_cast<LinkedCellArray::index>(Dimensions.at(0,0)); ++index[0]) {
|
---|
159 | for(index[1] = 0; index[1] != static_cast<LinkedCellArray::index>(Dimensions.at(1,1)); ++index[1]) {
|
---|
160 | for(index[2] = 0; index[2] != static_cast<LinkedCellArray::index>(Dimensions.at(2,2)); ++index[2]) {
|
---|
161 | LOG(5, "INFO: Creating cell at " << index[0] << " " << index[1] << " " << index[2] << ".");
|
---|
162 | N->setN()(index) = new LinkedCell(index);
|
---|
163 | }
|
---|
164 | }
|
---|
165 | }
|
---|
166 | }
|
---|
167 |
|
---|
168 | /** Frees all Linked Cell instances and sets array dimensions to (0,0,0).
|
---|
169 | *
|
---|
170 | */
|
---|
171 | void LinkedCell::LinkedCell_Model::Reset()
|
---|
172 | {
|
---|
173 | // free all LinkedCell instances
|
---|
174 | for(iterator3 iter3 = N->setN().begin(); iter3 != N->setN().end(); ++iter3) {
|
---|
175 | for(iterator2 iter2 = (*iter3).begin(); iter2 != (*iter3).end(); ++iter2) {
|
---|
176 | for(iterator1 iter1 = (*iter2).begin(); iter1 != (*iter2).end(); ++iter1) {
|
---|
177 | delete *iter1;
|
---|
178 | }
|
---|
179 | }
|
---|
180 | }
|
---|
181 | // set dimensions to zero
|
---|
182 | N->setN().resize(boost::extents[0][0][0]);
|
---|
183 | }
|
---|
184 |
|
---|
185 | /** Inserts all points contained in \a set.
|
---|
186 | *
|
---|
187 | * @param set set with points to insert into linked cell structure
|
---|
188 | */
|
---|
189 | void LinkedCell::LinkedCell_Model::insertPointCloud(IPointCloud &set)
|
---|
190 | {
|
---|
191 | if (set.IsEmpty()) {
|
---|
192 | ELOG(1, "set is NULL or contains no linked cell nodes!");
|
---|
193 | return;
|
---|
194 | }
|
---|
195 |
|
---|
196 | // put each atom into its respective cell
|
---|
197 | set.GoToFirst();
|
---|
198 | while (!set.IsEnd()) {
|
---|
199 | TesselPoint *Walker = set.GetPoint();
|
---|
200 | addNode(Walker);
|
---|
201 | set.GoToNext();
|
---|
202 | }
|
---|
203 | }
|
---|
204 |
|
---|
205 | /** Calculates the required edge length for the given desired distance.
|
---|
206 | *
|
---|
207 | * We need to make some matrix transformations in order to obtain the required
|
---|
208 | * edge lengths per axis. Goal is guarantee that whatever the shape of the
|
---|
209 | * domain that always all points at least up to \a distance away are contained
|
---|
210 | * in the nearest neighboring cells.
|
---|
211 | *
|
---|
212 | * @param distance distance of this linked cell array
|
---|
213 | */
|
---|
214 | void LinkedCell::LinkedCell_Model::setPartition(double distance)
|
---|
215 | {
|
---|
216 | // generate box matrix of desired edge length
|
---|
217 | RealSpaceMatrix neighborhood;
|
---|
218 | neighborhood.setIdentity();
|
---|
219 | neighborhood *= distance;
|
---|
220 |
|
---|
221 | // obtain refs to both domain matrix transformations
|
---|
222 | //const RealSpaceMatrix &M = domain.getM();
|
---|
223 | const RealSpaceMatrix &Minv = domain.getMinv();
|
---|
224 |
|
---|
225 | RealSpaceMatrix output = Minv * neighborhood;
|
---|
226 | std::cout << Minv << " * " << neighborhood << " = " << output << std::endl;
|
---|
227 |
|
---|
228 | Dimensions = output.invert();
|
---|
229 | std::cout << "Dimensions are then " << Dimensions << std::endl;
|
---|
230 |
|
---|
231 | // now dimensions is floating-point, but we need it to be integer (for allocation)
|
---|
232 | for (size_t col = 0; col < NDIM; ++col) {
|
---|
233 | for (size_t row = 0; row < NDIM; ++row) {
|
---|
234 | ASSERT(fabs(Dimensions.at(row,col) - Dimensions.at(col,row)) < 1.e+3*std::numeric_limits<double>::epsilon(),
|
---|
235 | "LinkedCell_Model::setPartition() - Dimensions is not symmetric by "
|
---|
236 | +toString(fabs(Dimensions.at(row,col) - Dimensions.at(col,row)))+ ".");
|
---|
237 | if (col != row) {
|
---|
238 | ASSERT(fabs(Dimensions.at(row,col)) < 1.e+3*std::numeric_limits<double>::epsilon(),
|
---|
239 | "LinkedCell_Model::setPartition() - Dimensions is not diagonal by "
|
---|
240 | +toString(fabs(Dimensions.at(row,col)))+".");
|
---|
241 | } else {
|
---|
242 | Dimensions.set(row,col, floor(Dimensions.at(row,col)));
|
---|
243 | }
|
---|
244 | }
|
---|
245 | }
|
---|
246 |
|
---|
247 |
|
---|
248 | Partition = Minv*Dimensions; //
|
---|
249 |
|
---|
250 | std::cout << "Partition matrix is then " << Partition << std::endl;
|
---|
251 | }
|
---|
252 |
|
---|
253 | /** Returns the number of required neighbor-shells to get all neighboring points
|
---|
254 | * in the given \a distance.
|
---|
255 | *
|
---|
256 | * @param distance radius of neighborhood sphere
|
---|
257 | * @return number of LinkedCell's per dimension to get all neighbors
|
---|
258 | */
|
---|
259 | const LinkedCell::tripleIndex LinkedCell::LinkedCell_Model::getStep(const double distance) const
|
---|
260 | {
|
---|
261 | tripleIndex index;
|
---|
262 | index[0] = index[1] = index[2] = 0;
|
---|
263 |
|
---|
264 | if (fabs(distance) < std::numeric_limits<double>::min())
|
---|
265 | return index;
|
---|
266 | // generate box matrix of desired edge length
|
---|
267 | RealSpaceMatrix neighborhood;
|
---|
268 | neighborhood.setIdentity();
|
---|
269 | neighborhood *= distance;
|
---|
270 |
|
---|
271 | const RealSpaceMatrix output = Partition * neighborhood;
|
---|
272 |
|
---|
273 | //std::cout << "GetSteps: " << Partition << " * " << neighborhood << " = " << output << std::endl;
|
---|
274 |
|
---|
275 | const RealSpaceMatrix steps = output;
|
---|
276 | for (size_t i =0; i<NDIM; ++i)
|
---|
277 | index[i] = ceil(steps.at(i,i));
|
---|
278 | LOG(2, "INFO: number of shells are ("+toString(index[0])+","+toString(index[1])+","+toString(index[2])+").");
|
---|
279 |
|
---|
280 | return index;
|
---|
281 | }
|
---|
282 |
|
---|
283 | /** Calculates the index of the cell \a position would belong to.
|
---|
284 | *
|
---|
285 | * @param position position whose associated cell to calculate
|
---|
286 | * @return index of the cell
|
---|
287 | */
|
---|
288 | const LinkedCell::tripleIndex LinkedCell::LinkedCell_Model::getIndexToVector(const Vector &position) const
|
---|
289 | {
|
---|
290 | tripleIndex index;
|
---|
291 | Vector x(Partition*position);
|
---|
292 | LOG(2, "INFO: Transformed position is " << x << ".");
|
---|
293 | for (int i=0;i<NDIM;i++) {
|
---|
294 | index[i] = static_cast<LinkedCellArray::index>(floor(x[i]));
|
---|
295 | }
|
---|
296 | return index;
|
---|
297 | }
|
---|
298 |
|
---|
299 | /** Adds an update to the list of lazy changes to add a node.
|
---|
300 | *
|
---|
301 | * @param Walker node to add
|
---|
302 | */
|
---|
303 | void LinkedCell::LinkedCell_Model::addNode(const TesselPoint *Walker)
|
---|
304 | {
|
---|
305 | LOG(2, "INFO: Requesting update to add node " << *Walker << ".");
|
---|
306 | Changes->addUpdate(Walker, 0, boost::bind(&LinkedCell_Model::addNode_internal, this, _1));
|
---|
307 | }
|
---|
308 |
|
---|
309 | /** Adds an update to the list of lazy changes to add remove a node.
|
---|
310 | *
|
---|
311 | * We do nothing of Walker is not found
|
---|
312 | *
|
---|
313 | * @param Walker node to remove
|
---|
314 | */
|
---|
315 | void LinkedCell::LinkedCell_Model::deleteNode(const TesselPoint *Walker)
|
---|
316 | {
|
---|
317 | LOG(2, "INFO: Requesting update to delete node " << *Walker << ".");
|
---|
318 | Changes->addUpdate(Walker, 0, boost::bind(&LinkedCell_Model::deleteNode_internal, this, _1));
|
---|
319 | }
|
---|
320 |
|
---|
321 | /** Adds an update to the list of lazy changes to move a node.
|
---|
322 | *
|
---|
323 | * @param Walker node who has moved.
|
---|
324 | */
|
---|
325 | void LinkedCell::LinkedCell_Model::moveNode(const TesselPoint *Walker)
|
---|
326 | {
|
---|
327 | LOG(2, "INFO: Requesting update to move node " << *Walker << " to position "
|
---|
328 | << Walker->getPosition() << ".");
|
---|
329 | Changes->addUpdate(Walker, 10, boost::bind(&LinkedCell_Model::moveNode_internal, this, _1));
|
---|
330 | }
|
---|
331 |
|
---|
332 | /** Internal function to add a node to the linked cell structure
|
---|
333 | *
|
---|
334 | * @param Walker node to add
|
---|
335 | */
|
---|
336 | void LinkedCell::LinkedCell_Model::addNode_internal(const TesselPoint *Walker)
|
---|
337 | {
|
---|
338 | // find index
|
---|
339 | tripleIndex index = getIndexToVector(Walker->getPosition());
|
---|
340 | LOG(2, "INFO: " << *Walker << " goes into cell " << index[0] << ", " << index[1] << ", " << index[2] << ".");
|
---|
341 | LOG(2, "INFO: Cell's indices are "
|
---|
342 | << (N->getN())(index)->getIndex(0) << " "
|
---|
343 | << (N->getN())(index)->getIndex(1) << " "
|
---|
344 | << (N->getN())(index)->getIndex(2) << ".");
|
---|
345 | // add to cell
|
---|
346 | (N->setN())(index)->addPoint(Walker);
|
---|
347 | // add to index with check for presence
|
---|
348 | std::pair<MapPointToCell::iterator, bool> inserter = CellLookup.insert( std::make_pair(Walker, (N->setN())(index)) );
|
---|
349 | ASSERT( inserter.second,
|
---|
350 | "LinkedCell_Model::addNode() - Walker "
|
---|
351 | +toString(*Walker)+" is already present in cell "
|
---|
352 | +toString((inserter.first)->second->getIndex(0))+" "
|
---|
353 | +toString((inserter.first)->second->getIndex(1))+" "
|
---|
354 | +toString((inserter.first)->second->getIndex(2))+".");
|
---|
355 | }
|
---|
356 |
|
---|
357 | /** Internal function to remove a node to the linked cell structure
|
---|
358 | *
|
---|
359 | * We do nothing of Walker is not found
|
---|
360 | *
|
---|
361 | * @param Walker node to remove
|
---|
362 | */
|
---|
363 | void LinkedCell::LinkedCell_Model::deleteNode_internal(const TesselPoint *Walker)
|
---|
364 | {
|
---|
365 | MapPointToCell::iterator iter = CellLookup.find(Walker);
|
---|
366 | ASSERT(iter != CellLookup.end(),
|
---|
367 | "LinkedCell_Model::deleteNode() - Walker not present in cell stored under CellLookup.");
|
---|
368 | if (iter != CellLookup.end()) {
|
---|
369 | // remove from lookup
|
---|
370 | CellLookup.erase(iter);
|
---|
371 | // remove from cell
|
---|
372 | iter->second->deletePoint(Walker);
|
---|
373 | }
|
---|
374 | }
|
---|
375 |
|
---|
376 | /** Internal function to move node from current cell to another on position change.
|
---|
377 | *
|
---|
378 | * @param Walker node who has moved.
|
---|
379 | */
|
---|
380 | void LinkedCell::LinkedCell_Model::moveNode_internal(const TesselPoint *Walker)
|
---|
381 | {
|
---|
382 | MapPointToCell::iterator iter = CellLookup.find(Walker);
|
---|
383 | ASSERT(iter != CellLookup.end(),
|
---|
384 | "LinkedCell_Model::deleteNode() - Walker not present in cell stored under CellLookup.");
|
---|
385 | if (iter != CellLookup.end()) {
|
---|
386 | tripleIndex index = getIndexToVector(Walker->getPosition());
|
---|
387 | if (index != iter->second->getIndices()) {
|
---|
388 | // remove in old cell
|
---|
389 | iter->second->deletePoint(Walker);
|
---|
390 | // add to new cell
|
---|
391 | N->setN()(index)->addPoint(Walker);
|
---|
392 | // update lookup
|
---|
393 | iter->second = N->setN()(index);
|
---|
394 | }
|
---|
395 | }
|
---|
396 | }
|
---|
397 |
|
---|
398 | /** Checks whether cell indicated by \a relative relative to LinkedCell_Model::internal_index
|
---|
399 | * is out of bounds.
|
---|
400 | *
|
---|
401 | * \note We do not check for boundary conditions of LinkedCell_Model::domain,
|
---|
402 | * we only look at the array sizes
|
---|
403 | *
|
---|
404 | * @param relative index relative to LinkedCell_Model::internal_index.
|
---|
405 | * @return true - relative index is still inside bounds, false - outside
|
---|
406 | */
|
---|
407 | bool LinkedCell::LinkedCell_Model::checkArrayBounds(const tripleIndex &index) const
|
---|
408 | {
|
---|
409 | bool status = true;
|
---|
410 | for (size_t i=0;i<NDIM;++i) {
|
---|
411 | status = status && (
|
---|
412 | (index[i] >= 0) &&
|
---|
413 | (index[i] < getSize(i))
|
---|
414 | );
|
---|
415 | }
|
---|
416 | return status;
|
---|
417 | }
|
---|
418 |
|
---|
419 | /** Corrects \a index according to boundary conditions of LinkedCell_Model::domain .
|
---|
420 | *
|
---|
421 | * @param index index to correct according to boundary conditions
|
---|
422 | */
|
---|
423 | void LinkedCell::LinkedCell_Model::applyBoundaryConditions(tripleIndex &index) const
|
---|
424 | {
|
---|
425 | for (size_t i=0;i<NDIM;++i) {
|
---|
426 | switch (domain.getConditions()[i]) {
|
---|
427 | case Box::Wrap:
|
---|
428 | if ((index[i] < 0) || (index[i] >= getSize(i)))
|
---|
429 | index[i] = (index[i] % getSize(i));
|
---|
430 | break;
|
---|
431 | case Box::Bounce:
|
---|
432 | if (index[i] < 0)
|
---|
433 | index[i] = 0;
|
---|
434 | if (index[i] >= getSize(i))
|
---|
435 | index[i] = getSize(i)-1;
|
---|
436 | break;
|
---|
437 | case Box::Ignore:
|
---|
438 | if (index[i] < 0)
|
---|
439 | index[i] = 0;
|
---|
440 | if (index[i] >= getSize(i))
|
---|
441 | index[i] = getSize(i)-1;
|
---|
442 | break;
|
---|
443 | default:
|
---|
444 | ASSERT(0, "LinkedCell_Model::checkBounds() - unknown boundary conditions.");
|
---|
445 | break;
|
---|
446 | }
|
---|
447 | }
|
---|
448 | }
|
---|
449 |
|
---|
450 | /** Calculates the interval bounds of the linked cell grid.
|
---|
451 | *
|
---|
452 | * \note we assume for index to allows be valid, i.e. within the range of LinkedCell_Model::N.
|
---|
453 | *
|
---|
454 | * \param index index to give relative bounds to
|
---|
455 | * \param step how deep to check the neighbouring cells (i.e. number of layers to check)
|
---|
456 | * \return pair of tripleIndex indicating lower and upper bounds
|
---|
457 | */
|
---|
458 | const LinkedCell::LinkedCell_Model::LinkedCellNeighborhoodBounds LinkedCell::LinkedCell_Model::getNeighborhoodBounds(
|
---|
459 | const tripleIndex &index,
|
---|
460 | const tripleIndex &step
|
---|
461 | ) const
|
---|
462 | {
|
---|
463 | LinkedCellNeighborhoodBounds neighbors;
|
---|
464 |
|
---|
465 | // calculate bounds
|
---|
466 | for (size_t i=0;i<NDIM;++i) {
|
---|
467 | ASSERT(index[i] >= 0,
|
---|
468 | "LinkedCell_Model::getNeighborhoodBounds() - index "+toString(index)+" out of lower bounds.");
|
---|
469 | ASSERT (index[i] < getSize(i),
|
---|
470 | "LinkedCell_Model::getNeighborhoodBounds() - index "+toString(index)+" out of upper bounds.");
|
---|
471 | switch (domain.getConditions()[i]) {
|
---|
472 | case Box::Wrap:
|
---|
473 | if ((index[i] - step[i]) < 0)
|
---|
474 | neighbors.first[i] = getSize(i) + (index[i] - step[i]);
|
---|
475 | else if ((index[i] - step[i]) >= getSize(i))
|
---|
476 | neighbors.first[i] = (index[i] - step[i]) - getSize(i);
|
---|
477 | else
|
---|
478 | neighbors.first[i] = index[i] - step[i];
|
---|
479 | neighbors.second[i] = 2*step[i]+1;
|
---|
480 | break;
|
---|
481 | case Box::Bounce:
|
---|
482 | neighbors.second[i] = 2*step[i]+1;
|
---|
483 | if (index[i] - step[i] >= 0) {
|
---|
484 | neighbors.first[i] = index[i] - step[i];
|
---|
485 | } else {
|
---|
486 | neighbors.first[i] = 0;
|
---|
487 | neighbors.second[i] = index[i] + step[i]+1;
|
---|
488 | }
|
---|
489 | if (index[i] + step[i] >= getSize(i)) {
|
---|
490 | neighbors.second[i] = getSize(i) - (index[i] - step[i]);
|
---|
491 | }
|
---|
492 | break;
|
---|
493 | case Box::Ignore:
|
---|
494 | if (index[i] - step[i] < 0)
|
---|
495 | neighbors.first[i] = 0;
|
---|
496 | else
|
---|
497 | neighbors.first[i] = index[i] - step[i];
|
---|
498 | if ((neighbors.first[i] + 2*step[i]+1) >= getSize(i))
|
---|
499 | neighbors.second[i] = getSize(i) - neighbors.first[i];
|
---|
500 | else
|
---|
501 | neighbors.second[i] = 2*step[i]+1;
|
---|
502 | break;
|
---|
503 | default:
|
---|
504 | ASSERT(0, "LinkedCell_Model::getNeighborhoodBounds() - unknown boundary conditions.");
|
---|
505 | break;
|
---|
506 | }
|
---|
507 | }
|
---|
508 | LOG(3, "INFO: Resulting neighborhood bounds are [" << neighbors.first << "]<->[" << neighbors.second << "].");
|
---|
509 |
|
---|
510 | return neighbors;
|
---|
511 | }
|
---|
512 |
|
---|
513 | /** Returns a reference to the cell indicated by LinkedCell_Model::internal_index.
|
---|
514 | *
|
---|
515 | * \return LinkedCell ref to current cell
|
---|
516 | */
|
---|
517 | const LinkedCell::LinkedCell& LinkedCell::LinkedCell_Model::getCell(const tripleIndex &index) const
|
---|
518 | {
|
---|
519 | return *(N->getN()(index));
|
---|
520 | }
|
---|
521 |
|
---|
522 |
|
---|
523 | /** Returns size of array for given \a dim.
|
---|
524 | *
|
---|
525 | * @param dim desired dimension
|
---|
526 | * @return size of array along dimension
|
---|
527 | */
|
---|
528 | LinkedCell::LinkedCellArray::index LinkedCell::LinkedCell_Model::getSize(const size_t dim) const
|
---|
529 | {
|
---|
530 | ASSERT((dim >= 0) && (dim < NDIM),
|
---|
531 | "LinkedCell_Model::getSize() - dimension "
|
---|
532 | +toString(dim)+" is out of bounds.");
|
---|
533 | return N->getN().shape()[dim];
|
---|
534 | }
|
---|
535 |
|
---|
536 | /** Callback function for Observer mechanism.
|
---|
537 | *
|
---|
538 | * @param publisher reference to the Observable that calls
|
---|
539 | */
|
---|
540 | void LinkedCell::LinkedCell_Model::update(Observable *publisher)
|
---|
541 | {
|
---|
542 | ELOG(2, "LinkedCell_Model received inconclusive general update from "
|
---|
543 | << publisher << ".");
|
---|
544 | }
|
---|
545 |
|
---|
546 | /** Callback function for the Notifications mechanism.
|
---|
547 | *
|
---|
548 | * @param publisher reference to the Observable that calls
|
---|
549 | * @param notification specific notification as cause of the call
|
---|
550 | */
|
---|
551 | void LinkedCell::LinkedCell_Model::recieveNotification(Observable *publisher, Notification_ptr notification)
|
---|
552 | {
|
---|
553 | // either it's the World or from the atom (through relay) itself
|
---|
554 | if (publisher == World::getPointer()) {
|
---|
555 | switch(notification->getChannelNo()) {
|
---|
556 | case World::AtomInserted:
|
---|
557 | addNode(World::getInstance().lastChanged<atom>());
|
---|
558 | break;
|
---|
559 | case World::AtomRemoved:
|
---|
560 | deleteNode(World::getInstance().lastChanged<atom>());
|
---|
561 | break;
|
---|
562 | }
|
---|
563 | } else {
|
---|
564 | switch(notification->getChannelNo()) {
|
---|
565 | case AtomObservable::PositionChanged:
|
---|
566 | {
|
---|
567 | moveNode(dynamic_cast<const TesselPoint *>(publisher));
|
---|
568 | break;
|
---|
569 | }
|
---|
570 | default:
|
---|
571 | LOG(2, "LinkedCell_Model received unwanted notification from AtomObserver's channel "
|
---|
572 | << notification->getChannelNo() << ".");
|
---|
573 | break;
|
---|
574 | }
|
---|
575 | }
|
---|
576 | }
|
---|
577 |
|
---|
578 | /** Callback function when an Observer dies.
|
---|
579 | *
|
---|
580 | * @param publisher reference to the Observable that calls
|
---|
581 | */
|
---|
582 | void LinkedCell::LinkedCell_Model::subjectKilled(Observable *publisher)
|
---|
583 | {}
|
---|
584 |
|
---|