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 00015 #include "Foundation/console.h" 00016 #include "Geometry/GeometryInfo.h" 00017 #include "Geometry/CircularNeighbourTable.h" 00018 00019 #include <stdexcept> 00020 #include <fstream> 00021 #include <sstream> 00022 #include <iomanip> 00023 00024 namespace esys 00025 { 00026 namespace lsm 00027 { 00028 //========================================================================== 00029 PackingInfo::PackingInfo( 00030 const BoundingBox &bBox, 00031 const BoolVector &periodicDimensions, 00032 Orientation orientation, 00033 double minRadius, 00034 double maxRadius 00035 ) : m_bBox(bBox), 00036 m_periodicDimensions(periodicDimensions), 00037 m_orientation(orientation), 00038 m_minRadius(minRadius), 00039 m_maxRadius(maxRadius) 00040 { 00041 initialiseFitPlaneVector(); 00042 } 00043 00044 bool PackingInfo::is3d() const 00045 { 00046 return (m_bBox.getSizes().Z() > 0.0); 00047 } 00048 00049 void PackingInfo::initialiseFitPlaneVector() 00050 { 00051 m_fitPlaneVector.clear(); 00052 if ((m_orientation != XZ) && (!getPeriodicDimensions()[1])) { 00053 m_fitPlaneVector.push_back( 00054 Plane(Vec3(0, 1, 0), getBBox().getMinPt()) 00055 ); 00056 m_fitPlaneVector.push_back( 00057 Plane(Vec3(0, -1, 0), getBBox().getMaxPt()) 00058 ); 00059 } 00060 if ((m_orientation != YZ) && (!getPeriodicDimensions()[0])) { 00061 m_fitPlaneVector.push_back( 00062 Plane(Vec3( 1, 0, 0), getBBox().getMinPt()) 00063 ); 00064 m_fitPlaneVector.push_back( 00065 Plane(Vec3(-1, 0, 0), getBBox().getMaxPt()) 00066 ); 00067 } 00068 if ( 00069 is3d() 00070 && 00071 (m_orientation != XY) 00072 && 00073 (!getPeriodicDimensions()[2]) 00074 ) { 00075 m_fitPlaneVector.push_back( 00076 Plane(Vec3(0, 0, 1), getBBox().getMinPt()) 00077 ); 00078 m_fitPlaneVector.push_back( 00079 Plane(Vec3(0, 0, -1), getBBox().getMaxPt()) 00080 ); 00081 } 00082 } 00083 00084 const BoundingBox &PackingInfo::getBBox() const 00085 { 00086 return m_bBox; 00087 } 00088 00089 const PlaneVector &PackingInfo::getFitPlaneVector() const 00090 { 00091 return m_fitPlaneVector; 00092 } 00093 00094 double PackingInfo::getMinParticleRadius() const 00095 { 00096 return m_minRadius; 00097 } 00098 00099 double PackingInfo::getMaxParticleRadius() const 00100 { 00101 return m_maxRadius; 00102 } 00103 00104 const BoolVector &PackingInfo::getPeriodicDimensions() const 00105 { 00106 return m_periodicDimensions; 00107 } 00108 00109 //========================================================================== 00110 00111 template <typename TGrainGen> 00112 GougePackingInfo<TGrainGen>::GougePackingInfo( 00113 const BoundingBox &bBox, 00114 const BoolVector &periodicDimensions, 00115 Orientation orientation, 00116 ParticleGrainGen &particleGrainGen 00117 ) : Inherited( 00118 bBox, 00119 periodicDimensions, 00120 orientation, 00121 particleGrainGen.getMinParticleRadius(), 00122 particleGrainGen.getMaxParticleRadius() 00123 ), 00124 m_pParticleGrainGen(&particleGrainGen) 00125 00126 { 00127 } 00128 00129 template <typename TGrainGen> 00130 typename GougePackingInfo<TGrainGen>::ParticleGrainGen & 00131 GougePackingInfo<TGrainGen>::getParticleGrainGen() const 00132 { 00133 return *m_pParticleGrainGen; 00134 } 00135 00136 #if 0 00137 template <typename TGrainGen> 00138 const typename GougePackingInfo<TGrainGen>::ParticleGrainGen & 00139 GougePackingInfo<TGrainGen>::getParticleGrainGen() const 00140 { 00141 return *m_pParticleGrainGen; 00142 } 00143 #endif 00144 00145 template <typename TGrainGen> 00146 double GougePackingInfo<TGrainGen>::getMinGrainRadius() const 00147 { 00148 return getParticleGrainGen().getMinGrainRadius(); 00149 } 00150 00151 template <typename TGrainGen> 00152 double GougePackingInfo<TGrainGen>::getMaxGrainRadius() const 00153 { 00154 return getParticleGrainGen().getMaxGrainRadius(); 00155 } 00156 00157 //========================================================================== 00158 00159 ParticleRndPackPrms::ParticleRndPackPrms() 00160 : m_size(0.0), 00161 m_minParticleRadius(0.0), 00162 m_maxParticleRadius(0.0) 00163 { 00164 } 00165 00166 ParticleRndPackPrms::ParticleRndPackPrms( 00167 double size, 00168 double minRadius, 00169 double maxRadius 00170 ) : m_size(size), 00171 m_minParticleRadius(minRadius), 00172 m_maxParticleRadius(maxRadius) 00173 { 00174 } 00175 00176 ParticleRndPackPrms::~ParticleRndPackPrms() 00177 { 00178 } 00179 00180 double ParticleRndPackPrms::getSize() const 00181 { 00182 return m_size; 00183 } 00184 00185 double ParticleRndPackPrms::getMinParticleRadius() const 00186 { 00187 return m_minParticleRadius; 00188 } 00189 00190 double ParticleRndPackPrms::getMaxParticleRadius() const 00191 { 00192 return m_maxParticleRadius; 00193 } 00194 00195 //========================================================================== 00196 00197 template <typename TPGrainGen> 00198 GrainRndPackPrms<TPGrainGen>::GrainRndPackPrms() 00199 : 00200 Inherited(), 00201 m_pParticleGrainGen(NULL) 00202 { 00203 } 00204 00205 template <typename TPGrainGen> 00206 GrainRndPackPrms<TPGrainGen>::GrainRndPackPrms( 00207 double size, 00208 ParticleGrainGen &particleGrainGen, 00209 int connectionTag 00210 ) : Inherited( 00211 size, 00212 particleGrainGen.getMinParticleRadius(), 00213 particleGrainGen.getMaxParticleRadius() 00214 ), 00215 m_pParticleGrainGen(&particleGrainGen), 00216 m_connectionTag(connectionTag) 00217 { 00218 } 00219 00220 template <typename TPGrainGen> 00221 typename GrainRndPackPrms<TPGrainGen>::ParticleGrainGen & 00222 GrainRndPackPrms<TPGrainGen>::getParticleGrainGen() const 00223 { 00224 return *m_pParticleGrainGen; 00225 } 00226 00227 template <typename TPGrainGen> 00228 int 00229 GrainRndPackPrms<TPGrainGen>::getConnectionTag() const 00230 { 00231 return m_connectionTag; 00232 } 00233 00234 #if 0 00235 template <typename TPGrainGen> 00236 const typename GrainRndPackPrms<TPGrainGen>::ParticleGrainGen & 00237 GrainRndPackPrms<TPGrainGen>::getParticleGrainGen() const 00238 { 00239 return *m_pParticleGrainGen; 00240 } 00241 #endif 00242 00243 template <typename TPGrainGen> 00244 double GrainRndPackPrms<TPGrainGen>::getMinGrainRadius() 00245 { 00246 return getParticleGrainGen().getMinGrainRadius(); 00247 } 00248 00249 template <typename TPGrainGen> 00250 double GrainRndPackPrms<TPGrainGen>::getMaxGrainRadius() 00251 { 00252 return getParticleGrainGen().getMaxGrainRadius(); 00253 } 00254 00255 //========================================================================== 00256 00257 template <typename TPGrainGen> 00258 GougeConfigPrms<TPGrainGen>::GougeConfigPrms() 00259 : m_bBox(Vec3::ZERO, Vec3::ZERO), 00260 m_padRadius(0.0), 00261 m_orientation(XZ), 00262 m_faultPrms(), 00263 m_gougePrms(), 00264 m_periodicDimensions(3, false), 00265 m_maxInsertionFailures(50), 00266 m_tolerance(DBL_EPSILON*128), 00267 m_connectionTolerance(DBL_EPSILON*128*10), 00268 m_blockConnectionTag(0) 00269 { 00270 } 00271 00272 template <typename TPGrainGen> 00273 GougeConfigPrms<TPGrainGen>::GougeConfigPrms( 00274 const BoundingBox &bBox, 00275 double padRadius, 00276 Orientation orientation, 00277 const ParticleRndPackPrms &faultRegionPrms, 00278 const GrainRPackPrms &gougeRegionPrms, 00279 const BoolVector &periodicDimensions, 00280 int maxInsertionFailures, 00281 double tolerance , 00282 double connectionTolerance, 00283 int blockConnectionTag 00284 ) : m_bBox(bBox), 00285 m_padRadius(padRadius), 00286 m_orientation(orientation), 00287 m_faultPrms(faultRegionPrms), 00288 m_gougePrms(gougeRegionPrms), 00289 m_periodicDimensions(periodicDimensions), 00290 m_maxInsertionFailures(maxInsertionFailures), 00291 m_tolerance(tolerance), 00292 m_connectionTolerance(connectionTolerance), 00293 m_blockConnectionTag(blockConnectionTag) 00294 { 00295 m_bBox = GridIterator(bBox, getMaxRadius()).getSphereBBox(); 00296 } 00297 00298 template <typename TPGrainGen> 00299 GougeConfigPrms<TPGrainGen>::~GougeConfigPrms() 00300 { 00301 } 00302 00303 template <typename TPGrainGen> 00304 const BoundingBox &GougeConfigPrms<TPGrainGen>::getBBox() const 00305 { 00306 return m_bBox; 00307 } 00308 00309 template <typename TPGrainGen> 00310 int GougeConfigPrms<TPGrainGen>::getGougeConnectionTag() const 00311 { 00312 return m_gougePrms.getConnectionTag(); 00313 } 00314 00315 template <typename TPGrainGen> 00316 int GougeConfigPrms<TPGrainGen>::getBlockConnectionTag() const 00317 { 00318 return m_blockConnectionTag; 00319 } 00320 00321 template <typename TPGrainGen> 00322 int GougeConfigPrms<TPGrainGen>::getMaxInsertionFailures() const 00323 { 00324 return m_maxInsertionFailures; 00325 } 00326 00327 template <typename TPGrainGen> 00328 const BoolVector &GougeConfigPrms<TPGrainGen>::getPeriodicDimensions() const 00329 { 00330 return m_periodicDimensions; 00331 } 00332 00333 template <typename TPGrainGen> 00334 Orientation GougeConfigPrms<TPGrainGen>::getOrientation() const 00335 { 00336 return m_orientation; 00337 } 00338 00339 template <typename TPGrainGen> 00340 double GougeConfigPrms<TPGrainGen>::getTolerance() const 00341 { 00342 return m_tolerance; 00343 } 00344 00345 template <typename TPGrainGen> 00346 double GougeConfigPrms<TPGrainGen>::getConnectionTolerance() const 00347 { 00348 return m_connectionTolerance; 00349 } 00350 00351 template <typename TPGrainGen> 00352 int GougeConfigPrms<TPGrainGen>::getOrientationIndex() const 00353 { 00354 int idx = 0; 00355 switch (m_orientation) { 00356 case XZ: 00357 { 00358 idx = 1; 00359 break; 00360 } 00361 case YZ: 00362 { 00363 idx = 0; 00364 break; 00365 } 00366 case XY: 00367 { 00368 idx = 2; 00369 break; 00370 } 00371 default: 00372 { 00373 std::stringstream msg; 00374 msg << "Invalid orientation: " << m_orientation; 00375 throw std::runtime_error(msg.str()); 00376 } 00377 } 00378 return idx; 00379 } 00380 00381 template <typename TPGrainGen> 00382 BoundingBox GougeConfigPrms<TPGrainGen>::cutFromCentre(double d1, double d2) const 00383 { 00384 const int idx = getOrientationIndex(); 00385 const BoundingBox bBox = getBBox(); 00386 00387 const double cmp1 = d1 + (bBox.getMaxPt()[idx] + bBox.getMinPt()[idx])/2.0; 00388 const double cmp2 = d2 + (bBox.getMaxPt()[idx] + bBox.getMinPt()[idx])/2.0; 00389 Vec3 minPt = bBox.getMinPt(); 00390 Vec3 maxPt = bBox.getMaxPt(); 00391 00392 minPt[idx] = std::min(cmp1, cmp2); 00393 maxPt[idx] = std::max(cmp1, cmp2); 00394 00395 Vec3 tmpPt = maxPt; 00396 tmpPt[idx] = minPt[idx]; 00397 00398 if ((tmpPt - bBox.getMinPt())[idx] > getTolerance()) { 00399 const BoundingBox minBBox = GridIterator(BoundingBox(bBox.getMinPt(), tmpPt), getMaxRadius()).getSphereBBox(); 00400 tmpPt = minPt; 00401 tmpPt[idx] = minBBox.getMaxPt()[idx]; 00402 } 00403 else { 00404 tmpPt = bBox.getMinPt(); 00405 } 00406 const BoundingBox maxBBox = GridIterator(BoundingBox(bBox.getMinPt(), maxPt), getMaxRadius()).getSphereBBox(); 00407 BoundingBox returnBBox = BoundingBox(tmpPt, maxBBox.getMaxPt()); 00408 00409 return returnBBox; 00410 } 00411 00412 template <typename TPGrainGen> 00413 double GougeConfigPrms<TPGrainGen>::getRegularBlockRadius() const 00414 { 00415 return m_padRadius; 00416 } 00417 00418 template <typename TPGrainGen> 00419 double GougeConfigPrms<TPGrainGen>::getFaultMinRadius() const 00420 { 00421 return m_faultPrms.getMinParticleRadius(); 00422 } 00423 00424 template <typename TPGrainGen> 00425 double GougeConfigPrms<TPGrainGen>::getFaultMaxRadius() const 00426 { 00427 return m_faultPrms.getMaxParticleRadius(); 00428 } 00429 00430 template <typename TPGrainGen> 00431 double GougeConfigPrms<TPGrainGen>::getGougeMinRadius() const 00432 { 00433 return m_gougePrms.getMinParticleRadius(); 00434 } 00435 00436 template <typename TPGrainGen> 00437 double GougeConfigPrms<TPGrainGen>::getGougeMaxRadius() const 00438 { 00439 return m_gougePrms.getMaxParticleRadius(); 00440 } 00441 00442 template <typename TPGrainGen> 00443 double GougeConfigPrms<TPGrainGen>::getOrientationSize() const 00444 { 00445 return 00446 ( 00447 getBBox().getMaxPt()[getOrientationIndex()] 00448 - 00449 getBBox().getMinPt()[getOrientationIndex()] 00450 ); 00451 } 00452 00453 template <typename TPGrainGen> 00454 BoundingBoxVector GougeConfigPrms<TPGrainGen>::getRegularBBoxVector() const 00455 { 00456 BoundingBoxVector bBoxVector; 00457 if ( 00458 (getOrientationSize() - (m_gougePrms.getSize() + 2*m_faultPrms.getSize())) 00459 > 00460 2.0*m_padRadius 00461 ) 00462 { 00463 bBoxVector.reserve(2); 00464 bBoxVector.push_back( 00465 cutFromCentre( 00466 -(m_gougePrms.getSize()/2.0 + m_faultPrms.getSize()), 00467 -getOrientationSize()/2.0 00468 ) 00469 ); 00470 bBoxVector.push_back( 00471 cutFromCentre( 00472 m_gougePrms.getSize()/2.0 + m_faultPrms.getSize(), 00473 getOrientationSize()/2.0 00474 ) 00475 ); 00476 } 00477 return bBoxVector; 00478 } 00479 00480 template <typename TPGrainGen> 00481 typename GougeConfigPrms<TPGrainGen>::GougePackingInfoVector 00482 GougeConfigPrms<TPGrainGen>::getGougePackingInfoVector() const 00483 { 00484 GougePackingInfoVector infoVec; 00485 if (m_gougePrms.getSize() > 0.0) { 00486 Vec3 overlap = Vec3::ZERO; 00487 overlap[getOrientationIndex()] = m_faultPrms.getMaxParticleRadius(); 00488 BoundingBox bBox = 00489 cutFromCentre( 00490 m_gougePrms.getSize()/2.0, 00491 -m_gougePrms.getSize()/2.0 00492 ); 00493 00494 infoVec.push_back( 00495 GougePackInfo( 00496 BoundingBox(bBox.getMinPt() - overlap, bBox.getMaxPt() + overlap), 00497 getPeriodicDimensions(), 00498 getOrientation(), 00499 m_gougePrms.getParticleGrainGen() 00500 ) 00501 ); 00502 } 00503 return infoVec; 00504 } 00505 00506 template <typename TPGrainGen> 00507 PackingInfoVector GougeConfigPrms<TPGrainGen>::getFaultPackingInfoVector() const 00508 { 00509 PackingInfoVector infoVec; 00510 if (m_faultPrms.getSize() > 0.0) 00511 { 00512 if ( 00513 (getOrientationSize() - (m_gougePrms.getSize() + 2.0*m_faultPrms.getSize())) 00514 > 00515 0.0 00516 ) 00517 { 00518 infoVec.reserve(2); 00519 const double roughnessSize = m_faultPrms.getSize(); 00520 Vec3 overlap = Vec3::ZERO; 00521 overlap[getOrientationIndex()] = m_padRadius; 00522 const BoundingBox bBox1 = 00523 cutFromCentre( 00524 -m_gougePrms.getSize()/2.0, 00525 -(m_gougePrms.getSize()/2.0 + roughnessSize) 00526 ); 00527 const BoundingBox bBox2 = 00528 cutFromCentre( 00529 m_gougePrms.getSize()/2.0, 00530 m_gougePrms.getSize()/2.0 + roughnessSize 00531 ); 00532 00533 infoVec.push_back( 00534 PackingInfo( 00535 BoundingBox(bBox1.getMinPt() - overlap, bBox1.getMaxPt()), 00536 getPeriodicDimensions(), 00537 getOrientation(), 00538 getFaultMinRadius(), 00539 getFaultMaxRadius() 00540 ) 00541 ); 00542 00543 infoVec.push_back( 00544 PackingInfo( 00545 BoundingBox(bBox2.getMinPt(), bBox2.getMaxPt() + overlap), 00546 getPeriodicDimensions(), 00547 getOrientation(), 00548 getFaultMinRadius(), 00549 getFaultMaxRadius() 00550 ) 00551 ); 00552 } 00553 else 00554 { 00555 std::stringstream msg; 00556 msg 00557 << "Roughness size plus gouge size is greater than block size: " 00558 << "2*" << m_faultPrms.getSize() << " + " << m_gougePrms.getSize() 00559 << " > " << getOrientationSize(); 00560 throw std::runtime_error(msg.str().c_str()); 00561 } 00562 } 00563 return infoVec; 00564 } 00565 00566 template <typename TPGrainGen> 00567 double GougeConfigPrms<TPGrainGen>::getMaxRadius() const 00568 { 00569 return 00570 std::max( 00571 m_padRadius, 00572 std::max(m_faultPrms.getMaxParticleRadius(), m_gougePrms.getMaxParticleRadius()) 00573 ); 00574 } 00575 00576 template <typename TPGrainGen> 00577 double GougeConfigPrms<TPGrainGen>::getMinRadius() const 00578 { 00579 return 00580 std::min( 00581 m_padRadius, 00582 std::min(m_faultPrms.getMinParticleRadius(), m_gougePrms.getMinParticleRadius()) 00583 ); 00584 } 00585 00586 template <typename TPGrainGen> 00587 bool GougeConfigPrms<TPGrainGen>::is2d() const 00588 { 00589 return (getBBox().getSizes()[2] == 0.0); 00590 } 00591 00592 //========================================================================== 00593 00594 template <typename TGPckr,typename TPPckr,typename TConn> 00595 int GougeConfig<TGPckr,TPPckr,TConn>::getNumParticles() const 00596 { 00597 int numParticles = 0; 00598 for ( 00599 typename GeneratorPtrVector::const_iterator it = m_genPtrVector.begin(); 00600 it != m_genPtrVector.end(); 00601 it++ 00602 ) 00603 { 00604 numParticles += (*it)->getNumParticles(); 00605 } 00606 return numParticles; 00607 } 00608 00609 template <typename TGPckr,typename TPPckr,typename TConn> 00610 int GougeConfig<TGPckr,TPPckr,TConn>::getNumGrains() const 00611 { 00612 int numGrains = 0; 00613 for ( 00614 typename GrainRndPackerPtrVector::const_iterator packerIt = 00615 getGougeGeneratorVector().begin(); 00616 packerIt != getGougeGeneratorVector().end(); 00617 packerIt++ 00618 ) 00619 { 00620 numGrains += (*packerIt)->getNumGrains(); 00621 } 00622 return numGrains; 00623 } 00624 00625 template <typename TGPckr,typename TPPckr,typename TConn> 00626 int GougeConfig<TGPckr,TPPckr,TConn>::getNumConnections() const 00627 { 00628 return m_connectionSet.size(); 00629 } 00630 00631 template <typename TGPckr,typename TPPckr,typename TConn> 00632 GougeConfig<TGPckr,TPPckr,TConn>::GougeConfig(const GougeConfPrms &prms) 00633 : m_nTablePtr(), 00634 m_prms(prms), 00635 m_connectionSet(), 00636 m_gougeGenPtrVector(), 00637 m_genPtrVector(), 00638 m_particlePoolPtr(new ParticlePool(8*4096)), 00639 m_grainPoolPtr(new GrainPool(4096)), 00640 m_regularGenPtrVector(), 00641 m_faultGenPtrVector() 00642 { 00643 /* 00644 * Adjust the size of the ntable bounding-box to accommodate circular 00645 * boundary conditions. 00646 */ 00647 const BoundingBox bBox = m_prms.getBBox(); 00648 Vec3 ntableAdjust = 00649 Vec3( 00650 m_prms.getPeriodicDimensions()[0] ? 1 : 0, 00651 m_prms.getPeriodicDimensions()[1] ? 1 : 0, 00652 m_prms.getPeriodicDimensions()[2] ? 1 : 0 00653 )*m_prms.getMaxRadius(); 00654 if (m_prms.getBBox().getSizes().Z() >= 4*m_prms.getMaxRadius()) { 00655 ntableAdjust += Vec3(m_prms.getMaxRadius(), 0, 0); 00656 } 00657 const BoundingBox nTableBBox(bBox.getMinPt(), bBox.getMaxPt() - ntableAdjust); 00658 m_nTablePtr = 00659 NTablePtr( 00660 new NTable( 00661 nTableBBox, 00662 (4.0*m_prms.getMinRadius()), // grid spacing 00663 m_prms.getPeriodicDimensions(), 00664 2.1*m_prms.getMaxRadius() // width of border-region in which 00665 // particles are duplicated 00666 // for circular boundry 00667 ) 00668 ); 00669 } 00670 00671 template <typename TGPckr,typename TPPckr,typename TConn> 00672 void GougeConfig<TGPckr,TPPckr,TConn>::createRegularBlockGenerators() 00673 { 00674 BoundingBoxVector bBoxVector = m_prms.getRegularBBoxVector(); 00675 for ( 00676 BoundingBoxVector::const_iterator it = bBoxVector.begin(); 00677 it != bBoxVector.end(); 00678 it++ 00679 ) { 00680 console.Debug() 00681 << "GougeConfig<TGPckr,TPPckr,TConn>::createRegularBlockGenerators:" 00682 << "Creating RegBoxPacker in box: " << StringUtil::toString(*it) 00683 << "\n"; 00684 00685 GeneratorPtr genPtr = 00686 GeneratorPtr( 00687 new RegBoxPacker( 00688 RegRadiusGenPtr(new RegRadiusGen(m_prms.getRegularBlockRadius())), 00689 m_particlePoolPtr, 00690 m_nTablePtr, 00691 *it, 00692 m_prms.getPeriodicDimensions(), 00693 m_prms.getTolerance(), 00694 m_prms.getRegularBlockRadius() 00695 ) 00696 ); 00697 m_genPtrVector.push_back(genPtr); 00698 m_regularGenPtrVector.push_back(genPtr); 00699 } 00700 } 00701 00702 template <typename TGPckr,typename TPPckr,typename TConn> 00703 void GougeConfig<TGPckr,TPPckr,TConn>::createFaultBlockGenerators() 00704 { 00705 PackingInfoVector infoVec = m_prms.getFaultPackingInfoVector(); 00706 for ( 00707 PackingInfoVector::const_iterator it = infoVec.begin(); 00708 it != infoVec.end(); 00709 it++ 00710 ) { 00711 console.Debug() 00712 << "GougeConfig<TGPckr,TPPckr,TConn>::createFaultBlockGenerators:" 00713 << "Creating RndBoxPacker in box: " << StringUtil::toString(it->getBBox()) 00714 << "\n"; 00715 GeneratorPtr genPtr = 00716 GeneratorPtr( 00717 new RndBoxPacker( 00718 RndRadiusGenPtr( 00719 new RndRadiusGen( 00720 it->getMinParticleRadius(), 00721 it->getMaxParticleRadius() 00722 ) 00723 ), 00724 m_particlePoolPtr, 00725 m_nTablePtr, 00726 it->getBBox(), 00727 it->getPeriodicDimensions(), 00728 m_prms.getTolerance(), 00729 it->getMaxParticleRadius(), 00730 m_prms.getMaxInsertionFailures(), 00731 it->getFitPlaneVector() 00732 ) 00733 ); 00734 m_genPtrVector.push_back(genPtr); 00735 m_faultGenPtrVector.push_back(genPtr); 00736 } 00737 } 00738 00739 template <typename TGPckr,typename TPPckr,typename TConn> 00740 void GougeConfig<TGPckr,TPPckr,TConn>::createGougeConfigGenerators() 00741 { 00742 GougePackingInfoVector infoVec = m_prms.getGougePackingInfoVector(); 00743 for ( 00744 typename GougePackingInfoVector::const_iterator it = infoVec.begin(); 00745 it != infoVec.end(); 00746 it++ 00747 ) { 00748 console.Debug() 00749 << "GougeConfig<TGPckr,TPPckr,TConn>::createGougeConfigGenerators:" 00750 << "Creating GrainRandomPacker in box: " << StringUtil::toString(it->getBBox()) 00751 << "\n"; 00752 GrainRandomPackerPtr genPtr = 00753 GrainRandomPackerPtr( 00754 new GrainRandomPacker( 00755 typename GrainRandomPacker::ParticleGrainGenPtr(), 00756 m_particlePoolPtr, 00757 m_nTablePtr, 00758 it->getBBox(), 00759 it->getPeriodicDimensions(), 00760 m_prms.getTolerance(), 00761 it->getParticleGrainGen().getMaxGrainRadius(), 00762 m_prms.getMaxInsertionFailures(), 00763 it->getFitPlaneVector(), 00764 m_grainPoolPtr 00765 ) 00766 ); 00767 genPtr->setParticleGrainGen(it->getParticleGrainGen()); 00768 m_genPtrVector.push_back(genPtr); 00769 m_gougeGenPtrVector.push_back(genPtr); 00770 } 00771 } 00772 00773 template <typename TGPckr,typename TPPckr,typename TConn> 00774 GougeConfig<TGPckr,TPPckr,TConn>::~GougeConfig() 00775 { 00776 } 00777 00778 template <typename TGPckr,typename TPPckr,typename TConn> 00779 void GougeConfig<TGPckr,TPPckr,TConn>::generate() 00780 { 00781 // setup block generators 00782 createRegularBlockGenerators(); 00783 createFaultBlockGenerators(); 00784 createGougeConfigGenerators(); 00785 00786 // use block generators 00787 console.Info() << "bbox = " << m_prms.getBBox().getMinPt() << " " << m_prms.getBBox().getMaxPt() << "\n"; 00788 for ( 00789 typename GeneratorPtrVector::iterator it = m_genPtrVector.begin(); 00790 it != m_genPtrVector.end(); 00791 it++ 00792 ) 00793 { 00794 (*it)->generate(); 00795 } 00796 createConnectionSet(); 00797 } 00798 00799 template <typename TGPckr,typename TPPckr,typename TConn> 00800 void GougeConfig<TGPckr,TPPckr,TConn>::writeToFile(const std::string &fileName) const 00801 { 00802 std::ofstream fStream(fileName.c_str()); 00803 write(fStream); 00804 } 00805 00806 template <typename TGPckr,typename TPPckr,typename TConn> 00807 const typename GougeConfig<TGPckr,TPPckr,TConn>::GrainRndPackerPtrVector & 00808 GougeConfig<TGPckr,TPPckr,TConn>::getGougeGeneratorVector() const 00809 { 00810 return m_gougeGenPtrVector; 00811 } 00812 00813 template <typename TGPckr,typename TPPckr,typename TConn> 00814 typename GougeConfig<TGPckr,TPPckr,TConn>::GrainRndPackerPtrVector & 00815 GougeConfig<TGPckr,TPPckr,TConn>::getGougeGeneratorVector() 00816 { 00817 return m_gougeGenPtrVector; 00818 } 00819 00820 template <typename TGPckr,typename TPPckr,typename TConn> 00821 const typename GougeConfig<TGPckr,TPPckr,TConn>::GeneratorPtrVector & 00822 GougeConfig<TGPckr,TPPckr,TConn>::getFaultGeneratorVector() const 00823 { 00824 return m_faultGenPtrVector; 00825 } 00826 00827 template <typename TGPckr,typename TPPckr,typename TConn> 00828 bool GougeConfig<TGPckr,TPPckr,TConn>::areInDifferentFaultBlocks( 00829 const Particle &p1, 00830 const Particle &p2 00831 ) const 00832 { 00833 const GeneratorPtrVector &generators = getFaultGeneratorVector(); 00834 if (generators.size() == 2) { 00835 return 00836 ( 00837 (generators[0]->contains(p1) && generators[1]->contains(p2)) 00838 || 00839 (generators[0]->contains(p2) && generators[1]->contains(p1)) 00840 ); 00841 } 00842 else if (generators.size() > 2) { 00843 throw 00844 std::runtime_error( 00845 "GougeConfig<TGPckr,TPPckr,TConn>::areInDifferentFaultBlocks: " 00846 "More than two fault blocks." 00847 ); 00848 } 00849 return false; 00850 } 00851 00852 template <typename TGPckr,typename TPPckr,typename TConn> 00853 bool GougeConfig<TGPckr,TPPckr,TConn>::isGougeParticle(const Particle &particle) const 00854 { 00855 const GrainRndPackerPtrVector &generators = getGougeGeneratorVector(); 00856 for ( 00857 typename GrainRndPackerPtrVector::const_iterator it = generators.begin(); 00858 it != generators.end(); 00859 it++ 00860 ) 00861 { 00862 if ((*it)->contains(particle)) { 00863 return true; 00864 } 00865 } 00866 return false; 00867 } 00868 00869 template <typename TGPckr,typename TPPckr,typename TConn> 00870 void GougeConfig<TGPckr,TPPckr,TConn>::createConnectionSet() 00871 { 00872 // 00873 // First created connections in the elastic blocks. 00874 // 00875 typename NTable::ParticleIterator particleIt = m_nTablePtr->getParticleIterator(); 00876 ConnectionValidator validator = ConnectionValidator(*this, m_prms.getConnectionTolerance()); 00877 while (particleIt.hasNext()) { 00878 const typename NTable::Particle *pParticle = particleIt.next(); 00879 const typename NTable::ParticleVector neighbours = 00880 m_nTablePtr->getNeighbourVector( 00881 pParticle->getPos(), 00882 pParticle->getRad() + m_prms.getConnectionTolerance() 00883 ); 00884 for ( 00885 typename NTable::ParticleVector::const_iterator it = neighbours.begin(); 00886 it != neighbours.end(); 00887 it++ 00888 ) 00889 { 00890 if (validator.isValid(*pParticle, *(*it))) { 00891 m_connectionSet.insert( 00892 typename ConnectionSet::value_type( 00893 pParticle->getID(), 00894 (*it)->getID(), 00895 m_prms.getBlockConnectionTag() 00896 ) 00897 ); 00898 } 00899 } 00900 } 00901 const int numBlockConns = m_connectionSet.size(); 00902 console.Info() 00903 << "Created " << numBlockConns << " connections in " 00904 << "bonded blocks.\n"; 00905 00906 // 00907 // Create connections with grains. 00908 // 00909 console.Debug() 00910 << "Prms BBox: " << StringUtil::toString(m_prms.getBBox()) << "\n"; 00911 console.Debug() 00912 << "NTbl BBox: " << StringUtil::toString(m_nTablePtr->getBBox()) << "\n"; 00913 for ( 00914 typename GrainRndPackerPtrVector::iterator packerIt = getGougeGeneratorVector().begin(); 00915 packerIt != getGougeGeneratorVector().end(); 00916 packerIt++ 00917 ) 00918 { 00919 typename GrainRandomPacker::GrainIterator grainIt = (*packerIt)->getGrainIterator(); 00920 while (grainIt.hasNext()) 00921 { 00922 ConnectionFinder connFinder = 00923 ConnectionFinder( 00924 m_prms.getConnectionTolerance(), 00925 m_prms.getGougeConnectionTag(), 00926 m_nTablePtr->getBBox(), 00927 m_nTablePtr->getPeriodicDimensions() 00928 ); 00929 Grain &g = grainIt.next(); 00930 connFinder.create(g.getParticleIterator()); 00931 typename ConnectionFinder::Iterator connIt = connFinder.getIterator(); 00932 while (connIt.hasNext()) 00933 { 00934 m_connectionSet.insert(connIt.next()); 00935 } 00936 if (connFinder.getNumConnections() == 0) 00937 { 00938 console.Info() 00939 << "Found no connections in grain " << g.getId() 00940 << ":\n"; 00941 typename Grain::ParticleIterator partIt = g.getParticleIterator(); 00942 while (partIt.hasNext()) 00943 { 00944 console.Info() << StringUtil::toString(partIt.next()) << "\n"; 00945 } 00946 } 00947 } 00948 } 00949 console.Info() 00950 << "Created " << m_connectionSet.size()-numBlockConns << " connections in " 00951 << "gouge region.\n"; 00952 } 00953 00954 template <typename TGPckr,typename TPPckr,typename TConn> 00955 typename GougeConfig<TGPckr,TPPckr,TConn>::ParticleCollection 00956 GougeConfig<TGPckr,TPPckr,TConn>::getParticleCollection() 00957 { 00958 ParticleCollection pCollection(m_particlePoolPtr); 00959 for ( 00960 typename GeneratorPtrVector::iterator it = m_genPtrVector.begin(); 00961 it != m_genPtrVector.end(); 00962 it++ 00963 ) 00964 { 00965 ParticleIterator particleIt = (*it)->getParticleIterator(); 00966 while (particleIt.hasNext()) { 00967 pCollection.insertRef(particleIt.next()); 00968 } 00969 } 00970 00971 return pCollection; 00972 } 00973 00974 template <typename TGPckr,typename TPPckr,typename TConn> 00975 typename GougeConfig<TGPckr,TPPckr,TConn>::GrainCollection 00976 GougeConfig<TGPckr,TPPckr,TConn>::getGrainCollection() 00977 { 00978 GrainCollection gCollection(m_particlePoolPtr, m_grainPoolPtr); 00979 for ( 00980 typename GrainRndPackerPtrVector::iterator packerIt = 00981 getGougeGeneratorVector().begin(); 00982 packerIt != getGougeGeneratorVector().end(); 00983 packerIt++ 00984 ) 00985 { 00986 GrainIterator grainIt = (*packerIt)->getGrainIterator(); 00987 while (grainIt.hasNext()) { 00988 gCollection.insertRef(grainIt.next()); 00989 } 00990 } 00991 00992 return gCollection; 00993 } 00994 00995 template <typename TGPckr,typename TPPckr,typename TConn> 00996 const typename GougeConfig<TGPckr,TPPckr,TConn>::ConnectionSet & 00997 GougeConfig<TGPckr,TPPckr,TConn>::getConnectionSet() const 00998 { 00999 return m_connectionSet; 01000 } 01001 01002 template <typename TGPckr,typename TPPckr,typename TConn> 01003 void GougeConfig<TGPckr,TPPckr,TConn>::write(std::ostream &oStream) const 01004 { 01005 Vec3 minPt = m_nTablePtr->getBBox().getMinPt(); 01006 Vec3 maxPt = m_nTablePtr->getBBox().getMaxPt(); 01007 if (fabs(maxPt.Z() - minPt.Z()) < (2*m_prms.getMaxRadius())) { 01008 minPt.Z() = minPt.Z() - m_prms.getMaxRadius() - m_prms.getTolerance(); 01009 maxPt.Z() = maxPt.Z() + m_prms.getMaxRadius() + m_prms.getTolerance(); 01010 } 01011 01012 const BoundingBox geoBBox(minPt, maxPt + m_prms.getTolerance()); 01013 01014 GeometryInfo info = 01015 GeometryInfo( 01016 1.2, 01017 geoBBox.getMinPt(), 01018 geoBBox.getMaxPt(), 01019 m_prms.getPeriodicDimensions(), 01020 (m_prms.getBBox().getSizes().Z() <= 0.0) 01021 ); 01022 info.write(oStream); 01023 01024 /* 01025 * Some ugliness to ensure that duplicated particles (because of 01026 * circular boundary conditions) get the correct tag. First, create 01027 * a set of particles which have already been tagged. 01028 */ 01029 typename NTable::ParticleIterator particleIt = m_nTablePtr->getParticleIterator(); 01030 typedef std::set<Particle *, IdCompare> ParticleSet; 01031 ParticleSet taggedParticleSet; 01032 while (particleIt.hasNext()) { 01033 taggedParticleSet.insert(particleIt.next()); 01034 } 01035 01036 /* 01037 * Now eliminate any particles which were generated with their centre-point 01038 * lying outside the geo bounding box. Also set the tag in case the particle 01039 * was a duplicate, created due to circular boundary. 01040 */ 01041 particleIt = m_nTablePtr->getParticleIterator(); 01042 ParticleSet particleSet; 01043 while (particleIt.hasNext()) { 01044 Particle *pParticle = particleIt.next(); 01045 if (geoBBox.contains(pParticle->getPos())) { 01046 pParticle->setTag((*(taggedParticleSet.find(pParticle)))->getTag()); 01047 particleSet.insert(pParticle); 01048 } 01049 } 01050 01051 /* 01052 * Write particles to the stream. 01053 */ 01054 oStream 01055 << "\n" 01056 << "BeginParticles" 01057 << "\n" 01058 << "Simple" 01059 << "\n" 01060 << particleSet.size() 01061 << "\n"; 01062 const int precision = 12; 01063 GeoParticleWriter particleVisitor(oStream, precision); 01064 for ( 01065 typename ParticleSet::const_iterator it = particleSet.begin(); 01066 it != particleSet.end(); 01067 it++ 01068 ) 01069 { 01070 particleVisitor.visitParticle(*(*it)); 01071 } 01072 01073 oStream << "EndParticles\n" << "BeginConnect\n"; 01074 oStream << getConnectionSet().size() << "\n"; 01075 oStream.flush(); 01076 01077 GeoConnectionWriter connectionVisitor(oStream); 01078 visitConnections(connectionVisitor); 01079 oStream << "EndConnect"; 01080 oStream.flush(); 01081 } 01082 01083 template <typename TGPckr,typename TPPckr,typename TConn> 01084 void GougeConfig<TGPckr,TPPckr,TConn>::tagGougeParticles(int tag) 01085 { 01086 for ( 01087 typename GrainRndPackerPtrVector::iterator it = m_gougeGenPtrVector.begin(); 01088 it != m_gougeGenPtrVector.end(); 01089 it++ 01090 ) 01091 { 01092 typename GrainRandomPacker::ParticleIterator particleIt = 01093 (*it)->getParticleIterator(); 01094 while (particleIt.hasNext()) { 01095 particleIt.next().setTag(tag); 01096 } 01097 } 01098 } 01099 01100 template <typename TGPckr,typename TPPckr,typename TConn> 01101 void GougeConfig<TGPckr,TPPckr,TConn>::tagRndBlockParticles(int tag) 01102 { 01103 for ( 01104 typename GeneratorPtrVector::iterator it = m_faultGenPtrVector.begin(); 01105 it != m_faultGenPtrVector.end(); 01106 it++ 01107 ) 01108 { 01109 ParticleIterator particleIt = (*it)->getParticleIterator(); 01110 while (particleIt.hasNext()) { 01111 particleIt.next().setTag(tag); 01112 } 01113 } 01114 } 01115 01116 template <typename TGPckr,typename TPPckr,typename TConn> 01117 void GougeConfig<TGPckr,TPPckr,TConn>::tagDrivingPlateParticles( 01118 int lowDrivingTag, 01119 int highDrivingTag, 01120 double distanceFromBBoxEdge 01121 ) 01122 { 01123 ParticleCollection particleCollection = getParticleCollection(); 01124 const BoundingBox bBox = particleCollection.getParticleBBox(); 01125 01126 const int idx = this->m_prms.getOrientationIndex(); 01127 const double maxLow = bBox.getMinPt()[idx] + distanceFromBBoxEdge; 01128 const double minHigh = bBox.getMaxPt()[idx] - distanceFromBBoxEdge; 01129 01130 int lowTagCount = 0; 01131 int highTagCount = 0; 01132 typename ParticleCollection::ParticleIterator particleIt = 01133 particleCollection.getParticleIterator(); 01134 while (particleIt.hasNext()) 01135 { 01136 Particle &particle = particleIt.next(); 01137 const double dimPos = particle.getPos()[idx]; 01138 const double radius = particle.getRad(); 01139 01140 if (dimPos - radius <= maxLow) { 01141 particle.setTag(lowDrivingTag); 01142 lowTagCount++; 01143 } 01144 if (dimPos + radius >= minHigh) { 01145 particle.setTag(highDrivingTag); 01146 highTagCount++; 01147 } 01148 } 01149 console.Info() << "Tagged " << lowTagCount << " particles with " << lowDrivingTag << "\n"; 01150 console.Info() << "Tagged " << highTagCount << " particles with " << highDrivingTag << "\n"; 01151 } 01152 } 01153 }