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_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