ESyS-Particle  2.3.2
SphereNeighbours.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 #include "Foundation/StringUtil.h"
15 #include <functional>
16 #include <limits>
17 
18 namespace esys
19 {
20  namespace lsm
21  {
22  template <typename TmplSphere, typename TmplIdPair>
24  double maxDist,
25  const BoundingBox &bBox,
26  const BoolVector &circDimensions
27  )
28  : m_connectionPoolPtr(new IdPairPool(4096)),
29  m_connectionSet(),
30  m_nTablePtr(),
31  m_minRadius(std::numeric_limits<double>::max()),
32  m_maxRadius(-std::numeric_limits<double>::max()),
33  m_maxDist(maxDist),
34  m_minPt(bBox.getMinPt()),
35  m_maxPt(bBox.getMaxPt())
36  {
37  const double gridSize =
38  (
39  max(
40  bBox.getSizes()[0],
41  max(
42  bBox.getSizes()[1],
43  bBox.getSizes()[2]
44  )
45  )
46  )/5.0;
47  m_nTablePtr =
48  NTablePtr(
49  new NTable(
50  bBox,
51  gridSize,
52  circDimensions,
53  2*gridSize
54  )
55  );
56  }
57 
58  template <typename TmplSphere, typename TmplIdPair>
60  {
61  }
62 
63  template <typename TmplSphere, typename TmplIdPair>
65  {
66  return m_nTablePtr->size();
67  }
68 
69  template <typename TmplSphere, typename TmplIdPair>
71  {
72  return m_connectionSet.size();
73  }
74 
75  template <typename TmplSphere, typename TmplIdPair>
77  {
78  return m_minRadius;
79  }
80 
81  template <typename TmplSphere, typename TmplIdPair>
83  {
84  return m_maxRadius;
85  }
86 
87  template <typename TmplSphere, typename TmplIdPair>
90  {
91  return m_nTablePtr->getIterator();
92  }
93 
94  template <typename TmplSphere, typename TmplIdPair>
97  const Sphere &p1,
98  const Sphere &p2
99  )
100  {
101  return
102  **(m_connectionSet.insert(
103  m_connectionPoolPtr->construct(p1.getId(), p2.getId())
104  ).first);
105  }
106 
107  template <typename TmplSphere>
109  {
110  public:
111  bool operator()(const TmplSphere &p1, const TmplSphere &p2) const
112  {
113  return (p1.getId() < p2.getId());
114  }
115 
116  bool operator()(const TmplSphere *p1, const TmplSphere *p2) const
117  {
118  return (p1->getId() < p2->getId());
119  }
120  };
121 
122  template <typename TmplType>
123  class Deref
124  {
125  public:
126  typedef const TmplType& result_type;
127  typedef const TmplType* argument_type;
128 
131  { return *a; }
132  };
133 
134  template <typename TmplSphere, typename TmplIdPair>
135  template <typename TmplSphereIterator>
136  typename SphereNeighbours<TmplSphere,TmplIdPair>::IdPairVector
138  TmplSphereIterator it
139  )
140  {
141  // Put the spheres in the neighbour table
142  typedef std::set<Sphere *, CmpSphereId<Sphere> > SphereSet;
143  SphereSet pSet;
144  while (it.hasNext())
145  {
146  Sphere &p = it.next();
147  insert(p);
148  pSet.insert(&p);
149  }
150  ConstIdPairSet idPairSet;
151  // Resize ntable according to min and max sphere radii.
152  m_nTablePtr->resize(
153  getSphereBBox(),
154  4.1*getMinRadius(),
155  2.1*getMaxRadius()
156  );
157 
158  // For each sphere in the iterator it, determine it's
159  // neighours within m_maxDist distance.
160  for (
161  typename SphereSet::const_iterator pIt = pSet.begin();
162  pIt != pSet.end();
163  pIt++
164  )
165  {
166  // get the vector of spheres which are contained in nearby cells.
167  typename NTable::ParticleVector nVector =
168  m_nTablePtr->getNeighbourVector(
169  (*pIt)->getPos(),
170  (*pIt)->getRad() + m_maxDist
171  );
172  // for each of these cell-spheres, determine if they
173  // are closer than m_maxDist.
174  for (
175  typename NTable::ParticleVector::const_iterator nIt = nVector.begin();
176  nIt != nVector.end();
177  nIt++
178  )
179  {
180  Sphere *p1 = (*pIt);
181  Sphere *p2 = (*nIt);
182 
183  if (
184  (
185  (p1->getId() < p2->getId())
186  &&
187  (pSet.find(p1) != pSet.end())
188  &&
189  (pSet.find(p2) != pSet.end())
190  )
191  ||
192  (
193  ((pSet.find(p1)==pSet.end()) && (pSet.find(p2)!= pSet.end()))
194  ||
195  ((pSet.find(p1)!=pSet.end()) && (pSet.find(p2)== pSet.end()))
196  )
197  )
198  {
199  p1 =
200  ((*pIt)->getId() < (*nIt)->getId())
201  ?
202  (*pIt)
203  :
204  (*nIt);
205  p2 =
206  ((*pIt)->getId() < (*nIt)->getId())
207  ?
208  (*nIt)
209  :
210  (*pIt);
211  const double radiusSumPlusTol =
212  m_maxDist + p1->getRad() + p2->getRad();
213  const double radiusSumPlusTolSqrd =
214  radiusSumPlusTol*radiusSumPlusTol;
215 
216  if (
217  (p1->getPos() - p2->getPos()).norm2()
218  <=
219  (radiusSumPlusTolSqrd)
220  )
221  {
222  idPairSet.insert(&createIdPair(*p1, *p2));
223  }
224  }
225  }
226  }
227  IdPairVector idPairVector;
228  idPairVector.reserve(idPairSet.size());
229  std::transform(
230  idPairSet.begin(),
231  idPairSet.end(),
232  std::back_insert_iterator<IdPairVector>(idPairVector),
233  Deref<IdPair>()
234  );
235  return idPairVector;
236  }
237 
238  template <typename TmplSphere, typename TmplIdPair>
239  void
241  {
242  if (p.getRad() < m_minRadius)
243  {
244  m_minRadius = p.getRad();
245  }
246  if (p.getRad() > m_maxRadius)
247  {
248  m_maxRadius = p.getRad();
249  }
250 
251  m_nTablePtr->insert(p);
252 
253  for (int i = 0; i < 3; i++)
254  {
255  if (!(m_nTablePtr->getPeriodicDimensions()[i]))
256  {
257  if (p.getPos()[i]-p.getRad() < m_minPt[i])
258  {
259  m_minPt[i] = p.getPos()[i]-p.getRad();
260  }
261  if (p.getPos()[i]+p.getRad() > m_maxPt[i])
262  {
263  m_maxPt[i] = p.getPos()[i]+p.getRad();
264  }
265  }
266  }
267  }
268 
269  template <typename TmplSphere, typename TmplIdPair>
272  {
273  return BoundingBox(m_minPt, m_maxPt);
274  }
275  }
276 }