ESyS-Particle  4.0.1
ClosePackIterator.hpp
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_LSMCLOSEPACKITERATOR_HPP
00015 #define ESYS_LSMCLOSEPACKITERATOR_HPP
00016 
00017 namespace esys
00018 {
00019   namespace lsm
00020   {
00021     template <int NI, int NJ, int NK>
00022     TmplMatrix<NI,NJ,NK>::TmplMatrix()
00023     {
00024       for (int i = 0; i < NI; i++)
00025       {
00026         for (int j = 0; j < NJ; j++)
00027         {
00028           for (int k = 0; k < NJ; k++)
00029           {
00030             m_matrix[i][j][k] = 0.0;
00031           }
00032         }
00033       }
00034     }
00035 
00036     template <int NI, int NJ, int NK>
00037     TmplMatrix<NI,NJ,NK>::TmplMatrix(const TmplMatrix &m)
00038     {
00039       for (int i = 0; i < NI; i++)
00040       {
00041         for (int j = 0; j < NJ; j++)
00042         {
00043           for (int k = 0; k < NJ; k++)
00044           {
00045             m_matrix[i][j][k] = m(i,j,k);
00046           }
00047         }
00048       }
00049     }
00050 
00051     template <int NI, int NJ, int NK>
00052     TmplMatrix<NI,NJ,NK> &
00053     TmplMatrix<NI,NJ,NK>::TmplMatrix::operator=(const TmplMatrix &m)
00054     {
00055       for (int i = 0; i < NI; i++)
00056       {
00057         for (int j = 0; j < NJ; j++)
00058         {
00059           for (int k = 0; k < NJ; k++)
00060           {
00061             m_matrix[i][j][k] = m(i,j,k);
00062           }
00063         }
00064       }
00065       return *this;
00066     }
00067 
00068     template <int NI, int NJ, int NK>
00069     const double &TmplMatrix<NI,NJ,NK>::operator()(int i, int j, int k) const
00070     {
00071       return m_matrix[i][j][k];
00072     }
00073 
00074     template <int NI, int NJ, int NK>
00075     double &TmplMatrix<NI,NJ,NK>::operator()(int i, int j, int k)
00076     {
00077       return m_matrix[i][j][k];
00078     }
00079 
00080     template <int NI, int NJ, int NK>
00081     int TmplMatrix<NI,NJ,NK>::getNumI() const
00082     {
00083       return NI;
00084     }
00085 
00086     template <int NI, int NJ, int NK>
00087     int TmplMatrix<NI,NJ,NK>::getNumJ() const
00088     {
00089       return NJ;
00090     }
00091 
00092     template <int NI, int NJ, int NK>
00093     int TmplMatrix<NI,NJ,NK>::getNumK() const
00094     {
00095       return NK;
00096     }
00097 
00098     ClosePackIterator::ClosePackIterator()
00099       : m_radius(),
00100         m_minPt(),
00101         m_offsetMatrix(),
00102         m_dimRepeat(),
00103         m_dimCount(),
00104         m_dimIdx(),
00105         m_dim()
00106     {
00107     }
00108 
00109     ClosePackIterator::ClosePackIterator(
00110       int numI,
00111       int numJ,
00112       int numK,
00113       double sphereRadius,
00114       ClosePackOrientation orientation
00115     )
00116       : m_radius(sphereRadius),
00117         m_minPt(),
00118         m_offsetMatrix(),
00119         m_dimRepeat(),
00120         m_dimCount(numI, numJ, numK),
00121         m_dimIdx(),
00122         m_dim(s_orientationDimMap[orientation])
00123     {
00124       for (int i = 0; i < 3; i++)
00125       {
00126         if (m_dimCount[i] <= 0)
00127         {
00128           m_dimCount[i] = 1;
00129         }
00130       }
00131     }
00132 
00133     void ClosePackIterator::setDimRepeat(const Vec3L &dimRepeat)
00134     {
00135       m_dimRepeat = dimRepeat;
00136     }
00137 
00138     void ClosePackIterator::setOffsetMatrix(const OffsetMatrix &offsetMatrix)
00139     {
00140       m_offsetMatrix = offsetMatrix;
00141     }
00142 
00143     double ClosePackIterator::getRadius() const
00144     {
00145       return m_radius;
00146     }
00147     
00148     bool ClosePackIterator::hasNext() const
00149     {
00150        return (m_dimIdx[2] < m_dimCount[2]);
00151     }
00152 
00153     const Vec3 &ClosePackIterator::getMinPt() const
00154     {
00155       return m_minPt;
00156     }
00157 
00158     double ClosePackIterator::getOffset(int i) const
00159     {
00160       return
00161         m_offsetMatrix(
00162           i,
00163           (m_dimIdx[(i+1) % 3]) % m_dimRepeat[i],
00164           (m_dimIdx[(i+2) % 3]) % m_dimRepeat[i]
00165         );
00166     }
00167 
00168     void ClosePackIterator::incrementDimIndex()
00169     {
00170       (m_dimIdx[0])++;
00171       if (m_dimIdx[0] >= m_dimCount[0])
00172       {
00173         m_dimIdx[0] = 0;
00174         (m_dimIdx[1])++;
00175         if (m_dimIdx[1] >= m_dimCount[1])
00176         {
00177           m_dimIdx[1] = 0;
00178           (m_dimIdx[2])++;
00179         }
00180       }
00181     }
00182     
00183     Vec3 ClosePackIterator::next()
00184     {
00185       Vec3 centre;
00186       centre[m_dim[0]] = getOffset(0) + m_dimIdx[0]*(m_radius*2.0);
00187       centre[m_dim[1]] = getOffset(1) + m_dimIdx[1]*(m_radius*SQRT_3);
00188       centre[m_dim[2]] = getOffset(2) + m_dimIdx[2]*(m_radius*SQRT_8_OVER_3);
00189       incrementDimIndex();
00190       return (getMinPt() + centre);
00191     }
00192   }
00193 }
00194 
00195 #endif