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