← Back to branch summary

~ubuntu-branches/ubuntu/quantal/python-demgengeo/quantal

« back to all changes in this revision

Viewing changes to src/HGrainGenerator.cc

  • Committer: Package Import Robot
  • Author(s): Anton Gladky
  • Date: 2011-11-18 21:47:18 UTC
  • Revision ID: package-import@ubuntu.com-20111118214718-4ysqm3dhpqwdd7gd
Tags: upstream-0.99~bzr106
Import upstream version 0.99~bzr106

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/HGrainGenerator.cc
 
1
/////////////////////////////////////////////////////////////
 
2
//                                                         //
 
3
// Copyright (c) 2007-2011 by The University of Queensland //
 
4
// Earth Systems Science Computational Centre (ESSCC)      //
 
5
// http://www.uq.edu.au/esscc                              //
 
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
//                                                         //
 
11
/////////////////////////////////////////////////////////////
 
12
 
 
13
// --- Project includes ---
 
14
#include "HGrainGenerator.h"
 
15
 
 
16
// --- system includes ---
 
17
#include <cmath>
 
18
 
 
19
using std::ceil;
 
20
using std::sqrt;
 
21
 
 
22
// --- IO includes ---
 
23
#include <iostream>
 
24
 
 
25
HGrainGenerator2D::HGrainGenerator2D()
 
26
{}
 
27
 
 
28
/*!
 
29
  constructor
 
30
 
 
31
  \param rad particle radius
 
32
*/
 
33
HGrainGenerator2D::HGrainGenerator2D(double rad)
 
34
{
 
35
  m_rad=rad;
 
36
}
 
37
 
 
38
/*!
 
39
  destructor - do nothing
 
40
*/
 
41
HGrainGenerator2D::~HGrainGenerator2D()
 
42
{}
 
43
 
 
44
/*!
 
45
  generate packing
 
46
  
 
47
  \param vol a pointer to the volume in which the packing is generated
 
48
  \param ntable a pointer to the neighbour table used 
 
49
  \param gid particle group id
 
50
  \param tag the particle tag
 
51
*/
 
52
void HGrainGenerator2D::generatePacking(AVolume2D* vol,MNTable2D* ntable,int gid, int tag) 
 
53
{
 
54
  const double tol=1e-5;
 
55
  // --- GENERATE REGULAR LATTICE ---
 
56
  // get bounding box
 
57
  pair<Vector3,Vector3> bbx=vol->getBoundingBox();
 
58
  double dx=(bbx.second.X()-bbx.first.X())-2.0*m_rad;
 
59
  double dy=(bbx.second.Y()-bbx.first.Y())-2.0*m_rad;
 
60
  // get index limits for seeding
 
61
  int imax=int(floor(dx/(m_rad*2.0)))+1;
 
62
  if(dx-(double(imax)*m_rad*2.0)>0.5*m_rad) imax++;
 
63
  int jmax=int(floor(dy/(m_rad*sqrt(3.0))))+1;
 
64
  // check if odd or even
 
65
  bool even=(dx-(double(imax)*m_rad*2.0)>0.5*m_rad);
 
66
  std::cerr << "imax, jmax, even  " << imax << " " << jmax << " " << even << std::endl; 
 
67
  // seed positions
 
68
  for(int i=0;i<imax-1;i++){
 
69
    for(int j=0;j<jmax;j++){
 
70
      // get position
 
71
      double px=bbx.first.X()+tol+m_rad+(double(i)+0.5*double(j%2))*m_rad*2.0;
 
72
      double py=bbx.first.Y()+tol+m_rad+double(j)*sqrt(3.0)*m_rad;
 
73
 
 
74
      Sphere S(Vector3(px,py,0.0),m_rad);
 
75
      S.setTag(tag);
 
76
      ntable->insert(S,gid);    
 
77
    }
 
78
  }
 
79
  for(int j=0;j<jmax;j++){
 
80
    if(even || ((j%2)==0)){
 
81
      // get position
 
82
      double px=bbx.first.X()+tol+m_rad+(double(imax-1)+0.5*double(j%2))*m_rad*2.0;
 
83
      double py=bbx.first.Y()+tol+m_rad+double(j)*sqrt(3.0)*m_rad;
 
84
 
 
85
      Sphere S(Vector3(px,py,0.0),m_rad);
 
86
      S.setTag(tag);
 
87
      ntable->insert(S,gid);    
 
88
    }
 
89
  }
 
90
  // pick hex grains
 
91
  if(!even){
 
92
    double x0=bbx.first.X()+tol+4.0*m_rad;
 
93
    double y0=bbx.first.Y()+tol+(1.0+sqrt(3.0))*m_rad;
 
94
    double dxi=5.0*m_rad;
 
95
    double dxi3=m_rad;
 
96
    double dyi=sqrt(3.0)*m_rad;
 
97
    int igmax=ceil(dx/5.0*m_rad);
 
98
    double dxj=m_rad;
 
99
    double dyj=3.0*sqrt(3.0)*m_rad;
 
100
    double dyj3=sqrt(3.0)*m_rad;
 
101
    int jgmax=ceil(dy/(3.0*sqrt(3.0))*m_rad);
 
102
    for(int i=0;i<igmax;i++){
 
103
      for(int j=0;j<jgmax;j++){
 
104
        double px=x0+double(i)*dxi-double(i/3)*dxi3+double(j%5)*dxj;
 
105
        double py=y0+double(i%3)*dyi+double(j)*dyj-double(j/5)*dyj3;
 
106
        // check if far enough inside
 
107
        if((px-bbx.first.X()>=3.0*m_rad) && 
 
108
           (bbx.second.X()-px>=3.0*m_rad)&& 
 
109
           (py-bbx.first.Y()>=(sqrt(3.0)+1.0)*m_rad) &&
 
110
           (bbx.second.Y()-py>=(sqrt(3.0)+1.0)*m_rad)){
 
111
          ntable->tagParticlesNear(Vector3(px,py,0.0),m_rad+tol,gid,2);
 
112
          ntable->generateBondsWithMask(gid,tol,2,2,2);
 
113
          ntable->tagParticlesNear(Vector3(px,py,0.0),m_rad+tol,gid,1);
 
114
        }
 
115
      }
 
116
    }
 
117
  }
 
118
  ntable->removeTagged(gid,0,7);
 
119
}
 
120
 
 
121
ostream& operator << (ostream& ost,const HGrainGenerator2D& T)
 
122
{
 
123
   return ost;
 
124
}
 
125