ESyS-Particle  2.3.2
GridIterator.h
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_LSMGRIDITERATOR_H
15 #define ESYS_LSMGRIDITERATOR_H
16 
17 #include "Foundation/BoundingBox.h"
18 #include "Foundation/vec3.h"
19 
20 namespace esys
21 {
22  namespace lsm
23  {
29  {
30  public:
31  inline GridIterator()
32  : m_sphereRadius(0.0),
33  m_minI(0),
34  m_minJ(0),
35  m_minK(0),
36  m_maxI(0),
37  m_maxJ(0),
38  m_maxK(0),
39  m_i(0),
40  m_j(0),
41  m_k(0),
43  {
44  }
45 
46  inline GridIterator(int numPtsX, int numPtsY, int numPtsZ, double
47 sphereRadius)
48  : m_sphereRadius(sphereRadius),
49  m_minI(0),
50  m_minJ(0),
51  m_minK(0),
52  m_maxI(numPtsX),
53  m_maxJ(numPtsZ),
54  m_maxK(numPtsY),
55  m_i(0),
56  m_j(0),
57  m_k(0),
59  {
60  Vec3 minPt;
61  Vec3 maxPt;
62  if (numPtsZ > 1)
63  {
64  minPt = Vec3(0,0,0) + sphereRadius;
65  maxPt =
66  Vec3(
67  (numPtsX-1)*2.0*sphereRadius
68  +
69  (
70  (numPtsZ > 1)
71  ?
72  sphereRadius
73  :
74  0.0
75  )
76  +
77  (
78  (numPtsY > 1)
79  ?
80  sphereRadius
81  :
82  0.0
83  ),
84  (numPtsY-1)*(sphereRadius*2.0*sqrt(2.0/3.0)),
85  (numPtsZ-1)*(sphereRadius*sqrt(3.0))
86  +
87  ((numPtsY > 1) ? sphereRadius*sqrt(3.0)/3.0 : 0.0)
88  );
89  } else {
90  minPt = Vec3(sphereRadius, sphereRadius, 0);
91  maxPt =
92  Vec3(
93  (numPtsX-1)*2.0*sphereRadius
94  +
95  (
96  (numPtsY > 1)
97  ?
98  sphereRadius :
99  0.0
100  ),
101  (numPtsY-1)*(sphereRadius*sqrt(3.0)),
102  0
103  );
104  m_minJ = m_maxJ;
105  }
106  m_centrePtBBox = BoundingBox(minPt, (minPt + maxPt));
107  }
108 
112  inline GridIterator(const BoundingBox &particleBBox, double sphereRadius, bool hard_limit=false)
113  : m_sphereRadius(sphereRadius),
114  m_minI(0),
115  m_minJ(0),
116  m_minK(0),
117  m_maxI(0),
118  m_maxJ(0),
119  m_maxK(0),
120  m_i(0),
121  m_j(0),
122  m_k(0),
124  {
125  int numPtsX ,numPtsY ,numPtsZ;
126  if(!hard_limit){ // find "best fit"
127  numPtsX = max(1, int(nearbyint((particleBBox.getSizes().X()-(sphereRadius/4.0))/(sphereRadius*2.0))));
128  numPtsY = max(1, int(nearbyint(particleBBox.getSizes().Y()/(sphereRadius*2.0*sqrt(2.0/3.0)))));
129  numPtsZ = max(1, int(nearbyint(particleBBox.getSizes().Z()/(sphereRadius*sqrt(3.0)))));
130  } else { // bounding box is hard limit
131  numPtsX = max(1, int(floor((particleBBox.getSizes().X()-(sphereRadius/4.0))/(sphereRadius*2.0))));
132  numPtsY = max(1, int(floor(particleBBox.getSizes().Y()/(sphereRadius*2.0*sqrt(2.0/3.0)))));
133  numPtsZ = max(1, int(floor(particleBBox.getSizes().Z()/(sphereRadius*sqrt(3.0)))));
134  }
135  if ((numPtsZ > 1) && (numPtsY > 1)) {
136  numPtsX--;
137  }
138  if (particleBBox.getSizes().Z() > 0.0) {
139  const Vec3 minPt = particleBBox.getMinPt() + sphereRadius;
140  const Vec3 maxPt =
141  Vec3(
142  (numPtsX-1)*2.0*sphereRadius + ((numPtsZ > 1) ? sphereRadius : 0.0) + ((numPtsY > 1) ? sphereRadius : 0.0),
143  (numPtsY-1)*(sphereRadius*2.0*sqrt(2.0/3.0)),
144  (numPtsZ-1)*(sphereRadius*sqrt(3.0))
145  +
146  ((numPtsY > 1) ? sphereRadius*sqrt(3.0)/3.0 : 0.0)
147  );
148  m_centrePtBBox = BoundingBox(minPt, (minPt + maxPt));
149  } else {
150  numPtsX = int(nearbyint((particleBBox.getSizes().X()-(sphereRadius/4.0))/(sphereRadius*2.0)));
151  numPtsY = int(nearbyint(particleBBox.getSizes().Y()/(sphereRadius*sqrt(3.0))));
152  numPtsZ = 0;
153 
154  Vec3 minPt = particleBBox.getMinPt() + sphereRadius;
155  minPt.Z() = particleBBox.getMinPt().Z();
156  const Vec3 maxPt =
157  Vec3(
158  (numPtsX-1)*2.0*sphereRadius + ((numPtsY > 1) ? sphereRadius : 0.0),
159  (numPtsY-1)*(sphereRadius*sqrt(3.0)),
160  particleBBox.getMaxPt().Z()
161  );
162  m_centrePtBBox = BoundingBox(minPt, (minPt + maxPt));
163  }
164  m_minI = 0;
165  m_minJ = 0;
166  m_minK = 0;
167  m_maxI = numPtsX;
168  m_maxK = numPtsY;
169  m_maxJ = numPtsZ;
170 
171  m_i = m_minI;
172  m_j = m_minJ;
173  m_k = m_minK;
174  }
175 
176  inline ~GridIterator()
177  {
178  }
179 
180  inline const BoundingBox &getBoundingBox() const
181  {
182  return m_centrePtBBox;
183  }
184 
185  inline BoundingBox getSphereBBox() const
186  {
187  return
188  (
189  is2d()
190  ?
191  BoundingBox(
194  )
195  :
196  BoundingBox(
199  )
200  );
201  }
202 
207  inline bool hasNext() const
208  {
209  return (m_i < m_maxI);
210  }
211 
212  inline bool is2d() const
213  {
214  return (m_minJ == m_maxJ);
215  }
216 
217  inline Vec3 getPoint() const
218  {
219  if (is2d()) {
220  return
221  Vec3(
222  m_centrePtBBox.getMinPt().X() + (double(m_i)+0.5*double(m_k%2))*m_sphereRadius*2.0,
223  m_centrePtBBox.getMinPt().Y() + double(m_k)*sqrt(3.0)*m_sphereRadius,
224  0.0
225  );
226  } else {
227  return
228  Vec3(
229  m_centrePtBBox.getMinPt().X() + (double(m_i)+0.5*double(m_j%2)+0.5*double(m_k%2))*m_sphereRadius*2.0,
230  m_centrePtBBox.getMinPt().Y() + (double(m_k)*2.0*sqrt(2.0/3.0))*m_sphereRadius,
231  m_centrePtBBox.getMinPt().Z() + ((double(m_j)+double(m_k%2)/3.0)*sqrt(3.0))*m_sphereRadius
232  );
233  }
234  }
235 
236  inline void increment()
237  {
238  if (m_k+1 < m_maxK) {
239  m_k++;
240  }
241  else if (m_j+1 < m_maxJ) {
242  m_k = m_minK;
243  m_j++;
244  }
245  else {
246  m_k = m_minK;
247  m_j = m_minJ;
248  m_i++;
249  }
250  }
251 
255  inline Vec3 next()
256  {
257  const Vec3 pt = getPoint();
258  increment();
259  return pt;
260  }
261 
262  private:
264 
265  int m_minI;
266  int m_minJ;
267  int m_minK;
268 
269  int m_maxI;
270  int m_maxJ;
271  int m_maxK;
272 
273  int m_i;
274  int m_j;
275  int m_k;
276 
278  };
279  }
280 }
281 
282 #endif