ESyS-Particle  4.0.1
GougeConfig.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 
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 }