trunk/examples/no_release/polycrystal_rahman.cc

   1 // Voronoi method to generate nanocrystalline grain boundaries
   2 // Nov 18, 2011
   3
   4 #include <cstdio>
   5 #include <cstdlib>
   6 #include <cmath>
   7 using namespace std;
   8
   9 #include "voro++.hh"
  10 using namespace voro;
  11
  12 const double pi=3.1415926535897932384626433832795;
  13
  14 // Box geometry
  15 const double x_min=0,x_max=81;
  16 const double y_min=0,y_max=81;
  17 const double z_min=0,z_max=81;
  18 const double cvol=(x_max-x_min)*(y_max-y_min)*(x_max-x_min);
  19
  20 // Total number of particles
  21 const int particles=24;
  22
  23 // Lattice size
  24 const double h=4;
  25
  26 // Function for random double between 0 and 1
  27 double rnd() {return double(rand())/RAND_MAX;}
  28
  29 int main() {
  30         int i,j=0,c_id;
  31         double x,y,z,x1,y1,z1,x2,y2,z2,x3,y3,z3,xt,yt,zt,rx,ry,rz;
  32         double theta,cth,sth,r,v,xc,yc,zc,vx,vy,vz;
  33
  34         // Create the container class
  35         container con(x_min,x_max,y_min,y_max,z_min,z_max,10,10,10,
  36                         false,false,false,8);
  37
  38         // Add generators to the box
  39         for(i=1;i<=particles;i++) {
  40                 x=x_min+rnd()*(x_max-x_min);
  41                 y=y_min+rnd()*(y_max-y_min);
  42                 z=z_min+rnd()*(z_max-z_min);
  43                 con.put(i,x,y,z);
  44         }
  45
  46         // Open the file for the particle positions
  47         FILE *fp=safe_fopen("lammps_input","w");
  48
  49         // Create a loop class to iterate over all of the generators
  50         // in the box
  51         c_loop_all cl(con);
  52         voronoicell c;
  53         if(cl.start()) do if(con.compute_cell(c,cl)) {
  54
  55                 // Generate the first vector of an orthonormal basis
  56                 x1=2*rnd()-1;y1=2*rnd()-1;z1=2*rnd()-1;r=1/sqrt(x1*x1+y1*y1+z1*z1);
  57                 x1*=r;y1*=r;z1*=r;
  58
  59                 // Construct a second perpendicular vector
  60                 if(abs(x1)>0.5||abs(y1)>0.5) {r=1/sqrt(x1*x1+y1*y1);x2=-y1*r;y2=x1*r;z2=0;}
  61                 else {r=1/sqrt(x1*x1+z1*z1);x2=-z1*r;y2=0;z2=x1*r;}
  62
  63                 // Construct a third perpendicular vector using the vector product
  64                 x=y2*z1-z2*y1;y=z2*x1-x2*z1;z=x2*y1-y2*x1;
  65
  66                 // Take a random rotation of the second and third vectors
  67                 theta=2*pi*rnd();cth=cos(theta);sth=sin(theta);
  68                 x3=x*cth+x2*sth;x2=-x*sth+x2*cth;
  69                 y3=y*cth+y2*sth;y2=-y*sth+y2*cth;
  70                 z3=z*cth+z2*sth;z2=-z*sth+z2*cth;
  71
  72                 // Get a bound on how far to search
  73                 r=sqrt(0.25*c.max_radius_squared());
  74                 v=(int(r/h)+2)*h;
  75
  76                 // Add small random displacement to lattice positioning,
  77                 // so that it's not always perfectly aligned with the generator
  78                 vx=-v+h*rnd();vy=-v+h*rnd();vz=-v+h*rnd();
  79
  80                 // Print diagnostic information about this generator
  81                 c_id=cl.pid();cl.pos(xc,yc,zc);
  82                 printf("Generator %d at (%g,%g,%g), random basis:\n",c_id,xc,yc,zc);
  83                 printf("%g %g %g\n",x1,y1,z1);
  84                 printf("%g %g %g\n",x2,y2,z2);
  85                 printf("%g %g %g\n\n",x3,y3,z3);
  86
  87                 // Loop over a local region of points
  88                 for(z=vx;z<=v;z+=h) for(y=vy;y<=v;y+=h) for(x=vz;x<=v;x+=h) {
  89
  90                         // Construct points rotated into the random basis
  91                         xt=xc+x*x1+y*x2+z*x3;
  92                         yt=yc+x*y1+y*y2+z*y3;
  93                         zt=zc+x*z1+y*z2+z*z3;
  94
  95                         // Skip if this lies outside the container
  96                         if(xt<x_min||xt>x_max||yt<y_min||yt>y_max||zt<z_min||zt>z_max) continue;
  97
  98                         // Find the nearest generator
  99                         con.find_voronoi_cell(xt,yt,zt,rx,ry,rz,i);
 100
 101                         // If the nearest generator matches, then save this point
 102                         if(i==c_id) {fprintf(fp,"%d %g %g %g\n",j,xt,yt,zt);j++;}
 103                 }
 104         } while(cl.inc());
 105
 106         // Close the output file
 107         fclose(fp);
 108
 109         // Output files for diagnosic purposes
 110         con.draw_particles("lammps_generators");
 111         con.draw_cells_gnuplot("lammps_cells");
 112 }