branches/exact/examples/no_release/sphere_mesh.cc

   1 // Voronoi calculation example code
   2 //
   3 // Author   : Chris H. Rycroft (LBL / UC Berkeley)
   4 // Email    : chr@alum.mit.edu
   5 // Date     : August 30th 2011
   6
   7 #include "voro++.hh"
   8 using namespace voro;
   9
  10 // Set up constants for the container geometry
  11 const double boxl=1.2;
  12
  13 // Set up the number of blocks that the container is divided into
  14 const int bl=24;
  15
  16 // Set the number of particles that are going to be randomly introduced
  17 const int particles=20000;
  18
  19 const int nface=11;
  20
  21 // This function returns a random double between 0 and 1
  22 double rnd() {return double(rand())/RAND_MAX;}
  23
  24 struct wall_shell : public wall {
  25         public:
  26                 wall_shell(double xc_,double yc_,double zc_,double rc,double sc,int w_id_=-99)
  27                         : w_id(w_id_), xc(xc_), yc(yc_), zc(zc_), lc(rc-sc), uc(rc+sc) {}
  28                 bool point_inside(double x,double y,double z) {
  29                         double rsq=(x-xc)*(x-xc)+(y-yc)*(y-yc)+(z-zc)*(z-zc);
  30                         return rsq>lc*lc&&rsq<uc*uc;
  31                 }
  32                 template<class v_cell>
  33                 bool cut_cell_base(v_cell &c,double x,double y,double z) {
  34                         double xd=x-xc,yd=y-yc,zd=z-zc,dq=xd*xd+yd*yd+zd*zd,dq2;
  35                         if (dq>1e-5) {
  36                                 dq2=2*(sqrt(dq)*lc-dq);
  37                                 dq=2*(sqrt(dq)*uc-dq);
  38                                 return c.nplane(xd,yd,zd,dq,w_id)&&c.nplane(-xd,-yd,-zd,-dq2,w_id);
  39                         }
  40                         return true;
  41                 }
  42                 bool cut_cell(voronoicell &c,double x,double y,double z) {return cut_cell_base(c,x,y,z);}
  43                 bool cut_cell(voronoicell_neighbor &c,double x,double y,double z) {return cut_cell_base(c,x,y,z);}
  44         private:
  45                 const int w_id;
  46                 const double xc,yc,zc,lc,uc;
  47 };
  48
  49
  50 int main() {
  51         int i=0,j,k,l,ll,o;
  52         double x,y,z,r,dx,dy,dz;
  53         int faces[nface],*fp;
  54         double p[3*particles];
  55
  56         // Create a container with the geometry given above, and make it
  57         // non-periodic in each of the three coordinates. Allocate space for
  58         // eight particles within each computational block
  59         container con(-boxl,boxl,-boxl,boxl,-boxl,boxl,bl,bl,bl,false,false,false,8);
  60
  61         wall_shell ws(0,0,0,1,0.00001);
  62         con.add_wall(ws);
  63
  64         // Randomly add particles into the container
  65         while(i<particles) {
  66                 x=boxl*(2*rnd()-1);
  67                 y=boxl*(2*rnd()-1);
  68                 z=boxl*(2*rnd()-1);
  69                 r=x*x+y*y+z*z;
  70                 if(r>1e-5) {
  71                         r=1/sqrt(r);x*=r;y*=r;z*=r;
  72                         con.put(i,x,y,z);
  73                         i++;
  74                 }
  75         }
  76
  77         for(l=0;l<1000;l++) {
  78                 c_loop_all vl(con);
  79                 voronoicell c;
  80                 for(fp=faces;fp<faces+nface;fp++) *fp=0;
  81                 if(vl.start()) do if(con.compute_cell(c,vl)) {
  82                         vl.pos(i,x,y,z,r);
  83                         c.centroid(dx,dy,dz);
  84                         p[3*i]=x+dx;
  85                         p[3*i+1]=y+dy;
  86                         p[3*i+2]=z+dz;
  87
  88                         i=c.number_of_faces()-4;
  89                         if(i<0) i=0;if(i>=nface) i=nface-1;
  90                         faces[i]++;
  91                 } while (vl.inc());
  92                 con.clear();
  93                 double fac=0;//l<9000?0.1/sqrt(double(l)):0;
  94                 for(i=0;i<particles;i++) con.put(i,p[3*i]+fac*(2*rnd()-1),p[3*i+1]+fac*(2*rnd()-1),p[3*i+2]+fac*(2*rnd()-1));
  95                 printf("%d",l);
  96                 for(fp=faces;fp<faces+nface;fp++) printf(" %d",*fp);
  97                 puts("");
  98         }
  99
 100         // Output the particle positions in gnuplot format
 101         con.draw_particles("sphere_mesh_p.gnu");
 102
 103         // Output the Voronoi cells in gnuplot format
 104         con.draw_cells_gnuplot("sphere_mesh_v.gnu");
 105
 106         // Allocate memory for neighbor relations
 107         int *q=new int[particles*nface],*qn=new int[particles],*qp;
 108         for(l=0;l<particles;l++) qn[l]=0;
 109
 110         // Create a table of all neighbor relations
 111         vector<int> vi;
 112         voronoicell_neighbor c;
 113         c_loop_all vl(con);
 114         if(vl.start()) do if(con.compute_cell(c,vl)) {
 115                 i=vl.pid();qp=q+i*nface;
 116                 c.neighbors(vi);
 117                 if(vi.size()>nface+2) voro_fatal_error("Too many faces; boost nface",5);
 118
 119                 for(l=0;l<(signed int) vi.size();l++) if(vi[l]>=0) qp[qn[i]++]=vi[l];
 120         } while (vl.inc());
 121
 122         // Sort the connections in anti-clockwise order
 123         bool connect;
 124         int tote=0;
 125         for(l=0;l<particles;l++) {
 126                 tote+=qn[l];
 127                 for(i=0;i<qn[l]-2;i++) {
 128                         o=q[l*nface+i];
 129                         //printf("---> %d,%d\n",i,o);
 130                         j=i+1;
 131                         while(j<qn[l]-1) {
 132                                 ll=q[l*nface+j];
 133                         //      printf("-> %d %d\n",j,ll);
 134                                 connect=false;
 135                                 for(k=0;k<qn[ll];k++) {
 136                         //              printf("%d %d %d\n",ll,k,q[ll*nface+k]);
 137                                         if(q[ll*nface+k]==o) {connect=true;break;}
 138                                 }
 139                                 if(connect) break;
 140                                 j++;
 141                         }
 142
 143                         // Swap the connected vertex into this location
 144                         //printf("%d %d\n",i+1,j);
 145                         o=q[l*nface+i+1];
 146                         q[l*nface+i+1]=q[l*nface+j];
 147                         q[l*nface+j]=o;
 148                 }
 149
 150                 // Reverse all connections if the have the wrong handedness
 151                 j=3*l;k=3*q[l*nface];o=3*q[l*nface+1];
 152                 x=p[j]-p[k];dx=p[j]-p[o];
 153                 y=p[j+1]-p[k+1];dy=p[j+1]-p[o+1];
 154                 z=p[j+2]-p[k+2];dz=p[j+2]-p[o+2];
 155                 if(p[j]*(y*dz-z*dy)+p[j+1]*(z*dx-x*dz)+p[j+2]*(x*dy-y*dx)<0) {
 156                         for(i=0;i<qn[l]/2;i++) {
 157                                 o=q[l*nface+i];
 158                                 q[l*nface+i]=q[l*nface+qn[l]-1-i];
 159                                 q[l*nface+qn[l]-1-i]=o;
 160                         }
 161                 }
 162         }
 163
 164         FILE *ff=safe_fopen("sphere_mesh.net","w");
 165         int *mp=new int[particles],*mpi=new int[particles];
 166         for(i=0;i<particles;i++) mp[i]=-1;
 167         *mpi=0;*mp=0;l=1;o=0;
 168         while(o<l) {
 169                 i=mpi[o];
 170                 for(j=0;j<qn[i];j++) {
 171                         k=q[i*nface+j];
 172                         if(mp[k]==-1) {
 173                                 mpi[l]=k;
 174                                 mp[k]=l++;
 175                         }
 176                         if(mp[i]<mp[k])
 177                                 fprintf(ff,"%g %g %g\n%g %g %g\n\n\n",p[3*i],p[3*i+1],p[3*i+2],p[3*k],p[3*k+1],p[3*k+2]);
 178                 }
 179                 o++;
 180         }
 181         fclose(ff);
 182
 183         // Save binary representation of the mesh
 184         FILE *fb=safe_fopen("sphere_mesh.bin","wb");
 185
 186         // Write header
 187         int kk[3],sz=tote+particles+2,*red(new int[sz]),*rp=red;
 188         *kk=1;kk[1]=sz;kk[2]=3*particles;
 189         fwrite(kk,sizeof(int),3,fb);
 190
 191         // Assemble the connections and write them
 192         *(rp++)=particles;*(rp++)=tote;
 193         for(l=0;l<particles;l++) *(rp++)=qn[mpi[l]];
 194         for(l=0;l<particles;l++) {
 195                 i=mpi[l];printf("%d",l);
 196                 for(j=0;j<qn[i];j++) {*(rp++)=mp[q[i*nface+j]];printf(" %d",*(rp-1));}
 197                 puts("");
 198         }
 199         fwrite(red,sizeof(int),sz,fb);
 200
 201
 202         double *pm=new double[3*particles],*a=pm,*b;
 203         for(i=0;i<particles;i++) {
 204                 b=p+3*mpi[i];
 205                 *(a++)=*(b++);*(a++)=*(b++);*(a++)=*b;
 206         }
 207         fwrite(pm,sizeof(double),3*particles,fb);
 208         delete [] pm;
 209
 210         // Free dynamically allocated arrays
 211         delete [] red;
 212         delete [] mpi;
 213         delete [] mp;
 214         delete [] qn;
 215         delete [] q;
 216 }