ESyS-Particle  2.3.2
GougeConfig.hpp
Go to the documentation of this file.
1 
2 // //
3 // Copyright (c) 2003-2014 by The University of Queensland //
4 // Centre for Geoscience Computing //
5 // http://earth.uq.edu.au/centre-geoscience-computing //
6 // //
7 // Primary Business: Brisbane, Queensland, Australia //
8 // Licensed under the Open Software License version 3.0 //
9 // http://www.opensource.org/licenses/osl-3.0.php //
10 // //
12 
13 
14 
15 #include "Foundation/console.h"
16 #include "Geometry/GeometryInfo.h"
18 
19 #include <stdexcept>
20 #include <fstream>
21 #include <sstream>
22 #include <iomanip>
23 
24 namespace esys
25 {
26  namespace lsm
27  {
28  //==========================================================================
30  const BoundingBox &bBox,
31  const BoolVector &periodicDimensions,
32  Orientation orientation,
33  double minRadius,
34  double maxRadius
35  ) : m_bBox(bBox),
36  m_periodicDimensions(periodicDimensions),
37  m_orientation(orientation),
38  m_minRadius(minRadius),
39  m_maxRadius(maxRadius)
40  {
41  initialiseFitPlaneVector();
42  }
43 
44  bool PackingInfo::is3d() const
45  {
46  return (m_bBox.getSizes().Z() > 0.0);
47  }
48 
49  void PackingInfo::initialiseFitPlaneVector()
50  {
51  m_fitPlaneVector.clear();
52  if ((m_orientation != XZ) && (!getPeriodicDimensions()[1])) {
53  m_fitPlaneVector.push_back(
54  Plane3D(Vec3(0, 1, 0), getBBox().getMinPt())
55  );
56  m_fitPlaneVector.push_back(
57  Plane3D(Vec3(0, -1, 0), getBBox().getMaxPt())
58  );
59  }
60  if ((m_orientation != YZ) && (!getPeriodicDimensions()[0])) {
61  m_fitPlaneVector.push_back(
62  Plane3D(Vec3( 1, 0, 0), getBBox().getMinPt())
63  );
64  m_fitPlaneVector.push_back(
65  Plane3D(Vec3(-1, 0, 0), getBBox().getMaxPt())
66  );
67  }
68  if (
69  is3d()
70  &&
71  (m_orientation != XY)
72  &&
73  (!getPeriodicDimensions()[2])
74  ) {
75  m_fitPlaneVector.push_back(
76  Plane3D(Vec3(0, 0, 1), getBBox().getMinPt())
77  );
78  m_fitPlaneVector.push_back(
79  Plane3D(Vec3(0, 0, -1), getBBox().getMaxPt())
80  );
81  }
82  }
83 
84  const BoundingBox &PackingInfo::getBBox() const
85  {
86  return m_bBox;
87  }
88 
89  const PlaneVector &PackingInfo::getFitPlaneVector() const
90  {
91  return m_fitPlaneVector;
92  }
93 
94  double PackingInfo::getMinParticleRadius() const
95  {
96  return m_minRadius;
97  }
98 
99  double PackingInfo::getMaxParticleRadius() const
100  {
101  return m_maxRadius;
102  }
103 
104  const BoolVector &PackingInfo::getPeriodicDimensions() const
105  {
106  return m_periodicDimensions;
107  }
108 
109  //==========================================================================
110 
111  template <typename TGrainGen>
113  const BoundingBox &bBox,
114  const BoolVector &periodicDimensions,
115  Orientation orientation,
116  ParticleGrainGen &particleGrainGen
117  ) : Inherited(
118  bBox,
119  periodicDimensions,
120  orientation,
121  particleGrainGen.getMinParticleRadius(),
122  particleGrainGen.getMaxParticleRadius()
123  ),
124  m_pParticleGrainGen(&particleGrainGen)
125 
126  {
127  }
128 
129  template <typename TGrainGen>
132  {
133  return *m_pParticleGrainGen;
134  }
135 
136 #if 0
137  template <typename TGrainGen>
140  {
141  return *m_pParticleGrainGen;
142  }
143 #endif
144 
145  template <typename TGrainGen>
147  {
148  return getParticleGrainGen().getMinGrainRadius();
149  }
150 
151  template <typename TGrainGen>
153  {
154  return getParticleGrainGen().getMaxGrainRadius();
155  }
156 
157  //==========================================================================
158 
160  : m_size(0.0),
161  m_minParticleRadius(0.0),
162  m_maxParticleRadius(0.0)
163  {
164  }
165 
167  double size,
168  double minRadius,
169  double maxRadius
170  ) : m_size(size),
171  m_minParticleRadius(minRadius),
172  m_maxParticleRadius(maxRadius)
173  {
174  }
175 
177  {
178  }
179 
181  {
182  return m_size;
183  }
184 
186  {
187  return m_minParticleRadius;
188  }
189 
191  {
192  return m_maxParticleRadius;
193  }
194 
195  //==========================================================================
196 
197  template <typename TPGrainGen>
199  :
200  Inherited(),
201  m_pParticleGrainGen(NULL)
202  {
203  }
204 
205  template <typename TPGrainGen>
207  double size,
208  ParticleGrainGen &particleGrainGen,
209  int connectionTag
210  ) : Inherited(
211  size,
212  particleGrainGen.getMinParticleRadius(),
213  particleGrainGen.getMaxParticleRadius()
214  ),
215  m_pParticleGrainGen(&particleGrainGen),
216  m_connectionTag(connectionTag)
217  {
218  }
219 
220  template <typename TPGrainGen>
223  {
224  return *m_pParticleGrainGen;
225  }
226 
227  template <typename TPGrainGen>
228  int
230  {
231  return m_connectionTag;
232  }
233 
234 #if 0
235  template <typename TPGrainGen>
238  {
239  return *m_pParticleGrainGen;
240  }
241 #endif
242 
243  template <typename TPGrainGen>
245  {
246  return getParticleGrainGen().getMinGrainRadius();
247  }
248 
249  template <typename TPGrainGen>
251  {
252  return getParticleGrainGen().getMaxGrainRadius();
253  }
254 
255  //==========================================================================
256 
257  template <typename TPGrainGen>
259  : m_bBox(Vec3::ZERO, Vec3::ZERO),
260  m_padRadius(0.0),
261  m_orientation(XZ),
262  m_faultPrms(),
263  m_gougePrms(),
264  m_periodicDimensions(3, false),
265  m_maxInsertionFailures(50),
266  m_tolerance(DBL_EPSILON*128),
267  m_connectionTolerance(DBL_EPSILON*128*10),
268  m_blockConnectionTag(0)
269  {
270  }
271 
272  template <typename TPGrainGen>
274  const BoundingBox &bBox,
275  double padRadius,
276  Orientation orientation,
277  const ParticleRndPackPrms &faultRegionPrms,
278  const GrainRPackPrms &gougeRegionPrms,
279  const BoolVector &periodicDimensions,
280  int maxInsertionFailures,
281  double tolerance ,
282  double connectionTolerance,
283  int blockConnectionTag
284  ) : m_bBox(bBox),
285  m_padRadius(padRadius),
286  m_orientation(orientation),
287  m_faultPrms(faultRegionPrms),
288  m_gougePrms(gougeRegionPrms),
289  m_periodicDimensions(periodicDimensions),
290  m_maxInsertionFailures(maxInsertionFailures),
291  m_tolerance(tolerance),
292  m_connectionTolerance(connectionTolerance),
293  m_blockConnectionTag(blockConnectionTag)
294  {
296  }
297 
298  template <typename TPGrainGen>
300  {
301  }
302 
303  template <typename TPGrainGen>
305  {
306  return m_bBox;
307  }
308 
309  template <typename TPGrainGen>
311  {
312  return m_gougePrms.getConnectionTag();
313  }
314 
315  template <typename TPGrainGen>
317  {
318  return m_blockConnectionTag;
319  }
320 
321  template <typename TPGrainGen>
323  {
324  return m_maxInsertionFailures;
325  }
326 
327  template <typename TPGrainGen>
329  {
330  return m_periodicDimensions;
331  }
332 
333  template <typename TPGrainGen>
335  {
336  return m_orientation;
337  }
338 
339  template <typename TPGrainGen>
341  {
342  return m_tolerance;
343  }
344 
345  template <typename TPGrainGen>
347  {
348  return m_connectionTolerance;
349  }
350 
351  template <typename TPGrainGen>
353  {
354  int idx = 0;
355  switch (m_orientation) {
356  case XZ:
357  {
358  idx = 1;
359  break;
360  }
361  case YZ:
362  {
363  idx = 0;
364  break;
365  }
366  case XY:
367  {
368  idx = 2;
369  break;
370  }
371  default:
372  {
373  std::stringstream msg;
374  msg << "Invalid orientation: " << m_orientation;
375  throw std::runtime_error(msg.str());
376  }
377  }
378  return idx;
379  }
380 
381  template <typename TPGrainGen>
383  {
384  const int idx = getOrientationIndex();
385  const BoundingBox bBox = getBBox();
386 
387  const double cmp1 = d1 + (bBox.getMaxPt()[idx] + bBox.getMinPt()[idx])/2.0;
388  const double cmp2 = d2 + (bBox.getMaxPt()[idx] + bBox.getMinPt()[idx])/2.0;
389  Vec3 minPt = bBox.getMinPt();
390  Vec3 maxPt = bBox.getMaxPt();
391 
392  minPt[idx] = std::min(cmp1, cmp2);
393  maxPt[idx] = std::max(cmp1, cmp2);
394 
395  Vec3 tmpPt = maxPt;
396  tmpPt[idx] = minPt[idx];
397 
398  if ((tmpPt - bBox.getMinPt())[idx] > getTolerance()) {
399  const BoundingBox minBBox = GridIterator(BoundingBox(bBox.getMinPt(), tmpPt), getMaxRadius()).getSphereBBox();
400  tmpPt = minPt;
401  tmpPt[idx] = minBBox.getMaxPt()[idx];
402  }
403  else {
404  tmpPt = bBox.getMinPt();
405  }
406  const BoundingBox maxBBox = GridIterator(BoundingBox(bBox.getMinPt(), maxPt), getMaxRadius()).getSphereBBox();
407  BoundingBox returnBBox = BoundingBox(tmpPt, maxBBox.getMaxPt());
408 
409  return returnBBox;
410  }
411 
412  template <typename TPGrainGen>
414  {
415  return m_padRadius;
416  }
417 
418  template <typename TPGrainGen>
420  {
421  return m_faultPrms.getMinParticleRadius();
422  }
423 
424  template <typename TPGrainGen>
426  {
427  return m_faultPrms.getMaxParticleRadius();
428  }
429 
430  template <typename TPGrainGen>
432  {
433  return m_gougePrms.getMinParticleRadius();
434  }
435 
436  template <typename TPGrainGen>
438  {
439  return m_gougePrms.getMaxParticleRadius();
440  }
441 
442  template <typename TPGrainGen>
444  {
445  return
446  (
447  getBBox().getMaxPt()[getOrientationIndex()]
448  -
449  getBBox().getMinPt()[getOrientationIndex()]
450  );
451  }
452 
453  template <typename TPGrainGen>
455  {
456  BoundingBoxVector bBoxVector;
457  if (
458  (getOrientationSize() - (m_gougePrms.getSize() + 2*m_faultPrms.getSize()))
459  >
460  2.0*m_padRadius
461  )
462  {
463  bBoxVector.reserve(2);
464  bBoxVector.push_back(
465  cutFromCentre(
466  -(m_gougePrms.getSize()/2.0 + m_faultPrms.getSize()),
467  -getOrientationSize()/2.0
468  )
469  );
470  bBoxVector.push_back(
471  cutFromCentre(
472  m_gougePrms.getSize()/2.0 + m_faultPrms.getSize(),
473  getOrientationSize()/2.0
474  )
475  );
476  }
477  return bBoxVector;
478  }
479 
480  template <typename TPGrainGen>
483  {
484  GougePackingInfoVector infoVec;
485  if (m_gougePrms.getSize() > 0.0) {
486  Vec3 overlap = Vec3::ZERO;
487  overlap[getOrientationIndex()] = m_faultPrms.getMaxParticleRadius();
488  BoundingBox bBox =
489  cutFromCentre(
490  m_gougePrms.getSize()/2.0,
491  -m_gougePrms.getSize()/2.0
492  );
493 
494  infoVec.push_back(
496  BoundingBox(bBox.getMinPt() - overlap, bBox.getMaxPt() + overlap),
497  getPeriodicDimensions(),
498  getOrientation(),
499  m_gougePrms.getParticleGrainGen()
500  )
501  );
502  }
503  return infoVec;
504  }
505 
506  template <typename TPGrainGen>
508  {
509  PackingInfoVector infoVec;
510  if (m_faultPrms.getSize() > 0.0)
511  {
512  if (
513  (getOrientationSize() - (m_gougePrms.getSize() + 2.0*m_faultPrms.getSize()))
514  >
515  0.0
516  )
517  {
518  infoVec.reserve(2);
519  const double roughnessSize = m_faultPrms.getSize();
520  Vec3 overlap = Vec3::ZERO;
521  overlap[getOrientationIndex()] = m_padRadius;
522  const BoundingBox bBox1 =
523  cutFromCentre(
524  -m_gougePrms.getSize()/2.0,
525  -(m_gougePrms.getSize()/2.0 + roughnessSize)
526  );
527  const BoundingBox bBox2 =
528  cutFromCentre(
529  m_gougePrms.getSize()/2.0,
530  m_gougePrms.getSize()/2.0 + roughnessSize
531  );
532 
533  infoVec.push_back(
534  PackingInfo(
535  BoundingBox(bBox1.getMinPt() - overlap, bBox1.getMaxPt()),
536  getPeriodicDimensions(),
537  getOrientation(),
538  getFaultMinRadius(),
539  getFaultMaxRadius()
540  )
541  );
542 
543  infoVec.push_back(
544  PackingInfo(
545  BoundingBox(bBox2.getMinPt(), bBox2.getMaxPt() + overlap),
546  getPeriodicDimensions(),
547  getOrientation(),
548  getFaultMinRadius(),
549  getFaultMaxRadius()
550  )
551  );
552  }
553  else
554  {
555  std::stringstream msg;
556  msg
557  << "Roughness size plus gouge size is greater than block size: "
558  << "2*" << m_faultPrms.getSize() << " + " << m_gougePrms.getSize()
559  << " > " << getOrientationSize();
560  throw std::runtime_error(msg.str().c_str());
561  }
562  }
563  return infoVec;
564  }
565 
566  template <typename TPGrainGen>
568  {
569  return
570  std::max(
571  m_padRadius,
572  std::max(m_faultPrms.getMaxParticleRadius(), m_gougePrms.getMaxParticleRadius())
573  );
574  }
575 
576  template <typename TPGrainGen>
578  {
579  return
580  std::min(
581  m_padRadius,
582  std::min(m_faultPrms.getMinParticleRadius(), m_gougePrms.getMinParticleRadius())
583  );
584  }
585 
586  template <typename TPGrainGen>
588  {
589  return (getBBox().getSizes()[2] == 0.0);
590  }
591 
592  //==========================================================================
593 
594  template <typename TGPckr,typename TPPckr,typename TConn>
596  {
597  int numParticles = 0;
598  for (
599  typename GeneratorPtrVector::const_iterator it = m_genPtrVector.begin();
600  it != m_genPtrVector.end();
601  it++
602  )
603  {
604  numParticles += (*it)->getNumParticles();
605  }
606  return numParticles;
607  }
608 
609  template <typename TGPckr,typename TPPckr,typename TConn>
611  {
612  int numGrains = 0;
613  for (
614  typename GrainRndPackerPtrVector::const_iterator packerIt =
615  getGougeGeneratorVector().begin();
616  packerIt != getGougeGeneratorVector().end();
617  packerIt++
618  )
619  {
620  numGrains += (*packerIt)->getNumGrains();
621  }
622  return numGrains;
623  }
624 
625  template <typename TGPckr,typename TPPckr,typename TConn>
627  {
628  return m_connectionSet.size();
629  }
630 
631  template <typename TGPckr,typename TPPckr,typename TConn>
633  : m_nTablePtr(),
634  m_prms(prms),
635  m_connectionSet(),
636  m_gougeGenPtrVector(),
637  m_genPtrVector(),
638  m_particlePoolPtr(new ParticlePool(8*4096)),
639  m_grainPoolPtr(new GrainPool(4096)),
640  m_regularGenPtrVector(),
641  m_faultGenPtrVector()
642  {
643  /*
644  * Adjust the size of the ntable bounding-box to accommodate circular
645  * boundary conditions.
646  */
647  const BoundingBox bBox = m_prms.getBBox();
648  Vec3 ntableAdjust =
649  Vec3(
650  m_prms.getPeriodicDimensions()[0] ? 1 : 0,
651  m_prms.getPeriodicDimensions()[1] ? 1 : 0,
652  m_prms.getPeriodicDimensions()[2] ? 1 : 0
653  )*m_prms.getMaxRadius();
654  if (m_prms.getBBox().getSizes().Z() >= 4*m_prms.getMaxRadius()) {
655  ntableAdjust += Vec3(m_prms.getMaxRadius(), 0, 0);
656  }
657  const BoundingBox nTableBBox(bBox.getMinPt(), bBox.getMaxPt() - ntableAdjust);
658  m_nTablePtr =
659  NTablePtr(
660  new NTable(
661  nTableBBox,
662  (4.0*m_prms.getMinRadius()), // grid spacing
664  2.1*m_prms.getMaxRadius() // width of border-region in which
665  // particles are duplicated
666  // for circular boundry
667  )
668  );
669  }
670 
671  template <typename TGPckr,typename TPPckr,typename TConn>
673  {
674  BoundingBoxVector bBoxVector = m_prms.getRegularBBoxVector();
675  for (
676  BoundingBoxVector::const_iterator it = bBoxVector.begin();
677  it != bBoxVector.end();
678  it++
679  ) {
680  console.Debug()
681  << "GougeConfig<TGPckr,TPPckr,TConn>::createRegularBlockGenerators:"
682  << "Creating RegBoxPacker in box: " << StringUtil::toString(*it)
683  << "\n";
684 
685  GeneratorPtr genPtr =
686  GeneratorPtr(
687  new RegBoxPacker(
688  RegRadiusGenPtr(new RegRadiusGen(m_prms.getRegularBlockRadius())),
689  m_particlePoolPtr,
690  m_nTablePtr,
691  *it,
692  m_prms.getPeriodicDimensions(),
693  m_prms.getTolerance(),
694  m_prms.getRegularBlockRadius()
695  )
696  );
697  m_genPtrVector.push_back(genPtr);
698  m_regularGenPtrVector.push_back(genPtr);
699  }
700  }
701 
702  template <typename TGPckr,typename TPPckr,typename TConn>
704  {
705  PackingInfoVector infoVec = m_prms.getFaultPackingInfoVector();
706  for (
707  PackingInfoVector::const_iterator it = infoVec.begin();
708  it != infoVec.end();
709  it++
710  ) {
711  console.Debug()
712  << "GougeConfig<TGPckr,TPPckr,TConn>::createFaultBlockGenerators:"
713  << "Creating RndBoxPacker in box: " << StringUtil::toString(it->getBBox())
714  << "\n";
715  GeneratorPtr genPtr =
716  GeneratorPtr(
717  new RndBoxPacker(
719  new RndRadiusGen(
720  it->getMinParticleRadius(),
721  it->getMaxParticleRadius()
722  )
723  ),
724  m_particlePoolPtr,
725  m_nTablePtr,
726  it->getBBox(),
727  it->getPeriodicDimensions(),
728  m_prms.getTolerance(),
729  it->getMaxParticleRadius(),
730  m_prms.getMaxInsertionFailures(),
731  it->getFitPlaneVector()
732  )
733  );
734  m_genPtrVector.push_back(genPtr);
735  m_faultGenPtrVector.push_back(genPtr);
736  }
737  }
738 
739  template <typename TGPckr,typename TPPckr,typename TConn>
741  {
742  GougePackingInfoVector infoVec = m_prms.getGougePackingInfoVector();
743  for (
744  typename GougePackingInfoVector::const_iterator it = infoVec.begin();
745  it != infoVec.end();
746  it++
747  ) {
748  console.Debug()
749  << "GougeConfig<TGPckr,TPPckr,TConn>::createGougeConfigGenerators:"
750  << "Creating GrainRandomPacker in box: " << StringUtil::toString(it->getBBox())
751  << "\n";
752  GrainRandomPackerPtr genPtr =
754  new GrainRandomPacker(
755  typename GrainRandomPacker::ParticleGrainGenPtr(),
756  m_particlePoolPtr,
757  m_nTablePtr,
758  it->getBBox(),
759  it->getPeriodicDimensions(),
760  m_prms.getTolerance(),
761  it->getParticleGrainGen().getMaxGrainRadius(),
762  m_prms.getMaxInsertionFailures(),
763  it->getFitPlaneVector(),
764  m_grainPoolPtr
765  )
766  );
767  genPtr->setParticleGrainGen(it->getParticleGrainGen());
768  m_genPtrVector.push_back(genPtr);
769  m_gougeGenPtrVector.push_back(genPtr);
770  }
771  }
772 
773  template <typename TGPckr,typename TPPckr,typename TConn>
775  {
776  }
777 
778  template <typename TGPckr,typename TPPckr,typename TConn>
780  {
781  // setup block generators
782  createRegularBlockGenerators();
783  createFaultBlockGenerators();
784  createGougeConfigGenerators();
785 
786  // use block generators
787  console.Info() << "bbox = " << m_prms.getBBox().getMinPt() << " " << m_prms.getBBox().getMaxPt() << "\n";
788  for (
789  typename GeneratorPtrVector::iterator it = m_genPtrVector.begin();
790  it != m_genPtrVector.end();
791  it++
792  )
793  {
794  (*it)->generate();
795  }
796  createConnectionSet();
797  }
798 
799  template <typename TGPckr,typename TPPckr,typename TConn>
800  void GougeConfig<TGPckr,TPPckr,TConn>::writeToFile(const std::string &fileName) const
801  {
802  std::ofstream fStream(fileName.c_str());
803  write(fStream);
804  }
805 
806  template <typename TGPckr,typename TPPckr,typename TConn>
809  {
810  return m_gougeGenPtrVector;
811  }
812 
813  template <typename TGPckr,typename TPPckr,typename TConn>
816  {
817  return m_gougeGenPtrVector;
818  }
819 
820  template <typename TGPckr,typename TPPckr,typename TConn>
823  {
824  return m_faultGenPtrVector;
825  }
826 
827  template <typename TGPckr,typename TPPckr,typename TConn>
829  const Particle &p1,
830  const Particle &p2
831  ) const
832  {
833  const GeneratorPtrVector &generators = getFaultGeneratorVector();
834  if (generators.size() == 2) {
835  return
836  (
837  (generators[0]->contains(p1) && generators[1]->contains(p2))
838  ||
839  (generators[0]->contains(p2) && generators[1]->contains(p1))
840  );
841  }
842  else if (generators.size() > 2) {
843  throw
844  std::runtime_error(
845  "GougeConfig<TGPckr,TPPckr,TConn>::areInDifferentFaultBlocks: "
846  "More than two fault blocks."
847  );
848  }
849  return false;
850  }
851 
852  template <typename TGPckr,typename TPPckr,typename TConn>
854  {
855  const GrainRndPackerPtrVector &generators = getGougeGeneratorVector();
856  for (
857  typename GrainRndPackerPtrVector::const_iterator it = generators.begin();
858  it != generators.end();
859  it++
860  )
861  {
862  if ((*it)->contains(particle)) {
863  return true;
864  }
865  }
866  return false;
867  }
868 
869  template <typename TGPckr,typename TPPckr,typename TConn>
871  {
872  //
873  // First created connections in the elastic blocks.
874  //
875  typename NTable::ParticleIterator particleIt = m_nTablePtr->getParticleIterator();
876  ConnectionValidator validator = ConnectionValidator(*this, m_prms.getConnectionTolerance());
877  while (particleIt.hasNext()) {
878  const typename NTable::Particle *pParticle = particleIt.next();
879  const typename NTable::ParticleVector neighbours =
880  m_nTablePtr->getNeighbourVector(
881  pParticle->getPos(),
882  pParticle->getRad() + m_prms.getConnectionTolerance()
883  );
884  for (
885  typename NTable::ParticleVector::const_iterator it = neighbours.begin();
886  it != neighbours.end();
887  it++
888  )
889  {
890  if (validator.isValid(*pParticle, *(*it))) {
891  m_connectionSet.insert(
892  typename ConnectionSet::value_type(
893  pParticle->getID(),
894  (*it)->getID(),
895  m_prms.getBlockConnectionTag()
896  )
897  );
898  }
899  }
900  }
901  const int numBlockConns = m_connectionSet.size();
902  console.Info()
903  << "Created " << numBlockConns << " connections in "
904  << "bonded blocks.\n";
905 
906  //
907  // Create connections with grains.
908  //
909  console.Debug()
910  << "Prms BBox: " << StringUtil::toString(m_prms.getBBox()) << "\n";
911  console.Debug()
912  << "NTbl BBox: " << StringUtil::toString(m_nTablePtr->getBBox()) << "\n";
913  for (
914  typename GrainRndPackerPtrVector::iterator packerIt = getGougeGeneratorVector().begin();
915  packerIt != getGougeGeneratorVector().end();
916  packerIt++
917  )
918  {
919  typename GrainRandomPacker::GrainIterator grainIt = (*packerIt)->getGrainIterator();
920  while (grainIt.hasNext())
921  {
922  ConnectionFinder connFinder =
924  m_prms.getConnectionTolerance(),
925  m_prms.getGougeConnectionTag(),
926  m_nTablePtr->getBBox(),
927  m_nTablePtr->getPeriodicDimensions()
928  );
929  Grain &g = grainIt.next();
930  connFinder.create(g.getParticleIterator());
931  typename ConnectionFinder::Iterator connIt = connFinder.getIterator();
932  while (connIt.hasNext())
933  {
934  m_connectionSet.insert(connIt.next());
935  }
936  if (connFinder.getNumConnections() == 0)
937  {
938  console.Info()
939  << "Found no connections in grain " << g.getId()
940  << ":\n";
941  typename Grain::ParticleIterator partIt = g.getParticleIterator();
942  while (partIt.hasNext())
943  {
944  console.Info() << StringUtil::toString(partIt.next()) << "\n";
945  }
946  }
947  }
948  }
949  console.Info()
950  << "Created " << m_connectionSet.size()-numBlockConns << " connections in "
951  << "gouge region.\n";
952  }
953 
954  template <typename TGPckr,typename TPPckr,typename TConn>
957  {
958  ParticleCollection pCollection(m_particlePoolPtr);
959  for (
960  typename GeneratorPtrVector::iterator it = m_genPtrVector.begin();
961  it != m_genPtrVector.end();
962  it++
963  )
964  {
965  ParticleIterator particleIt = (*it)->getParticleIterator();
966  while (particleIt.hasNext()) {
967  pCollection.insertRef(particleIt.next());
968  }
969  }
970 
971  return pCollection;
972  }
973 
974  template <typename TGPckr,typename TPPckr,typename TConn>
977  {
978  GrainCollection gCollection(m_particlePoolPtr, m_grainPoolPtr);
979  for (
980  typename GrainRndPackerPtrVector::iterator packerIt =
981  getGougeGeneratorVector().begin();
982  packerIt != getGougeGeneratorVector().end();
983  packerIt++
984  )
985  {
986  GrainIterator grainIt = (*packerIt)->getGrainIterator();
987  while (grainIt.hasNext()) {
988  gCollection.insertRef(grainIt.next());
989  }
990  }
991 
992  return gCollection;
993  }
994 
995  template <typename TGPckr,typename TPPckr,typename TConn>
998  {
999  return m_connectionSet;
1000  }
1001 
1002  template <typename TGPckr,typename TPPckr,typename TConn>
1003  void GougeConfig<TGPckr,TPPckr,TConn>::write(std::ostream &oStream) const
1004  {
1005  Vec3 minPt = m_nTablePtr->getBBox().getMinPt();
1006  Vec3 maxPt = m_nTablePtr->getBBox().getMaxPt();
1007  if (fabs(maxPt.Z() - minPt.Z()) < (2*m_prms.getMaxRadius())) {
1008  minPt.Z() = minPt.Z() - m_prms.getMaxRadius() - m_prms.getTolerance();
1009  maxPt.Z() = maxPt.Z() + m_prms.getMaxRadius() + m_prms.getTolerance();
1010  }
1011 
1012  const BoundingBox geoBBox(minPt, maxPt + m_prms.getTolerance());
1013 
1014  GeometryInfo info =
1015  GeometryInfo(
1016  1.2,
1017  geoBBox.getMinPt(),
1018  geoBBox.getMaxPt(),
1019  m_prms.getPeriodicDimensions(),
1020  (m_prms.getBBox().getSizes().Z() <= 0.0)
1021  );
1022  info.write(oStream);
1023 
1024  /*
1025  * Some ugliness to ensure that duplicated particles (because of
1026  * circular boundary conditions) get the correct tag. First, create
1027  * a set of particles which have already been tagged.
1028  */
1029  typename NTable::ParticleIterator particleIt = m_nTablePtr->getParticleIterator();
1030  typedef std::set<Particle *, IdCompare> ParticleSet;
1031  ParticleSet taggedParticleSet;
1032  while (particleIt.hasNext()) {
1033  taggedParticleSet.insert(particleIt.next());
1034  }
1035 
1036  /*
1037  * Now eliminate any particles which were generated with their centre-point
1038  * lying outside the geo bounding box. Also set the tag in case the particle
1039  * was a duplicate, created due to circular boundary.
1040  */
1041  particleIt = m_nTablePtr->getParticleIterator();
1042  ParticleSet particleSet;
1043  while (particleIt.hasNext()) {
1044  Particle *pParticle = particleIt.next();
1045  if (geoBBox.contains(pParticle->getPos())) {
1046  pParticle->setTag((*(taggedParticleSet.find(pParticle)))->getTag());
1047  particleSet.insert(pParticle);
1048  }
1049  }
1050 
1051  /*
1052  * Write particles to the stream.
1053  */
1054  oStream
1055  << "\n"
1056  << "BeginParticles"
1057  << "\n"
1058  << "Simple"
1059  << "\n"
1060  << particleSet.size()
1061  << "\n";
1062  const int precision = 12;
1063  GeoParticleWriter particleVisitor(oStream, precision);
1064  for (
1065  typename ParticleSet::const_iterator it = particleSet.begin();
1066  it != particleSet.end();
1067  it++
1068  )
1069  {
1070  particleVisitor.visitParticle(*(*it));
1071  }
1072 
1073  oStream << "EndParticles\n" << "BeginConnect\n";
1074  oStream << getConnectionSet().size() << "\n";
1075  oStream.flush();
1076 
1077  GeoConnectionWriter connectionVisitor(oStream);
1078  visitConnections(connectionVisitor);
1079  oStream << "EndConnect";
1080  oStream.flush();
1081  }
1082 
1083  template <typename TGPckr,typename TPPckr,typename TConn>
1085  {
1086  for (
1087  typename GrainRndPackerPtrVector::iterator it = m_gougeGenPtrVector.begin();
1088  it != m_gougeGenPtrVector.end();
1089  it++
1090  )
1091  {
1092  typename GrainRandomPacker::ParticleIterator particleIt =
1093  (*it)->getParticleIterator();
1094  while (particleIt.hasNext()) {
1095  particleIt.next().setTag(tag);
1096  }
1097  }
1098  }
1099 
1100  template <typename TGPckr,typename TPPckr,typename TConn>
1102  {
1103  for (
1104  typename GeneratorPtrVector::iterator it = m_faultGenPtrVector.begin();
1105  it != m_faultGenPtrVector.end();
1106  it++
1107  )
1108  {
1109  ParticleIterator particleIt = (*it)->getParticleIterator();
1110  while (particleIt.hasNext()) {
1111  particleIt.next().setTag(tag);
1112  }
1113  }
1114  }
1115 
1116  template <typename TGPckr,typename TPPckr,typename TConn>
1118  int lowDrivingTag,
1119  int highDrivingTag,
1120  double distanceFromBBoxEdge
1121  )
1122  {
1123  ParticleCollection particleCollection = getParticleCollection();
1124  const BoundingBox bBox = particleCollection.getParticleBBox();
1125 
1126  const int idx = this->m_prms.getOrientationIndex();
1127  const double maxLow = bBox.getMinPt()[idx] + distanceFromBBoxEdge;
1128  const double minHigh = bBox.getMaxPt()[idx] - distanceFromBBoxEdge;
1129 
1130  int lowTagCount = 0;
1131  int highTagCount = 0;
1132  typename ParticleCollection::ParticleIterator particleIt =
1133  particleCollection.getParticleIterator();
1134  while (particleIt.hasNext())
1135  {
1136  Particle &particle = particleIt.next();
1137  const double dimPos = particle.getPos()[idx];
1138  const double radius = particle.getRad();
1139 
1140  if (dimPos - radius <= maxLow) {
1141  particle.setTag(lowDrivingTag);
1142  lowTagCount++;
1143  }
1144  if (dimPos + radius >= minHigh) {
1145  particle.setTag(highDrivingTag);
1146  highTagCount++;
1147  }
1148  }
1149  console.Info() << "Tagged " << lowTagCount << " particles with " << lowDrivingTag << "\n";
1150  console.Info() << "Tagged " << highTagCount << " particles with " << highDrivingTag << "\n";
1151  }
1152  }
1153 }