ESyS-Particle  4.0.1
GridIterator.h
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