14 #ifndef ESYS_LSMCIRCULARNEIGHBOURTABLE_HPP
15 #define ESYS_LSMCIRCULARNEIGHBOURTABLE_HPP
17 #include "Geometry/NeighbourTable.h"
18 #include <boost/pool/object_pool.hpp>
19 #include <boost/shared_ptr.hpp>
32 template <
class TmplParticle>
33 CircularNeighbourTable<TmplParticle>::CircularNeighbourTable(
34 const BoundingBox &bBox,
36 const BoolVector &periodicDimensions,
37 double circBorderWidth
39 : Inherited(bBox, gridSpacing),
40 m_periodicDimensions(periodicDimensions),
41 m_particlePoolPtr(new ParticlePool(4096)),
42 m_clonedParticleSet(),
44 m_periodicDimIndex(-1)
46 checkPeriodicDimensions();
47 if (circBorderWidth <= 0.0) {
48 circBorderWidth = this->getGridSpacing();
50 setCircularBorderWidth(circBorderWidth, this->getGridSpacing());
53 template <
class TmplParticle>
54 CircularNeighbourTable<TmplParticle>::CircularNeighbourTable(
55 const BoundingBox &bBox,
57 ParticlePoolPtr particlePoolPtr,
58 const BoolVector &periodicDimensions,
59 double circBorderWidth
61 : Inherited(bBox, gridSpacing),
62 m_periodicDimensions(periodicDimensions),
63 m_particlePoolPtr(particlePoolPtr),
64 m_clonedParticleSet(),
66 m_periodicDimIndex(-1)
68 checkPeriodicDimensions();
69 if (circBorderWidth <= 0.0) {
70 circBorderWidth = this->getGridSpacing();
72 setCircularBorderWidth(circBorderWidth, this->getGridSpacing());
75 template <
class TmplParticle>
76 void CircularNeighbourTable<TmplParticle>::checkPeriodicDimensions()
78 if (m_periodicDimensions.size() != 3) {
79 std::stringstream msg;
81 <<
"CircularNeighbourTable::CircularNeighbourTable -"
82 <<
" size of periodic dimensions argument ("
83 << m_periodicDimensions.size()
84 <<
") is not equal to 3";
85 throw std::runtime_error(msg.str());
88 for (
int i = 0; i < 3; i++) {
89 if (m_periodicDimensions[i])
91 m_periodicDimIndex = i;
96 if (numPeriodic > 1) {
97 std::stringstream msg;
99 <<
"CircularNeighbourTable::CircularNeighbourTable -"
100 <<
" only a single dimension may be periodic.";
101 throw std::runtime_error(msg.str());
105 template <
class TmplParticle>
106 CircularNeighbourTable<TmplParticle>::~CircularNeighbourTable()
110 template <
class TmplParticle>
111 void CircularNeighbourTable<TmplParticle>::setCircularBorderWidth(
112 double circBorderWidth,
116 m_circGridWidth = int(ceil(circBorderWidth/gridSpacing));
119 template <
class TmplParticle>
120 void CircularNeighbourTable<TmplParticle>::setCircularBorderWidth(
121 double circBorderWidth
124 setCircularBorderWidth(circBorderWidth, this->getGridSpacing());
127 template <
class TmplParticle>
128 void CircularNeighbourTable<TmplParticle>::resize(
129 const BoundingBox &bBox,
131 double circBorderWidth
134 if (havePeriodicDimensions())
138 (bBox.getSizes()[m_periodicDimIndex])/2.0,
142 setCircularBorderWidth(circBorderWidth, gridSpacing);
143 ParticleVector particles = getNonClonedParticles();
144 clearClonedParticles();
145 this->clearAndRecomputeGrid(bBox, gridSpacing);
147 typename ParticleVector::iterator it = particles.begin();
148 it != particles.end();
156 template <
class TmplParticle>
157 void CircularNeighbourTable<TmplParticle>::resize(
158 const BoundingBox &bBox,
162 this->resize(bBox, gridSpacing, gridSpacing);
165 template <
class TmplParticle>
166 void CircularNeighbourTable<TmplParticle>::insertClone(
168 const Vec3 &newPosition
171 Particle *pNewParticle = m_particlePoolPtr->construct(*pParticle);
172 pNewParticle->moveTo(newPosition);
173 Inherited::insert(pNewParticle);
174 m_clonedParticleSet.insert(pNewParticle);
177 template <
class TmplParticle>
178 bool CircularNeighbourTable<TmplParticle>::havePeriodicDimensions()
const
180 return (m_periodicDimIndex >= 0);
183 template <
class TmplParticle>
184 Vec3 CircularNeighbourTable<TmplParticle>::getModdedPosn(
189 havePeriodicDimensions()
192 const int &dimIdx = m_periodicDimIndex;
194 (posn[dimIdx] < this->getBBox().getMinPt()[dimIdx])
196 (posn[dimIdx] > this->getBBox().getMaxPt()[dimIdx])
199 Vec3 moddedPosn = posn;
200 const double dimSize = this->getBBox().getSizes()[dimIdx];
201 moddedPosn[dimIdx] -= this->getBBox().getMinPt()[dimIdx];
203 (moddedPosn[dimIdx] > 0.0)
206 this->getBBox().getMinPt()[dimIdx]
208 moddedPosn[dimIdx] - (floor(moddedPosn[dimIdx]/dimSize)*dimSize)
212 this->getBBox().getMaxPt()[dimIdx]
214 (fabs(moddedPosn[dimIdx]) - (floor(fabs(moddedPosn[dimIdx])/dimSize)*dimSize))
223 template <
class TmplParticle>
224 void CircularNeighbourTable<TmplParticle>::insert(Particle *pParticle)
226 pParticle->moveTo(getModdedPosn(pParticle->getPos()));
227 const Vec3L minIdx = this->getVecIndex(pParticle->getPos() - pParticle->getRad());
228 const Vec3L maxIdx = this->getVecIndex(pParticle->getPos() + pParticle->getRad());
230 this->insertInTable(pParticle, minIdx, maxIdx);
231 this->addInserted(pParticle);
233 if (havePeriodicDimensions())
235 for (
int dimIdx = 0; dimIdx < 3; dimIdx++) {
236 if (m_periodicDimensions[dimIdx]) {
237 if (minIdx[dimIdx] < (this->getMinVecIndex()[dimIdx] + m_circGridWidth)) {
238 Vec3 shift = Vec3::ZERO;
239 shift[dimIdx] = this->getBBox().getSizes()[dimIdx];
240 this->insertClone(pParticle, pParticle->getPos() + shift);
242 if (maxIdx[dimIdx] > (this->getMaxVecIndex()[dimIdx] - m_circGridWidth)) {
243 Vec3 shift = Vec3::ZERO;
244 shift[dimIdx] = this->getBBox().getSizes()[dimIdx];
245 this->insertClone(pParticle, pParticle->getPos() - shift);
252 template <
class TmplParticle>
253 void CircularNeighbourTable<TmplParticle>::insert(Particle &particle)
255 this->insert(&particle);
258 template <
class TmplParticle>
259 size_t CircularNeighbourTable<TmplParticle>::getNumClonedParticles()
const
261 return m_clonedParticleSet.size();
264 template <
class TmplParticle>
265 size_t CircularNeighbourTable<TmplParticle>::getNumParticles()
const
267 return this->size() - getNumClonedParticles();
270 template <
class TmplParticle>
271 const typename CircularNeighbourTable<TmplParticle>::BoolVector &
272 CircularNeighbourTable<TmplParticle>::getPeriodicDimensions()
const
274 return m_periodicDimensions;
277 template <
class TmplParticle>
278 bool CircularNeighbourTable<TmplParticle>::isClone(
282 return (m_clonedParticleSet.find(p) != m_clonedParticleSet.end());
285 template <
class TmplParticle>
286 typename CircularNeighbourTable<TmplParticle>::ParticleVector
287 CircularNeighbourTable<TmplParticle>::getNonClonedParticles()
289 ParticleVector all = this->getInsertedParticles();
290 ParticleVector nonCloned;
291 nonCloned.reserve(all.size()/2);
293 typename ParticleVector::iterator it = all.begin();
300 nonCloned.push_back(*it);
306 template <
class TmplParticle>
307 void CircularNeighbourTable<TmplParticle>::clearClonedParticles()
310 typename ParticleSet::iterator it = m_clonedParticleSet.begin();
311 it != m_clonedParticleSet.end();
315 m_particlePoolPtr->destroy(*it);
317 m_clonedParticleSet.clear();