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 #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 };