ESyS-Particle  4.0.1
RandomBoxPacker.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 #include "Foundation/console.h"
00015 #include "Foundation/StringUtil.h"
00016 #include "Geometry/RandomBoxPacker.h"
00017 #include "Geometry/ParticleComparer.h"
00018 #include "Geometry/SphereFitter.h"
00019 
00020 #include <algorithm>
00021 #include <stdexcept>
00022 #include <float.h>
00023 
00024 namespace esys
00025 {
00026   namespace lsm
00027   {
00028     template <typename TmplParticle>
00029     class PlaneComparer
00030     {
00031     public:
00032       typedef TmplParticle Particle;
00033       PlaneComparer(const Particle &particle)
00034         : m_pParticle(&particle)
00035       {
00036       }
00037       
00038       inline bool operator()(const Plane &plane1, const Plane &plane2) const
00039       {
00040         return 
00041           (
00042             plane1.sep(m_pParticle->getPos())
00043             <
00044             plane2.sep(m_pParticle->getPos())
00045           );
00046       }
00047     private:
00048       const Particle *m_pParticle;
00049     };
00050 
00051     template<typename TmplTraits>
00052     FittedParticleIterator<TmplTraits>::FittedParticleIterator(
00053       Packer            &packer,
00054       int               maxInsertionFailures,
00055       const PlaneVector &fitPlaneVector
00056     ) :
00057       m_pPacker(&packer),
00058       m_fitPlaneVector(fitPlaneVector),
00059       m_maxInsertionFailures(maxInsertionFailures),
00060       m_lastFailCount(0),
00061       m_successCount(0),
00062       m_next(Fitter::getInvalidParticle()),
00063       m_fitterPtrVector()
00064     {
00065       if (m_next.isValid())
00066       {
00067         throw
00068           std::runtime_error(
00069             std::string("isValid method returns true for INVALID particle:")
00070             +
00071             StringUtil::toString(m_next)
00072           );
00073       }
00074       initialiseFitterPtrVector();
00075     }
00076 
00077     template<typename TmplTraits>
00078     const typename FittedParticleIterator<TmplTraits>::PlaneVector &
00079     FittedParticleIterator<TmplTraits>::getFitPlaneVector() const
00080     {
00081       return m_fitPlaneVector;
00082     }
00083 
00084     template<typename TmplTraits>
00085     int FittedParticleIterator<TmplTraits>::getMaxInsertionFailures() const
00086     {
00087       return m_maxInsertionFailures;
00088     }
00089 
00090     template<typename TmplTraits>
00091     const typename FittedParticleIterator<TmplTraits>::Packer &
00092     FittedParticleIterator<TmplTraits>::getPacker() const
00093     {
00094       return *m_pPacker;
00095     }
00096 
00097     template<typename TmplTraits>
00098     typename FittedParticleIterator<TmplTraits>::Packer &
00099     FittedParticleIterator<TmplTraits>::getPacker()
00100     {
00101       return *m_pPacker;
00102     }
00103 
00104     template<typename TmplTraits>
00105     void
00106     FittedParticleIterator<TmplTraits>::initialiseFitterPtrVector()
00107     {
00108       FitterPtrVector fitters;
00109       fitters.push_back(FitterPtr(new Move2SurfaceFitter(getPacker())));
00110       if (getPacker().is2d())
00111       {
00112         fitters.push_back(FitterPtr(new TwoDFitter(getPacker())));
00113         if (getFitPlaneVector().size() > 0) {
00114           fitters.push_back(FitterPtr(new TwoDPlaneFitter(getPacker())));
00115         }
00116       }
00117       else
00118       {
00119         fitters.push_back(FitterPtr(new ThreeDFitter(getPacker())));
00120         if (getFitPlaneVector().size() > 0) {
00121           fitters.push_back(FitterPtr(new ThreeDPlaneFitter(getPacker())));
00122         }
00123       }
00124 
00125       m_fitterPtrVector = fitters;
00126     }
00127 
00128     template<typename TmplTraits>
00129     typename FittedParticleIterator<TmplTraits>::FitterPtrVector &
00130     FittedParticleIterator<TmplTraits>::getFitterPtrVector()
00131     {
00132       return m_fitterPtrVector;
00133     }
00134 
00135     template<typename TmplTraits>
00136     const typename FittedParticleIterator<TmplTraits>::FitterPtrVector &
00137     FittedParticleIterator<TmplTraits>::getFitterPtrVector() const
00138     {
00139       return m_fitterPtrVector;
00140     }
00141 
00142     template<typename TmplTraits>
00143     typename FittedParticleIterator<TmplTraits>::Plane
00144     FittedParticleIterator<TmplTraits>::getClosestFitPlane(
00145       const Particle &particle
00146     ) const
00147     {
00148       PlaneVector fitPlanes = getFitPlaneVector();
00149       if (fitPlanes.size() > 0) {
00150         std::sort(fitPlanes.begin(), fitPlanes.end(), PlaneComparer<Particle>(particle));
00151 
00152         return fitPlanes[0];
00153       }
00154       return Plane(Vec3(-1.0, 0, 0), Vec3(DBL_MAX/2, DBL_MAX/2, DBL_MAX/2));
00155     }
00156 
00157     template<typename TmplTraits>
00158     Vec3 FittedParticleIterator<TmplTraits>::getRandomPoint() const
00159     {
00160       return getPacker().getRandomPoint();
00161     }
00162 
00163     template<typename TmplTraits>
00164     typename FittedParticleIterator<TmplTraits>::Particle
00165     FittedParticleIterator<TmplTraits>::getCandidateParticle(
00166       const Vec3 &point
00167     )
00168     {
00169       return getPacker().getCandidateParticle(point);
00170     }
00171 
00172     template<typename TmplTraits>
00173     typename FittedParticleIterator<TmplTraits>::ParticleVector
00174     FittedParticleIterator<TmplTraits>::getClosestNeighbours(
00175       const Particle& particle,
00176       int numClosest
00177     )
00178     {
00179       return getPacker().getClosestNeighbours(particle, numClosest);
00180     }
00181 
00182     template<typename TmplTraits>
00183     typename FittedParticleIterator<TmplTraits>::Particle &
00184     FittedParticleIterator<TmplTraits>::generateNext()
00185     {
00186       m_next = Fitter::getInvalidParticle();
00187       if (m_lastFailCount < getMaxInsertionFailures())
00188       {
00189         int numFails    = 0;
00190         //int insertCount = 0;
00191         FitterPtrVector fitters = getFitterPtrVector();
00192         Particle particle       = Fitter::getInvalidParticle();
00193         Particle fitParticle    = particle;
00194         Plane plane;
00195         ParticleVector neighbours;
00196         while ((!fitParticle.isValid()) && (numFails < getMaxInsertionFailures()))
00197         {
00198           particle = getCandidateParticle(getRandomPoint());
00199           neighbours = getClosestNeighbours(particle, 4);
00200           plane = getClosestFitPlane(particle);
00201           
00202           for (
00203             typename FitterPtrVector::iterator it = getFitterPtrVector().begin();
00204             (
00205               (it != getFitterPtrVector().end())
00206               &&
00207               (!fitParticle.isValid())
00208             );
00209             it++
00210           )
00211           {
00212             fitParticle = (*it)->getFitParticle(particle, neighbours, plane);
00213           }
00214           numFails++;
00215         }
00216         m_lastFailCount = numFails;
00217         console.Info()
00218           << "FittedParticleIterator<T>::generateNext: numFails="
00219           << numFails << "\n";
00220         m_successCount += ((fitParticle.isValid()) ? 1 : 0);
00221         m_next = fitParticle;
00222       }
00223       return m_next;
00224     }
00225 
00226     template<typename TmplTraits>
00227     bool FittedParticleIterator<TmplTraits>::hasNext()
00228     {
00229       return (m_next.isValid() || generateNext().isValid());
00230     }
00231 
00232     template<typename TmplTraits>
00233     typename FittedParticleIterator<TmplTraits>::Particle
00234     FittedParticleIterator<TmplTraits>::next()
00235     {
00236       if (!(m_next.isValid()))
00237       {
00238         generateNext();
00239       }
00240       const Particle next = m_next;
00241       m_next = Fitter::getInvalidParticle();
00242       return next;
00243     }
00244 
00245     template<typename TmplTraits>
00246     void FittedParticleIterator<TmplTraits>::logInfo()
00247     {
00248       console.Info()
00249         << "BoundingBox: minPt = " << getPacker().getBBox().getMinPt()
00250         << ", maxPt = " << getPacker().getBBox().getMaxPt() << "\n"
00251         << "BoundingBox: sizes = " << getPacker().getBBox().getSizes() << "\n";
00252 
00253       for (
00254           typename FitterPtrVector::iterator it = getFitterPtrVector().begin();
00255           (it != getFitterPtrVector().end());
00256           it++      
00257       ) {
00258         console.Info() << (*(*it)).toString() << "\n";
00259       }
00260       console.Info() << "Total successful fits = " << m_successCount << "\n";
00261     }
00262 
00263 
00264 
00265 
00266 
00267 
00268 
00269 
00270 
00271 
00272 
00273 
00274 
00275 
00276 
00277 
00278 
00279 
00280 
00281 
00282 
00283 
00284 
00285 
00286 
00287 
00288   
00289     template <typename TPartGen,template <typename TTPartGen> class TCubicPackWrap>
00290     RandomBoxPacker<TPartGen,TCubicPackWrap>::RandomBoxPacker(
00291       ParticleGeneratorPtr particleGeneratorPtr,
00292       ParticlePoolPtr      particlePoolPtr,
00293       NTablePtr            nTablePtr,
00294       const BoundingBox    &bBox,
00295       const BoolVector     &periodicDimensions,
00296       double               tolerance,
00297       double               cubicPackRadius,
00298       int                  maxInsertionFailures,
00299       const PlaneVector    &fitPlaneVector
00300     )
00301      : Inherited(
00302          particleGeneratorPtr,
00303          particlePoolPtr,
00304          nTablePtr,
00305          bBox,
00306          periodicDimensions,
00307          tolerance,
00308          cubicPackRadius
00309        ),
00310        m_fitPlaneVector(fitPlaneVector),
00311        m_maxInsertionFailures(maxInsertionFailures)
00312     {
00313     }
00314 
00315     template <typename TPartGen,template <typename TTPartGen> class TCubicPackWrap>
00316     RandomBoxPacker<TPartGen,TCubicPackWrap>::RandomBoxPacker(
00317       ParticleGeneratorPtr particleGeneratorPtr,
00318       ParticlePoolPtr      particlePoolPtr,
00319       NTablePtr            nTablePtr,
00320       const BoundingBox    &bBox,
00321       const BoolVector     &periodicDimensions,
00322       double               tolerance,
00323       double               cubicPackRadius,
00324       int                  maxInsertionFailures
00325     )
00326      : Inherited(
00327          particleGeneratorPtr,
00328          particlePoolPtr,
00329          nTablePtr,
00330          bBox,
00331          periodicDimensions,
00332          tolerance,
00333          cubicPackRadius
00334        ),
00335        m_fitPlaneVector(),
00336        m_maxInsertionFailures(maxInsertionFailures)
00337     {
00338       m_fitPlaneVector = getDefaultFitPlaneVector();
00339     }
00340 
00341     template <typename TPartGen,template <typename TTPartGen> class TCubicPackWrap>
00342     typename RandomBoxPacker<TPartGen,TCubicPackWrap>::PlaneVector
00343     RandomBoxPacker<TPartGen,TCubicPackWrap>::getDefaultFitPlaneVector() const
00344     {
00345       PlaneVector fitPlaneVector;
00346       if ((!this->getPeriodicDimensions()[1])) {
00347         fitPlaneVector.push_back(
00348           Plane(Vec3(0,  1, 0), this->getBBox().getMinPt())
00349         );
00350         fitPlaneVector.push_back(
00351           Plane(Vec3(0, -1, 0), this->getBBox().getMaxPt())
00352         );
00353       }
00354       if ((!this->getPeriodicDimensions()[0])) {
00355         fitPlaneVector.push_back(
00356           Plane(Vec3( 1, 0, 0), this->getBBox().getMinPt())
00357         );
00358         fitPlaneVector.push_back(
00359           Plane(Vec3(-1, 0, 0), this->getBBox().getMaxPt())
00360         );
00361       }
00362       if (
00363           !this->is2d()
00364           &&
00365           (!this->getPeriodicDimensions()[2])
00366       ) {
00367         fitPlaneVector.push_back(
00368           Plane(Vec3(0, 0,  1), this->getBBox().getMinPt())
00369         );
00370         fitPlaneVector.push_back(
00371           Plane(Vec3(0, 0, -1), this->getBBox().getMaxPt())
00372         );
00373       }
00374       return fitPlaneVector;
00375     }
00376 
00377     template <typename TPartGen,template <typename TTPartGen> class TCubicPackWrap>
00378     RandomBoxPacker<TPartGen,TCubicPackWrap>::~RandomBoxPacker()
00379     {
00380     }
00381 
00382     template <typename TPartGen,template <typename TTPartGen> class TCubicPackWrap>
00383     double
00384     RandomBoxPacker<TPartGen,TCubicPackWrap>::getRandom(
00385       double minVal,
00386       double maxVal
00387     ) const
00388     {
00389       return minVal + (maxVal-minVal)*rng::s_zeroOneUniform();
00390     }
00391 
00392     template <typename TPartGen,template <typename TTPartGen> class TCubicPackWrap>
00393     const typename RandomBoxPacker<TPartGen,TCubicPackWrap>::PlaneVector &
00394     RandomBoxPacker<TPartGen,TCubicPackWrap>::getFitPlaneVector() const
00395     {
00396       return m_fitPlaneVector;
00397     }
00398 
00399 
00400     template <typename TPartGen,template <typename TTPartGen> class TCubicPackWrap>
00401     Vec3 RandomBoxPacker<TPartGen,TCubicPackWrap>::getRandomPoint() const
00402     {
00403       return 
00404         Vec3(
00405           getRandom(this->getBBox().getMinPt().X(), this->getBBox().getMaxPt().X()),
00406           getRandom(this->getBBox().getMinPt().Y(), this->getBBox().getMaxPt().Y()),
00407           getRandom(this->getBBox().getMinPt().Z(), this->getBBox().getMaxPt().Z())
00408         );
00409     }
00410 
00411     template <typename TPartGen,template <typename TTPartGen> class TCubicPackWrap>
00412     typename RandomBoxPacker<TPartGen,TCubicPackWrap>::ParticleVector
00413     RandomBoxPacker<TPartGen,TCubicPackWrap>::getClosestNeighbours(
00414       const Particle& particle,
00415       int numClosest
00416     )
00417     {
00418       ParticleVector neighbourVector =
00419         this->getNTable().getUniqueNeighbourVector(
00420           particle.getPos(),
00421           this->getTolerance()
00422         );
00423 
00424       if (static_cast<int>(neighbourVector.size()) < numClosest)
00425       {
00426         neighbourVector =
00427           this->getNTable().getUniqueNeighbourVector(
00428             particle.getPos(),
00429             particle.getRad() + this->getTolerance()
00430           );
00431         if (static_cast<int>(neighbourVector.size()) < numClosest)
00432         {
00433           neighbourVector =
00434             this->getNTable().getUniqueNeighbourVector(
00435               particle.getPos(),
00436               this->getParticleGenerator().getMaxFitRadius() + this->getTolerance()
00437             );
00438         }
00439       }
00440 
00441       std::sort(
00442         neighbourVector.begin(),
00443         neighbourVector.end(),
00444         ParticleComparer<Particle>(particle)
00445       );
00446 
00447       if (static_cast<int>(neighbourVector.size()) > numClosest) {
00448         neighbourVector.erase(
00449           neighbourVector.begin() + numClosest,
00450           neighbourVector.end()
00451         );
00452       }
00453 
00454       return neighbourVector;
00455     }
00456 
00457     template <typename TPartGen,template <typename TTPartGen> class TCubicPackWrap>
00458     bool RandomBoxPacker<TPartGen,TCubicPackWrap >::particleIsValid(
00459       const Particle &particle
00460     ) const
00461     {
00462       return 
00463         (
00464           this->getParticleGenerator().isValidFitRadius(particle.getRad())
00465           &&
00466           Inherited::particleFitsInBBoxWithNeighbours(particle)
00467         );
00468     }
00469 
00470     template <typename TPartGen,template <typename TTPartGen> class TCubicPackWrap>
00471     int RandomBoxPacker<TPartGen,TCubicPackWrap>::getMaxInsertionFailures() const
00472     {
00473       return m_maxInsertionFailures;
00474     }
00475 
00476     template <typename TPartGen,template <typename TTPartGen> class TCubicPackWrap>
00477     void RandomBoxPacker<TPartGen,TCubicPackWrap>::generateRandomFill()
00478     {
00479       StuffedParticleIterator it =
00480         StuffedParticleIterator(
00481           *this,
00482           getMaxInsertionFailures(),
00483           getFitPlaneVector()
00484         );
00485       while (it.hasNext())
00486       {
00487         const Particle p = it.next();
00488         createAndInsertParticle(p);
00489       }
00490       it.logInfo();
00491     }
00492 
00493     template <typename TPartGen,template <typename TTPartGen> class TCubicPackWrap>
00494     void RandomBoxPacker<TPartGen,TCubicPackWrap>::generate()
00495     {
00496       this->generateCubicPacking();
00497       this->generateRandomFill();
00498     }
00499 
00500   };
00501 };