ESyS-Particle  4.0.1
CircularNeighbourTable.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_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 = getVecIndex(pParticle->getPos() - pParticle->getRad());
00228       const Vec3L maxIdx = 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