ESyS-Particle  2.3.2
CircularNeighbourTable.hpp
Go to the documentation of this file.
1 
2 // //
3 // Copyright (c) 2003-2014 by The University of Queensland //
4 // Centre for Geoscience Computing //
5 // http://earth.uq.edu.au/centre-geoscience-computing //
6 // //
7 // Primary Business: Brisbane, Queensland, Australia //
8 // Licensed under the Open Software License version 3.0 //
9 // http://www.opensource.org/licenses/osl-3.0.php //
10 // //
12 
13 
14 #ifndef ESYS_LSMCIRCULARNEIGHBOURTABLE_HPP
15 #define ESYS_LSMCIRCULARNEIGHBOURTABLE_HPP
16 
18 #include <boost/pool/object_pool.hpp>
19 #include <boost/shared_ptr.hpp>
20 
21 #include <sstream>
22 #include <stdexcept>
23 #include <set>
24 
25 namespace esys
26 {
27  namespace lsm
28  {
32  template <class TmplParticle>
34  const BoundingBox &bBox,
35  double gridSpacing,
36  const BoolVector &periodicDimensions,
37  double circBorderWidth
38  )
39  : Inherited(bBox, gridSpacing),
40  m_periodicDimensions(periodicDimensions),
41  m_particlePoolPtr(new ParticlePool(4096)),
42  m_clonedParticleSet(),
43  m_circGridWidth(1),
44  m_periodicDimIndex(-1)
45  {
47  if (circBorderWidth <= 0.0) {
48  circBorderWidth = this->getGridSpacing();
49  }
50  setCircularBorderWidth(circBorderWidth, this->getGridSpacing());
51  }
52 
53  template <class TmplParticle>
55  const BoundingBox &bBox,
56  double gridSpacing,
57  ParticlePoolPtr particlePoolPtr,
58  const BoolVector &periodicDimensions,
59  double circBorderWidth
60  )
61  : Inherited(bBox, gridSpacing),
62  m_periodicDimensions(periodicDimensions),
63  m_particlePoolPtr(particlePoolPtr),
64  m_clonedParticleSet(),
65  m_circGridWidth(1),
66  m_periodicDimIndex(-1)
67  {
69  if (circBorderWidth <= 0.0) {
70  circBorderWidth = this->getGridSpacing();
71  }
72  setCircularBorderWidth(circBorderWidth, this->getGridSpacing());
73  }
74 
75  template <class TmplParticle>
77  {
78  if (m_periodicDimensions.size() != 3) {
79  std::stringstream msg;
80  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());
86  }
87  int numPeriodic = 0;
88  for (int i = 0; i < 3; i++) {
89  if (m_periodicDimensions[i])
90  {
91  m_periodicDimIndex = i;
92  numPeriodic++;
93  }
94  }
95 
96  if (numPeriodic > 1) {
97  std::stringstream msg;
98  msg
99  << "CircularNeighbourTable::CircularNeighbourTable -"
100  << " only a single dimension may be periodic.";
101  throw std::runtime_error(msg.str());
102  }
103  }
104 
105  template <class TmplParticle>
107  {
108  }
109 
110  template <class TmplParticle>
112  double circBorderWidth,
113  double gridSpacing
114  )
115  {
116  m_circGridWidth = int(ceil(circBorderWidth/gridSpacing));
117  }
118 
119  template <class TmplParticle>
121  double circBorderWidth
122  )
123  {
124  setCircularBorderWidth(circBorderWidth, this->getGridSpacing());
125  }
126 
127  template <class TmplParticle>
129  const BoundingBox &bBox,
130  double gridSpacing,
131  double circBorderWidth
132  )
133  {
134  if (havePeriodicDimensions())
135  {
136  circBorderWidth =
137  min(
138  (bBox.getSizes()[m_periodicDimIndex])/2.0,
139  circBorderWidth
140  );
141  }
142  setCircularBorderWidth(circBorderWidth, gridSpacing);
143  ParticleVector particles = getNonClonedParticles();
144  clearClonedParticles();
145  this->clearAndRecomputeGrid(bBox, gridSpacing);
146  for (
147  typename ParticleVector::iterator it = particles.begin();
148  it != particles.end();
149  it++
150  )
151  {
152  this->insert(*it);
153  }
154  }
155 
156  template <class TmplParticle>
158  const BoundingBox &bBox,
159  double gridSpacing
160  )
161  {
162  this->resize(bBox, gridSpacing, gridSpacing);
163  }
164 
165  template <class TmplParticle>
167  Particle *pParticle,
168  const Vec3 &newPosition
169  )
170  {
171  Particle *pNewParticle = m_particlePoolPtr->construct(*pParticle);
172  pNewParticle->moveTo(newPosition);
173  Inherited::insert(pNewParticle);
174  m_clonedParticleSet.insert(pNewParticle);
175  }
176 
177  template <class TmplParticle>
179  {
180  return (m_periodicDimIndex >= 0);
181  }
182 
183  template <class TmplParticle>
185  const Vec3 &posn
186  ) const
187  {
188  if (
189  havePeriodicDimensions()
190  )
191  {
192  const int &dimIdx = m_periodicDimIndex;
193  if (
194  (posn[dimIdx] < this->getBBox().getMinPt()[dimIdx])
195  ||
196  (posn[dimIdx] > this->getBBox().getMaxPt()[dimIdx])
197  )
198  {
199  Vec3 moddedPosn = posn;
200  const double dimSize = this->getBBox().getSizes()[dimIdx];
201  moddedPosn[dimIdx] -= this->getBBox().getMinPt()[dimIdx];
202  moddedPosn[dimIdx] =
203  (moddedPosn[dimIdx] > 0.0)
204  ?
205  (
206  this->getBBox().getMinPt()[dimIdx]
207  +
208  moddedPosn[dimIdx] - (floor(moddedPosn[dimIdx]/dimSize)*dimSize)
209  )
210  :
211  (
212  this->getBBox().getMaxPt()[dimIdx]
213  -
214  (fabs(moddedPosn[dimIdx]) - (floor(fabs(moddedPosn[dimIdx])/dimSize)*dimSize))
215  );
216 
217  return moddedPosn;
218  }
219  }
220  return posn;
221  }
222 
223  template <class TmplParticle>
225  {
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());
229 
230  this->insertInTable(pParticle, minIdx, maxIdx);
231  this->addInserted(pParticle);
232 
233  if (havePeriodicDimensions())
234  {
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);
241  }
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);
246  }
247  }
248  }
249  }
250  }
251 
252  template <class TmplParticle>
254  {
255  this->insert(&particle);
256  }
257 
258  template <class TmplParticle>
260  {
261  return m_clonedParticleSet.size();
262  }
263 
264  template <class TmplParticle>
266  {
267  return this->size() - getNumClonedParticles();
268  }
269 
270  template <class TmplParticle>
273  {
274  return m_periodicDimensions;
275  }
276 
277  template <class TmplParticle>
279  Particle *p
280  ) const
281  {
282  return (m_clonedParticleSet.find(p) != m_clonedParticleSet.end());
283  }
284 
285  template <class TmplParticle>
288  {
289  ParticleVector all = this->getInsertedParticles();
290  ParticleVector nonCloned;
291  nonCloned.reserve(all.size()/2);
292  for (
293  typename ParticleVector::iterator it = all.begin();
294  it != all.end();
295  it++
296  )
297  {
298  if (!isClone(*it))
299  {
300  nonCloned.push_back(*it);
301  }
302  }
303  return nonCloned;
304  }
305 
306  template <class TmplParticle>
308  {
309  for (
310  typename ParticleSet::iterator it = m_clonedParticleSet.begin();
311  it != m_clonedParticleSet.end();
312  it++
313  )
314  {
315  m_particlePoolPtr->destroy(*it);
316  }
317  m_clonedParticleSet.clear();
318  }
319  }
320 }
321 
322 #endif