ESyS-Particle
4.0.1
|
00001 00002 // // 00003 // Copyright (c) 2003-2011 by The University of Queensland // 00004 // Earth Systems Science Computational Centre (ESSCC) // 00005 // http://www.uq.edu.au/esscc // 00006 // // 00007 // Primary Business: Brisbane, Queensland, Australia // 00008 // Licensed under the Open Software License version 3.0 // 00009 // http://www.opensource.org/licenses/osl-3.0.php // 00010 // // 00012 00013 00014 #ifndef ESYS_LSMGRIDITERATOR_H 00015 #define ESYS_LSMGRIDITERATOR_H 00016 00017 #include "Foundation/BoundingBox.h" 00018 #include "Foundation/vec3.h" 00019 00020 namespace esys 00021 { 00022 namespace lsm 00023 { 00028 class GridIterator 00029 { 00030 public: 00031 inline GridIterator() 00032 : m_sphereRadius(0.0), 00033 m_minI(0), 00034 m_minJ(0), 00035 m_minK(0), 00036 m_maxI(0), 00037 m_maxJ(0), 00038 m_maxK(0), 00039 m_i(0), 00040 m_j(0), 00041 m_k(0), 00042 m_centrePtBBox() 00043 { 00044 } 00045 00046 inline GridIterator(int numPtsX, int numPtsY, int numPtsZ, double 00047 sphereRadius) 00048 : m_sphereRadius(sphereRadius), 00049 m_minI(0), 00050 m_minJ(0), 00051 m_minK(0), 00052 m_maxI(numPtsX), 00053 m_maxJ(numPtsZ), 00054 m_maxK(numPtsY), 00055 m_i(0), 00056 m_j(0), 00057 m_k(0), 00058 m_centrePtBBox() 00059 { 00060 Vec3 minPt; 00061 Vec3 maxPt; 00062 if (numPtsZ > 1) 00063 { 00064 minPt = Vec3(0,0,0) + sphereRadius; 00065 maxPt = 00066 Vec3( 00067 (numPtsX-1)*2.0*sphereRadius 00068 + 00069 ( 00070 (numPtsZ > 1) 00071 ? 00072 sphereRadius 00073 : 00074 0.0 00075 ) 00076 + 00077 ( 00078 (numPtsY > 1) 00079 ? 00080 sphereRadius 00081 : 00082 0.0 00083 ), 00084 (numPtsY-1)*(sphereRadius*2.0*sqrt(2.0/3.0)), 00085 (numPtsZ-1)*(sphereRadius*sqrt(3.0)) 00086 + 00087 ((numPtsY > 1) ? sphereRadius*sqrt(3.0)/3.0 : 0.0) 00088 ); 00089 } else { 00090 minPt = Vec3(sphereRadius, sphereRadius, 0); 00091 maxPt = 00092 Vec3( 00093 (numPtsX-1)*2.0*sphereRadius 00094 + 00095 ( 00096 (numPtsY > 1) 00097 ? 00098 sphereRadius : 00099 0.0 00100 ), 00101 (numPtsY-1)*(sphereRadius*sqrt(3.0)), 00102 0 00103 ); 00104 m_minJ = m_maxJ; 00105 } 00106 m_centrePtBBox = BoundingBox(minPt, (minPt + maxPt)); 00107 } 00108 00112 inline GridIterator(const BoundingBox &particleBBox, double sphereRadius, bool hard_limit=false) 00113 : m_sphereRadius(sphereRadius), 00114 m_minI(0), 00115 m_minJ(0), 00116 m_minK(0), 00117 m_maxI(0), 00118 m_maxJ(0), 00119 m_maxK(0), 00120 m_i(0), 00121 m_j(0), 00122 m_k(0), 00123 m_centrePtBBox() 00124 { 00125 int numPtsX ,numPtsY ,numPtsZ; 00126 if(!hard_limit){ // find "best fit" 00127 numPtsX = max(1, int(nearbyint((particleBBox.getSizes().X()-(sphereRadius/4.0))/(sphereRadius*2.0)))); 00128 numPtsY = max(1, int(nearbyint(particleBBox.getSizes().Y()/(sphereRadius*2.0*sqrt(2.0/3.0))))); 00129 numPtsZ = max(1, int(nearbyint(particleBBox.getSizes().Z()/(sphereRadius*sqrt(3.0))))); 00130 } else { // bounding box is hard limit 00131 numPtsX = max(1, int(floor((particleBBox.getSizes().X()-(sphereRadius/4.0))/(sphereRadius*2.0)))); 00132 numPtsY = max(1, int(floor(particleBBox.getSizes().Y()/(sphereRadius*2.0*sqrt(2.0/3.0))))); 00133 numPtsZ = max(1, int(floor(particleBBox.getSizes().Z()/(sphereRadius*sqrt(3.0))))); 00134 } 00135 if ((numPtsZ > 1) && (numPtsY > 1)) { 00136 numPtsX--; 00137 } 00138 if (particleBBox.getSizes().Z() > 0.0) { 00139 const Vec3 minPt = particleBBox.getMinPt() + sphereRadius; 00140 const Vec3 maxPt = 00141 Vec3( 00142 (numPtsX-1)*2.0*sphereRadius + ((numPtsZ > 1) ? sphereRadius : 0.0) + ((numPtsY > 1) ? sphereRadius : 0.0), 00143 (numPtsY-1)*(sphereRadius*2.0*sqrt(2.0/3.0)), 00144 (numPtsZ-1)*(sphereRadius*sqrt(3.0)) 00145 + 00146 ((numPtsY > 1) ? sphereRadius*sqrt(3.0)/3.0 : 0.0) 00147 ); 00148 m_centrePtBBox = BoundingBox(minPt, (minPt + maxPt)); 00149 } else { 00150 numPtsX = int(nearbyint((particleBBox.getSizes().X()-(sphereRadius/4.0))/(sphereRadius*2.0))); 00151 numPtsY = int(nearbyint(particleBBox.getSizes().Y()/(sphereRadius*sqrt(3.0)))); 00152 numPtsZ = 0; 00153 00154 Vec3 minPt = particleBBox.getMinPt() + sphereRadius; 00155 minPt.Z() = particleBBox.getMinPt().Z(); 00156 const Vec3 maxPt = 00157 Vec3( 00158 (numPtsX-1)*2.0*sphereRadius + ((numPtsY > 1) ? sphereRadius : 0.0), 00159 (numPtsY-1)*(sphereRadius*sqrt(3.0)), 00160 particleBBox.getMaxPt().Z() 00161 ); 00162 m_centrePtBBox = BoundingBox(minPt, (minPt + maxPt)); 00163 } 00164 m_minI = 0; 00165 m_minJ = 0; 00166 m_minK = 0; 00167 m_maxI = numPtsX; 00168 m_maxK = numPtsY; 00169 m_maxJ = numPtsZ; 00170 00171 m_i = m_minI; 00172 m_j = m_minJ; 00173 m_k = m_minK; 00174 } 00175 00176 inline ~GridIterator() 00177 { 00178 } 00179 00180 inline const BoundingBox &getBoundingBox() const 00181 { 00182 return m_centrePtBBox; 00183 } 00184 00185 inline BoundingBox getSphereBBox() const 00186 { 00187 return 00188 ( 00189 is2d() 00190 ? 00191 BoundingBox( 00192 m_centrePtBBox.getMinPt() - Vec3(m_sphereRadius, m_sphereRadius, 0.0), 00193 m_centrePtBBox.getMaxPt() + Vec3(m_sphereRadius, m_sphereRadius, 0.0) 00194 ) 00195 : 00196 BoundingBox( 00197 m_centrePtBBox.getMinPt() - m_sphereRadius, 00198 m_centrePtBBox.getMaxPt() + m_sphereRadius 00199 ) 00200 ); 00201 } 00202 00207 inline bool hasNext() const 00208 { 00209 return (m_i < m_maxI); 00210 } 00211 00212 inline bool is2d() const 00213 { 00214 return (m_minJ == m_maxJ); 00215 } 00216 00217 inline Vec3 getPoint() const 00218 { 00219 if (is2d()) { 00220 return 00221 Vec3( 00222 m_centrePtBBox.getMinPt().X() + (double(m_i)+0.5*double(m_k%2))*m_sphereRadius*2.0, 00223 m_centrePtBBox.getMinPt().Y() + double(m_k)*sqrt(3.0)*m_sphereRadius, 00224 0.0 00225 ); 00226 } else { 00227 return 00228 Vec3( 00229 m_centrePtBBox.getMinPt().X() + (double(m_i)+0.5*double(m_j%2)+0.5*double(m_k%2))*m_sphereRadius*2.0, 00230 m_centrePtBBox.getMinPt().Y() + (double(m_k)*2.0*sqrt(2.0/3.0))*m_sphereRadius, 00231 m_centrePtBBox.getMinPt().Z() + ((double(m_j)+double(m_k%2)/3.0)*sqrt(3.0))*m_sphereRadius 00232 ); 00233 } 00234 } 00235 00236 inline void increment() 00237 { 00238 if (m_k+1 < m_maxK) { 00239 m_k++; 00240 } 00241 else if (m_j+1 < m_maxJ) { 00242 m_k = m_minK; 00243 m_j++; 00244 } 00245 else { 00246 m_k = m_minK; 00247 m_j = m_minJ; 00248 m_i++; 00249 } 00250 } 00251 00255 inline Vec3 next() 00256 { 00257 const Vec3 pt = getPoint(); 00258 increment(); 00259 return pt; 00260 } 00261 00262 private: 00263 double m_sphereRadius; 00264 00265 int m_minI; 00266 int m_minJ; 00267 int m_minK; 00268 00269 int m_maxI; 00270 int m_maxJ; 00271 int m_maxK; 00272 00273 int m_i; 00274 int m_j; 00275 int m_k; 00276 00277 BoundingBox m_centrePtBBox; 00278 }; 00279 }; 00280 }; 00281 00282 #endif