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_LSMCIRCULARNEIGHBOURTABLE_HPP 00015 #define ESYS_LSMCIRCULARNEIGHBOURTABLE_HPP 00016 00017 #include "Geometry/NeighbourTable.h" 00018 #include <boost/pool/object_pool.hpp> 00019 #include <boost/shared_ptr.hpp> 00020 00021 #include <sstream> 00022 #include <stdexcept> 00023 #include <set> 00024 00025 namespace esys 00026 { 00027 namespace lsm 00028 { 00032 template <class TmplParticle> 00033 CircularNeighbourTable<TmplParticle>::CircularNeighbourTable( 00034 const BoundingBox &bBox, 00035 double gridSpacing, 00036 const BoolVector &periodicDimensions, 00037 double circBorderWidth 00038 ) 00039 : Inherited(bBox, gridSpacing), 00040 m_periodicDimensions(periodicDimensions), 00041 m_particlePoolPtr(new ParticlePool(4096)), 00042 m_clonedParticleSet(), 00043 m_circGridWidth(1), 00044 m_periodicDimIndex(-1) 00045 { 00046 checkPeriodicDimensions(); 00047 if (circBorderWidth <= 0.0) { 00048 circBorderWidth = this->getGridSpacing(); 00049 } 00050 setCircularBorderWidth(circBorderWidth, this->getGridSpacing()); 00051 } 00052 00053 template <class TmplParticle> 00054 CircularNeighbourTable<TmplParticle>::CircularNeighbourTable( 00055 const BoundingBox &bBox, 00056 double gridSpacing, 00057 ParticlePoolPtr particlePoolPtr, 00058 const BoolVector &periodicDimensions, 00059 double circBorderWidth 00060 ) 00061 : Inherited(bBox, gridSpacing), 00062 m_periodicDimensions(periodicDimensions), 00063 m_particlePoolPtr(particlePoolPtr), 00064 m_clonedParticleSet(), 00065 m_circGridWidth(1), 00066 m_periodicDimIndex(-1) 00067 { 00068 checkPeriodicDimensions(); 00069 if (circBorderWidth <= 0.0) { 00070 circBorderWidth = this->getGridSpacing(); 00071 } 00072 setCircularBorderWidth(circBorderWidth, this->getGridSpacing()); 00073 } 00074 00075 template <class TmplParticle> 00076 void CircularNeighbourTable<TmplParticle>::checkPeriodicDimensions() 00077 { 00078 if (m_periodicDimensions.size() != 3) { 00079 std::stringstream msg; 00080 msg 00081 << "CircularNeighbourTable::CircularNeighbourTable -" 00082 << " size of periodic dimensions argument (" 00083 << m_periodicDimensions.size() 00084 << ") is not equal to 3"; 00085 throw std::runtime_error(msg.str()); 00086 } 00087 int numPeriodic = 0; 00088 for (int i = 0; i < 3; i++) { 00089 if (m_periodicDimensions[i]) 00090 { 00091 m_periodicDimIndex = i; 00092 numPeriodic++; 00093 } 00094 } 00095 00096 if (numPeriodic > 1) { 00097 std::stringstream msg; 00098 msg 00099 << "CircularNeighbourTable::CircularNeighbourTable -" 00100 << " only a single dimension may be periodic."; 00101 throw std::runtime_error(msg.str()); 00102 } 00103 } 00104 00105 template <class TmplParticle> 00106 CircularNeighbourTable<TmplParticle>::~CircularNeighbourTable() 00107 { 00108 } 00109 00110 template <class TmplParticle> 00111 void CircularNeighbourTable<TmplParticle>::setCircularBorderWidth( 00112 double circBorderWidth, 00113 double gridSpacing 00114 ) 00115 { 00116 m_circGridWidth = int(ceil(circBorderWidth/gridSpacing)); 00117 } 00118 00119 template <class TmplParticle> 00120 void CircularNeighbourTable<TmplParticle>::setCircularBorderWidth( 00121 double circBorderWidth 00122 ) 00123 { 00124 setCircularBorderWidth(circBorderWidth, this->getGridSpacing()); 00125 } 00126 00127 template <class TmplParticle> 00128 void CircularNeighbourTable<TmplParticle>::resize( 00129 const BoundingBox &bBox, 00130 double gridSpacing, 00131 double circBorderWidth 00132 ) 00133 { 00134 if (havePeriodicDimensions()) 00135 { 00136 circBorderWidth = 00137 min( 00138 (bBox.getSizes()[m_periodicDimIndex])/2.0, 00139 circBorderWidth 00140 ); 00141 } 00142 setCircularBorderWidth(circBorderWidth, gridSpacing); 00143 ParticleVector particles = getNonClonedParticles(); 00144 clearClonedParticles(); 00145 this->clearAndRecomputeGrid(bBox, gridSpacing); 00146 for ( 00147 typename ParticleVector::iterator it = particles.begin(); 00148 it != particles.end(); 00149 it++ 00150 ) 00151 { 00152 this->insert(*it); 00153 } 00154 } 00155 00156 template <class TmplParticle> 00157 void CircularNeighbourTable<TmplParticle>::resize( 00158 const BoundingBox &bBox, 00159 double gridSpacing 00160 ) 00161 { 00162 this->resize(bBox, gridSpacing, gridSpacing); 00163 } 00164 00165 template <class TmplParticle> 00166 void CircularNeighbourTable<TmplParticle>::insertClone( 00167 Particle *pParticle, 00168 const Vec3 &newPosition 00169 ) 00170 { 00171 Particle *pNewParticle = m_particlePoolPtr->construct(*pParticle); 00172 pNewParticle->moveTo(newPosition); 00173 Inherited::insert(pNewParticle); 00174 m_clonedParticleSet.insert(pNewParticle); 00175 } 00176 00177 template <class TmplParticle> 00178 bool CircularNeighbourTable<TmplParticle>::havePeriodicDimensions() const 00179 { 00180 return (m_periodicDimIndex >= 0); 00181 } 00182 00183 template <class TmplParticle> 00184 Vec3 CircularNeighbourTable<TmplParticle>::getModdedPosn( 00185 const Vec3 &posn 00186 ) const 00187 { 00188 if ( 00189 havePeriodicDimensions() 00190 ) 00191 { 00192 const int &dimIdx = m_periodicDimIndex; 00193 if ( 00194 (posn[dimIdx] < this->getBBox().getMinPt()[dimIdx]) 00195 || 00196 (posn[dimIdx] > this->getBBox().getMaxPt()[dimIdx]) 00197 ) 00198 { 00199 Vec3 moddedPosn = posn; 00200 const double dimSize = this->getBBox().getSizes()[dimIdx]; 00201 moddedPosn[dimIdx] -= this->getBBox().getMinPt()[dimIdx]; 00202 moddedPosn[dimIdx] = 00203 (moddedPosn[dimIdx] > 0.0) 00204 ? 00205 ( 00206 this->getBBox().getMinPt()[dimIdx] 00207 + 00208 moddedPosn[dimIdx] - (floor(moddedPosn[dimIdx]/dimSize)*dimSize) 00209 ) 00210 : 00211 ( 00212 this->getBBox().getMaxPt()[dimIdx] 00213 - 00214 (fabs(moddedPosn[dimIdx]) - (floor(fabs(moddedPosn[dimIdx])/dimSize)*dimSize)) 00215 ); 00216 00217 return moddedPosn; 00218 } 00219 } 00220 return posn; 00221 } 00222 00223 template <class TmplParticle> 00224 void CircularNeighbourTable<TmplParticle>::insert(Particle *pParticle) 00225 { 00226 pParticle->moveTo(getModdedPosn(pParticle->getPos())); 00227 const Vec3L minIdx = this->getVecIndex(pParticle->getPos() - pParticle->getRad()); 00228 const Vec3L maxIdx = this->getVecIndex(pParticle->getPos() + pParticle->getRad()); 00229 00230 this->insertInTable(pParticle, minIdx, maxIdx); 00231 this->addInserted(pParticle); 00232 00233 if (havePeriodicDimensions()) 00234 { 00235 for (int dimIdx = 0; dimIdx < 3; dimIdx++) { 00236 if (m_periodicDimensions[dimIdx]) { 00237 if (minIdx[dimIdx] < (this->getMinVecIndex()[dimIdx] + m_circGridWidth)) { 00238 Vec3 shift = Vec3::ZERO; 00239 shift[dimIdx] = this->getBBox().getSizes()[dimIdx]; 00240 this->insertClone(pParticle, pParticle->getPos() + shift); 00241 } 00242 if (maxIdx[dimIdx] > (this->getMaxVecIndex()[dimIdx] - m_circGridWidth)) { 00243 Vec3 shift = Vec3::ZERO; 00244 shift[dimIdx] = this->getBBox().getSizes()[dimIdx]; 00245 this->insertClone(pParticle, pParticle->getPos() - shift); 00246 } 00247 } 00248 } 00249 } 00250 } 00251 00252 template <class TmplParticle> 00253 void CircularNeighbourTable<TmplParticle>::insert(Particle &particle) 00254 { 00255 this->insert(&particle); 00256 } 00257 00258 template <class TmplParticle> 00259 size_t CircularNeighbourTable<TmplParticle>::getNumClonedParticles() const 00260 { 00261 return m_clonedParticleSet.size(); 00262 } 00263 00264 template <class TmplParticle> 00265 size_t CircularNeighbourTable<TmplParticle>::getNumParticles() const 00266 { 00267 return this->size() - getNumClonedParticles(); 00268 } 00269 00270 template <class TmplParticle> 00271 const typename CircularNeighbourTable<TmplParticle>::BoolVector & 00272 CircularNeighbourTable<TmplParticle>::getPeriodicDimensions() const 00273 { 00274 return m_periodicDimensions; 00275 } 00276 00277 template <class TmplParticle> 00278 bool CircularNeighbourTable<TmplParticle>::isClone( 00279 Particle *p 00280 ) const 00281 { 00282 return (m_clonedParticleSet.find(p) != m_clonedParticleSet.end()); 00283 } 00284 00285 template <class TmplParticle> 00286 typename CircularNeighbourTable<TmplParticle>::ParticleVector 00287 CircularNeighbourTable<TmplParticle>::getNonClonedParticles() 00288 { 00289 ParticleVector all = this->getInsertedParticles(); 00290 ParticleVector nonCloned; 00291 nonCloned.reserve(all.size()/2); 00292 for ( 00293 typename ParticleVector::iterator it = all.begin(); 00294 it != all.end(); 00295 it++ 00296 ) 00297 { 00298 if (!isClone(*it)) 00299 { 00300 nonCloned.push_back(*it); 00301 } 00302 } 00303 return nonCloned; 00304 } 00305 00306 template <class TmplParticle> 00307 void CircularNeighbourTable<TmplParticle>::clearClonedParticles() 00308 { 00309 for ( 00310 typename ParticleSet::iterator it = m_clonedParticleSet.begin(); 00311 it != m_clonedParticleSet.end(); 00312 it++ 00313 ) 00314 { 00315 m_particlePoolPtr->destroy(*it); 00316 } 00317 m_clonedParticleSet.clear(); 00318 } 00319 } 00320 } 00321 00322 #endif