branches/exact/src/container_prd.cc

   1 // Voro++, a 3D cell-based Voronoi library
   2 //
   3 // Author   : Chris H. Rycroft (LBL / UC Berkeley)
   4 // Email    : chr@alum.mit.edu
   5 // Date     : August 30th 2011
   6
   7 /** \file container_prd.cc
   8  * \brief Function implementations for the container_periodic_base and
   9  * related classes. */
  10
  11 #include "container_prd.hh"
  12
  13 namespace voro {
  14
  15 /** The class constructor sets up the geometry of container, initializing the
  16  * minimum and maximum coordinates in each direction, and setting whether each
  17  * direction is periodic or not. It divides the container into a rectangular
  18  * grid of blocks, and allocates memory for each of these for storing particle
  19  * positions and IDs.
  20  * \param[in] (bx_) The x coordinate of the first unit vector.
  21  * \param[in] (bxy_,by_) The x and y coordinates of the second unit vector.
  22  * \param[in] (bxz_,byz_,bz_) The x, y, and z coordinates of the third unit
  23  *                            vector.
  24  * \param[in] (nx_,ny_,nz_) the number of grid blocks in each of the three
  25  *                       coordinate directions.
  26  * \param[in] init_mem_ the initial memory allocation for each block.
  27  * \param[in] ps_ the number of floating point entries to store for each
  28  *                particle. */
  29 container_periodic_base::container_periodic_base(double bx_,double bxy_,double by_,
  30                 double bxz_,double byz_,double bz_,int nx_,int ny_,int nz_,int init_mem_,int ps_)
  31         : unitcell(bx_,bxy_,by_,bxz_,byz_,bz_), voro_base(nx_,ny_,nz_,bx_/nx_,by_/ny_,bz_/nz_),
  32         ey(int(max_uv_y*ysp+1)), ez(int(max_uv_z*zsp+1)), wy(ny+ey), wz(nz+ez),
  33         oy(ny+2*ey), oz(nz+2*ez), oxyz(nx*oy*oz), id(new int*[oxyz]), p(new double*[oxyz]),
  34         co(new int[oxyz]), mem(new int[oxyz]), img(new char[oxyz]), init_mem(init_mem_), ps(ps_) {
  35         int i,j,k,l;
  36
  37         // Clear the global arrays
  38         int *pp=co;while(pp<co+oxyz) *(pp++)=0;
  39         pp=mem;while(pp<mem+oxyz) *(pp++)=0;
  40         char *cp=img;while(cp<img+oxyz) *(cp++)=0;
  41
  42         // Set up memory for the blocks in the primary domain
  43         for(k=ez;k<wz;k++) for(j=ey;j<wy;j++) for(i=0;i<nx;i++) {
  44                 l=i+nx*(j+oy*k);
  45                 mem[l]=init_mem;
  46                 id[l]=new int[init_mem];
  47                 p[l]=new double[ps*init_mem];
  48         }
  49 }
  50
  51 /** The container destructor frees the dynamically allocated memory. */
  52 container_periodic_base::~container_periodic_base() {
  53         int l;
  54         for(l=0;l<oxyz;l++) if(mem[l]>0) delete [] p[l];
  55         for(l=0;l<oxyz;l++) if(mem[l]>0) delete [] id[l];
  56         delete [] p;
  57         delete [] id;
  58         delete [] mem;
  59         delete [] co;
  60 }
  61
  62 /** The class constructor sets up the geometry of container.
  63  * \param[in] (bx_) The x coordinate of the first unit vector.
  64  * \param[in] (bxy_,by_) The x and y coordinates of the second unit vector.
  65  * \param[in] (bxz_,byz_,bz_) The x, y, and z coordinates of the third unit
  66  *                            vector.
  67  * \param[in] (nx_,ny_,nz_) the number of grid blocks in each of the three
  68  *                          coordinate directions.
  69  * \param[in] init_mem_ the initial memory allocation for each block. */
  70 container_periodic::container_periodic(double bx_,double bxy_,double by_,double bxz_,double byz_,double bz_,
  71         int nx_,int ny_,int nz_,int init_mem_)
  72         : container_periodic_base(bx_,bxy_,by_,bxz_,byz_,bz_,nx_,ny_,nz_,init_mem_,3),
  73         vc(*this,2*nx_+1,2*ey+1,2*ez+1) {}
  74
  75 /** The class constructor sets up the geometry of container.
  76  * \param[in] (bx_) The x coordinate of the first unit vector.
  77  * \param[in] (bxy_,by_) The x and y coordinates of the second unit vector.
  78  * \param[in] (bxz_,byz_,bz_) The x, y, and z coordinates of the third unit
  79  *                            vector.
  80  * \param[in] (nx_,ny_,nz_) the number of grid blocks in each of the three
  81  *                          coordinate directions.
  82  * \param[in] init_mem_ the initial memory allocation for each block. */
  83 container_periodic_poly::container_periodic_poly(double bx_,double bxy_,double by_,double bxz_,double byz_,double bz_,
  84         int nx_,int ny_,int nz_,int init_mem_)
  85         : container_periodic_base(bx_,bxy_,by_,bxz_,byz_,bz_,nx_,ny_,nz_,init_mem_,4),
  86         vc(*this,2*nx_+1,2*ey+1,2*ez+1) {ppr=p;}
  87
  88 /** Put a particle into the correct region of the container.
  89  * \param[in] n the numerical ID of the inserted particle.
  90  * \param[in] (x,y,z) the position vector of the inserted particle. */
  91 void container_periodic::put(int n,double x,double y,double z) {
  92         int ijk;
  93         put_locate_block(ijk,x,y,z);
  94         id[ijk][co[ijk]]=n;
  95         double *pp=p[ijk]+3*co[ijk]++;
  96         *(pp++)=x;*(pp++)=y;*pp=z;
  97 }
  98
  99 /** Put a particle into the correct region of the container.
 100  * \param[in] n the numerical ID of the inserted particle.
 101  * \param[in] (x,y,z) the position vector of the inserted particle.
 102  * \param[in] r the radius of the particle. */
 103 void container_periodic_poly::put(int n,double x,double y,double z,double r) {
 104         int ijk;
 105         put_locate_block(ijk,x,y,z);
 106         id[ijk][co[ijk]]=n;
 107         double *pp=p[ijk]+4*co[ijk]++;
 108         *(pp++)=x;*(pp++)=y;*(pp++)=z;*pp=r;
 109         if(max_radius<r) max_radius=r;
 110 }
 111
 112 /** Put a particle into the correct region of the container.
 113  * \param[in] n the numerical ID of the inserted particle.
 114  * \param[in] (x,y,z) the position vector of the inserted particle.
 115  * \param[out] (ai,aj,ak) the periodic image displacement that the particle is
 116  *                        in, with (0,0,0) corresponding to the primary domain.
 117  */
 118 void container_periodic::put(int n,double x,double y,double z,int &ai,int &aj,int &ak) {
 119         int ijk;
 120         put_locate_block(ijk,x,y,z,ai,aj,ak);
 121         id[ijk][co[ijk]]=n;
 122         double *pp=p[ijk]+3*co[ijk]++;
 123         *(pp++)=x;*(pp++)=y;*pp=z;
 124 }
 125
 126 /** Put a particle into the correct region of the container.
 127  * \param[in] n the numerical ID of the inserted particle.
 128  * \param[in] (x,y,z) the position vector of the inserted particle.
 129  * \param[in] r the radius of the particle.
 130  * \param[out] (ai,aj,ak) the periodic image displacement that the particle is
 131  *                        in, with (0,0,0) corresponding to the primary domain.
 132  */
 133 void container_periodic_poly::put(int n,double x,double y,double z,double r,int &ai,int &aj,int &ak) {
 134         int ijk;
 135         put_locate_block(ijk,x,y,z,ai,aj,ak);
 136         id[ijk][co[ijk]]=n;
 137         double *pp=p[ijk]+4*co[ijk]++;
 138         *(pp++)=x;*(pp++)=y;*(pp++)=z;*pp=r;
 139         if(max_radius<r) max_radius=r;
 140 }
 141
 142 /** Put a particle into the correct region of the container, also recording
 143  * into which region it was stored.
 144  * \param[in] vo the ordering class in which to record the region.
 145  * \param[in] n the numerical ID of the inserted particle.
 146  * \param[in] (x,y,z) the position vector of the inserted particle. */
 147 void container_periodic::put(particle_order &vo,int n,double x,double y,double z) {
 148         int ijk;
 149         put_locate_block(ijk,x,y,z);
 150         id[ijk][co[ijk]]=n;
 151         vo.add(ijk,co[ijk]);
 152         double *pp=p[ijk]+3*co[ijk]++;
 153         *(pp++)=x;*(pp++)=y;*pp=z;
 154 }
 155
 156 /** Put a particle into the correct region of the container, also recording
 157  * into which region it was stored.
 158  * \param[in] vo the ordering class in which to record the region.
 159  * \param[in] n the numerical ID of the inserted particle.
 160  * \param[in] (x,y,z) the position vector of the inserted particle.
 161  * \param[in] r the radius of the particle. */
 162 void container_periodic_poly::put(particle_order &vo,int n,double x,double y,double z,double r) {
 163         int ijk;
 164         put_locate_block(ijk,x,y,z);
 165         id[ijk][co[ijk]]=n;
 166         vo.add(ijk,co[ijk]);
 167         double *pp=p[ijk]+4*co[ijk]++;
 168         *(pp++)=x;*(pp++)=y;*(pp++)=z;*pp=r;
 169         if(max_radius<r) max_radius=r;
 170 }
 171
 172 /** Takes a particle position vector and computes the region index into which
 173  * it should be stored. If the container is periodic, then the routine also
 174  * maps the particle position to ensure it is in the primary domain. If the
 175  * container is not periodic, the routine bails out.
 176  * \param[out] ijk the region index.
 177  * \param[in,out] (x,y,z) the particle position, remapped into the primary
 178  *                        domain if necessary.
 179  * \return True if the particle can be successfully placed into the container,
 180  * false otherwise. */
 181 void container_periodic_base::put_locate_block(int &ijk,double &x,double &y,double &z) {
 182
 183         // Remap particle in the z direction if necessary
 184         int k=step_int(z*zsp);
 185         if(k<0||k>=nz) {
 186                 int ak=step_div(k,nz);
 187                 z-=ak*bz;y-=ak*byz;x-=ak*bxz;k-=ak*nz;
 188         }
 189
 190         // Remap particle in the y direction if necessary
 191         int j=step_int(y*ysp);
 192         if(j<0||j>=ny) {
 193                 int aj=step_div(j,ny);
 194                 y-=aj*by;x-=aj*bxy;j-=aj*ny;
 195         }
 196
 197         // Remap particle in the x direction if necessary
 198         ijk=step_int(x*xsp);
 199         if(ijk<0||ijk>=nx) {
 200                 int ai=step_div(ijk,nx);
 201                 x-=ai*bx;ijk-=ai*nx;
 202         }
 203
 204         // Compute the block index and check memory allocation
 205         j+=ey;k+=ez;
 206         ijk+=nx*(j+oy*k);
 207         if(co[ijk]==mem[ijk]) add_particle_memory(ijk);
 208 }
 209
 210 /** Takes a particle position vector and computes the region index into which
 211  * it should be stored. If the container is periodic, then the routine also
 212  * maps the particle position to ensure it is in the primary domain. If the
 213  * container is not periodic, the routine bails out.
 214  * \param[out] ijk the region index.
 215  * \param[in,out] (x,y,z) the particle position, remapped into the primary
 216  *                        domain if necessary.
 217  * \param[out] (ai,aj,ak) the periodic image displacement that the particle is
 218  *                        in, with (0,0,0) corresponding to the primary domain.
 219  * \return True if the particle can be successfully placed into the container,
 220  * false otherwise. */
 221 void container_periodic_base::put_locate_block(int &ijk,double &x,double &y,double &z,int &ai,int &aj,int &ak) {
 222
 223         // Remap particle in the z direction if necessary
 224         int k=step_int(z*zsp);
 225         if(k<0||k>=nz) {
 226                 ak=step_div(k,nz);
 227                 z-=ak*bz;y-=ak*byz;x-=ak*bxz;k-=ak*nz;
 228         } else ak=0;
 229
 230         // Remap particle in the y direction if necessary
 231         int j=step_int(y*ysp);
 232         if(j<0||j>=ny) {
 233                 aj=step_div(j,ny);
 234                 y-=aj*by;x-=aj*bxy;j-=aj*ny;
 235         } else aj=0;
 236
 237         // Remap particle in the x direction if necessary
 238         ijk=step_int(x*xsp);
 239         if(ijk<0||ijk>=nx) {
 240                 ai=step_div(ijk,nx);
 241                 x-=ai*bx;ijk-=ai*nx;
 242         } else ai=0;
 243
 244         // Compute the block index and check memory allocation
 245         j+=ey;k+=ez;
 246         ijk+=nx*(j+oy*k);
 247         if(co[ijk]==mem[ijk]) add_particle_memory(ijk);
 248 }
 249
 250 /** Takes a position vector and remaps it into the primary domain.
 251  * \param[out] (ai,aj,ak) the periodic image displacement that the vector is in,
 252  *                        with (0,0,0) corresponding to the primary domain.
 253  * \param[out] (ci,cj,ck) the index of the block that the position vector is
 254  *                        within, once it has been remapped.
 255  * \param[in,out] (x,y,z) the position vector to consider, which is remapped
 256  *                        into the primary domain during the routine.
 257  * \param[out] ijk the block index that the vector is within. */
 258 inline void container_periodic_base::remap(int &ai,int &aj,int &ak,int &ci,int &cj,int &ck,double &x,double &y,double &z,int &ijk) {
 259
 260         // Remap particle in the z direction if necessary
 261         ck=step_int(z*zsp);
 262         if(ck<0||ck>=nz) {
 263                 ak=step_div(ck,nz);
 264                 z-=ak*bz;y-=ak*byz;x-=ak*bxz;ck-=ak*nz;
 265         } else ak=0;
 266
 267         // Remap particle in the y direction if necessary
 268         cj=step_int(y*ysp);
 269         if(cj<0||cj>=ny) {
 270                 aj=step_div(cj,ny);
 271                 y-=aj*by;x-=aj*bxy;cj-=aj*ny;
 272         } else aj=0;
 273
 274         // Remap particle in the x direction if necessary
 275         ci=step_int(x*xsp);
 276         if(ci<0||ci>=nx) {
 277                 ai=step_div(ci,nx);
 278                 x-=ai*bx;ci-=ai*nx;
 279         } else ai=0;
 280
 281         cj+=ey;ck+=ez;
 282         ijk=ci+nx*(cj+oy*ck);
 283 }
 284
 285 /** Takes a vector and finds the particle whose Voronoi cell contains that
 286  * vector. This is equivalent to finding the particle which is nearest to the
 287  * vector.
 288  * \param[in] (x,y,z) the vector to test.
 289  * \param[out] (rx,ry,rz) the position of the particle whose Voronoi cell
 290  *                        contains the vector. This may point to a particle in
 291  *                        a periodic image of the primary domain.
 292  * \param[out] pid the ID of the particle.
 293  * \return True if a particle was found. If the container has no particles,
 294  * then the search will not find a Voronoi cell and false is returned. */
 295 bool container_periodic::find_voronoi_cell(double x,double y,double z,double &rx,double &ry,double &rz,int &pid) {
 296         int ai,aj,ak,ci,cj,ck,ijk;
 297         particle_record w;
 298         double mrs;
 299
 300         // Remap the vector into the primary domain and then search for the
 301         // Voronoi cell that it is within
 302         remap(ai,aj,ak,ci,cj,ck,x,y,z,ijk);
 303         vc.find_voronoi_cell(x,y,z,ci,cj,ck,ijk,w,mrs);
 304
 305         if(w.ijk!=-1) {
 306
 307                 // Assemble the position vector of the particle to be returned,
 308                 // applying a periodic remapping if necessary
 309                 ci+=w.di;if(ci<0||ci>=nx) ai+=step_div(ci,nx);
 310                 rx=p[w.ijk][3*w.l]+ak*bxz+aj*bxy+ai*bx;
 311                 ry=p[w.ijk][3*w.l+1]+ak*byz+aj*by;
 312                 rz=p[w.ijk][3*w.l+2]+ak*bz;
 313                 pid=id[w.ijk][w.l];
 314                 return true;
 315         }
 316         return false;
 317 }
 318
 319 /** Takes a vector and finds the particle whose Voronoi cell contains that
 320  * vector. Additional wall classes are not considered by this routine.
 321  * \param[in] (x,y,z) the vector to test.
 322  * \param[out] (rx,ry,rz) the position of the particle whose Voronoi cell
 323  *                        contains the vector. If the container is periodic,
 324  *                        this may point to a particle in a periodic image of
 325  *                        the primary domain.
 326  * \param[out] pid the ID of the particle.
 327  * \return True if a particle was found. If the container has no particles,
 328  * then the search will not find a Voronoi cell and false is returned. */
 329 bool container_periodic_poly::find_voronoi_cell(double x,double y,double z,double &rx,double &ry,double &rz,int &pid) {
 330         int ai,aj,ak,ci,cj,ck,ijk;
 331         particle_record w;
 332         double mrs;
 333
 334         // Remap the vector into the primary domain and then search for the
 335         // Voronoi cell that it is within
 336         remap(ai,aj,ak,ci,cj,ck,x,y,z,ijk);
 337         vc.find_voronoi_cell(x,y,z,ci,cj,ck,ijk,w,mrs);
 338
 339         if(w.ijk!=-1) {
 340
 341                 // Assemble the position vector of the particle to be returned,
 342                 // applying a periodic remapping if necessary
 343                 ci+=w.di;if(ci<0||ci>=nx) ai+=step_div(ci,nx);
 344                 rx=p[w.ijk][4*w.l]+ak*bxz+aj*bxy+ai*bx;
 345                 ry=p[w.ijk][4*w.l+1]+ak*byz+aj*by;
 346                 rz=p[w.ijk][4*w.l+2]+ak*bz;
 347                 pid=id[w.ijk][w.l];
 348                 return true;
 349         }
 350         return false;
 351 }
 352
 353 /** Increase memory for a particular region.
 354  * \param[in] i the index of the region to reallocate. */
 355 void container_periodic_base::add_particle_memory(int i) {
 356
 357         // Handle the case when no memory has been allocated for this block
 358         if(mem[i]==0) {
 359                 mem[i]=init_mem;
 360                 id[i]=new int[init_mem];
 361                 p[i]=new double[ps*init_mem];
 362                 return;
 363         }
 364
 365         // Otherwise, double the memory allocation for this block. Carry out a
 366         // check on the memory allocation size, and print a status message if
 367         // requested.
 368         int l,nmem(mem[i]<<1);
 369         if(nmem>max_particle_memory)
 370                 voro_fatal_error("Absolute maximum memory allocation exceeded",VOROPP_MEMORY_ERROR);
 371 #if VOROPP_VERBOSE >=3
 372         fprintf(stderr,"Particle memory in region %d scaled up to %d\n",i,nmem);
 373 #endif
 374
 375         // Allocate new memory and copy in the contents of the old arrays
 376         int *idp=new int[nmem];
 377         for(l=0;l<co[i];l++) idp[l]=id[i][l];
 378         double *pp=new double[ps*nmem];
 379         for(l=0;l<ps*co[i];l++) pp[l]=p[i][l];
 380
 381         // Update pointers and delete old arrays
 382         mem[i]=nmem;
 383         delete [] id[i];id[i]=idp;
 384         delete [] p[i];p[i]=pp;
 385 }
 386
 387 /** Import a list of particles from an open file stream into the container.
 388  * Entries of four numbers (Particle ID, x position, y position, z position)
 389  * are searched for. If the file cannot be successfully read, then the routine
 390  * causes a fatal error.
 391  * \param[in] fp the file handle to read from. */
 392 void container_periodic::import(FILE *fp) {
 393         int i,j;
 394         double x,y,z;
 395         while((j=fscanf(fp,"%d %lg %lg %lg",&i,&x,&y,&z))==4) put(i,x,y,z);
 396         if(j!=EOF) voro_fatal_error("File import error",VOROPP_FILE_ERROR);
 397 }
 398
 399 /** Import a list of particles from an open file stream, also storing the order
 400  * of that the particles are read. Entries of four numbers (Particle ID, x
 401  * position, y position, z position) are searched for. If the file cannot be
 402  * successfully read, then the routine causes a fatal error.
 403  * \param[in,out] vo a reference to an ordering class to use.
 404  * \param[in] fp the file handle to read from. */
 405 void container_periodic::import(particle_order &vo,FILE *fp) {
 406         int i,j;
 407         double x,y,z;
 408         while((j=fscanf(fp,"%d %lg %lg %lg",&i,&x,&y,&z))==4) put(vo,i,x,y,z);
 409         if(j!=EOF) voro_fatal_error("File import error",VOROPP_FILE_ERROR);
 410 }
 411
 412 /** Import a list of particles from an open file stream into the container.
 413  * Entries of five numbers (Particle ID, x position, y position, z position,
 414  * radius) are searched for. If the file cannot be successfully read, then the
 415  * routine causes a fatal error.
 416  * \param[in] fp the file handle to read from. */
 417 void container_periodic_poly::import(FILE *fp) {
 418         int i,j;
 419         double x,y,z,r;
 420         while((j=fscanf(fp,"%d %lg %lg %lg %lg",&i,&x,&y,&z,&r))==5) put(i,x,y,z,r);
 421         if(j!=EOF) voro_fatal_error("File import error",VOROPP_FILE_ERROR);
 422 }
 423
 424 /** Import a list of particles from an open file stream, also storing the order
 425  * of that the particles are read. Entries of four numbers (Particle ID, x
 426  * position, y position, z position, radius) are searched for. If the file
 427  * cannot be successfully read, then the routine causes a fatal error.
 428  * \param[in,out] vo a reference to an ordering class to use.
 429  * \param[in] fp the file handle to read from. */
 430 void container_periodic_poly::import(particle_order &vo,FILE *fp) {
 431         int i,j;
 432         double x,y,z,r;
 433         while((j=fscanf(fp,"%d %lg %lg %lg %lg",&i,&x,&y,&z,&r))==5) put(vo,i,x,y,z,r);
 434         if(j!=EOF) voro_fatal_error("File import error",VOROPP_FILE_ERROR);
 435 }
 436
 437 /** Outputs the a list of all the container regions along with the number of
 438  * particles stored within each. */
 439 void container_periodic_base::region_count() {
 440         int i,j,k,*cop=co;
 441         for(k=0;k<nz;k++) for(j=0;j<ny;j++) for(i=0;i<nx;i++)
 442                 printf("Region (%d,%d,%d): %d particles\n",i,j,k,*(cop++));
 443 }
 444
 445 /** Clears a container of particles. */
 446 void container_periodic::clear() {
 447         for(int *cop=co;cop<co+nxyz;cop++) *cop=0;
 448 }
 449
 450 /** Clears a container of particles, also clearing resetting the maximum radius
 451  * to zero. */
 452 void container_periodic_poly::clear() {
 453         for(int *cop=co;cop<co+nxyz;cop++) *cop=0;
 454         max_radius=0;
 455 }
 456
 457 /** Computes all the Voronoi cells and saves customized information about them.
 458  * \param[in] format the custom output string to use.
 459  * \param[in] fp a file handle to write to. */
 460 void container_periodic::print_custom(const char *format,FILE *fp) {
 461         c_loop_all_periodic vl(*this);
 462         print_custom(vl,format,fp);
 463 }
 464
 465 /** Computes all the Voronoi cells and saves customized
 466  * information about them.
 467  * \param[in] format the custom output string to use.
 468  * \param[in] fp a file handle to write to. */
 469 void container_periodic_poly::print_custom(const char *format,FILE *fp) {
 470         c_loop_all_periodic vl(*this);
 471         print_custom(vl,format,fp);
 472 }
 473
 474 /** Computes all the Voronoi cells and saves customized information about them.
 475  * \param[in] format the custom output string to use.
 476  * \param[in] filename the name of the file to write to. */
 477 void container_periodic::print_custom(const char *format,const char *filename) {
 478         FILE *fp=safe_fopen(filename,"w");
 479         print_custom(format,fp);
 480         fclose(fp);
 481 }
 482
 483 /** Computes all the Voronoi cells and saves customized
 484  * information about them
 485  * \param[in] format the custom output string to use.
 486  * \param[in] filename the name of the file to write to. */
 487 void container_periodic_poly::print_custom(const char *format,const char *filename) {
 488         FILE *fp=safe_fopen(filename,"w");
 489         print_custom(format,fp);
 490         fclose(fp);
 491 }
 492
 493 /** Computes all of the Voronoi cells in the container, but does nothing
 494  * with the output. It is useful for measuring the pure computation time
 495  * of the Voronoi algorithm, without any additional calculations such as
 496  * volume evaluation or cell output. */
 497 void container_periodic::compute_all_cells() {
 498         voronoicell c;
 499         c_loop_all_periodic vl(*this);
 500         if(vl.start()) do compute_cell(c,vl);
 501         while(vl.inc());
 502 }
 503
 504 /** Computes all of the Voronoi cells in the container, but does nothing
 505  * with the output. It is useful for measuring the pure computation time
 506  * of the Voronoi algorithm, without any additional calculations such as
 507  * volume evaluation or cell output. */
 508 void container_periodic_poly::compute_all_cells() {
 509         voronoicell c;
 510         c_loop_all_periodic vl(*this);
 511         if(vl.start()) do compute_cell(c,vl);while(vl.inc());
 512 }
 513
 514 /** Calculates all of the Voronoi cells and sums their volumes. In most cases
 515  * without walls, the sum of the Voronoi cell volumes should equal the volume
 516  * of the container to numerical precision.
 517  * \return The sum of all of the computed Voronoi volumes. */
 518 double container_periodic::sum_cell_volumes() {
 519         voronoicell c;
 520         double vol=0;
 521         c_loop_all_periodic vl(*this);
 522         if(vl.start()) do if(compute_cell(c,vl)) vol+=c.volume();while(vl.inc());
 523         return vol;
 524 }
 525
 526 /** Calculates all of the Voronoi cells and sums their volumes. In most cases
 527  * without walls, the sum of the Voronoi cell volumes should equal the volume
 528  * of the container to numerical precision.
 529  * \return The sum of all of the computed Voronoi volumes. */
 530 double container_periodic_poly::sum_cell_volumes() {
 531         voronoicell c;
 532         double vol=0;
 533         c_loop_all_periodic vl(*this);
 534         if(vl.start()) do if(compute_cell(c,vl)) vol+=c.volume();while(vl.inc());
 535         return vol;
 536 }
 537
 538 /** This routine creates all periodic images of the particles. It is meant for
 539  * diagnostic purposes only, since usually periodic images are dynamically
 540  * created in when they are referenced. */
 541 void container_periodic_base::create_all_images() {
 542         int i,j,k;
 543         for(k=0;k<oz;k++) for(j=0;j<oy;j++) for(i=0;i<nx;i++) create_periodic_image(i,j,k);
 544 }
 545
 546 /** Checks that the particles within each block lie within that block's bounds.
 547  * This is useful for diagnosing problems with periodic image computation. */
 548 void container_periodic_base::check_compartmentalized() {
 549         int c,l,i,j,k;
 550         double mix,miy,miz,max,may,maz,*pp;
 551         for(k=l=0;k<oz;k++) for(j=0;j<oy;j++) for(i=0;i<nx;i++,l++) if(mem[l]>0) {
 552
 553                 // Compute the block's bounds, adding in a small tolerance
 554                 mix=i*boxx-tolerance;max=mix+boxx+tolerance;
 555                 miy=(j-ey)*boxy-tolerance;may=miy+boxy+tolerance;
 556                 miz=(k-ez)*boxz-tolerance;maz=miz+boxz+tolerance;
 557
 558                 // Print entries for any particles that lie outside the block's
 559                 // bounds
 560                 for(pp=p[l],c=0;c<co[l];c++,pp+=ps) if(*pp<mix||*pp>max||pp[1]<miy||pp[1]>may||pp[2]<miz||pp[2]>maz)
 561                         printf("%d %d %d %d %f %f %f %f %f %f %f %f %f\n",
 562                                id[l][c],i,j,k,*pp,pp[1],pp[2],mix,max,miy,may,miz,maz);
 563         }
 564 }
 565
 566 /** Creates particles within an image block that is aligned with the primary
 567  * domain in the z axis. In this case, the image block may be comprised of
 568  * particles from two primary blocks. The routine considers these two primary
 569  * blocks, and adds the needed particles to the image. The remaining particles
 570  * from the primary blocks are also filled into the neighboring images.
 571  * \param[in] (di,dj,dk) the index of the block to consider. The z index must
 572  *                       satisfy ez<=dk<wz. */
 573 void container_periodic_base::create_side_image(int di,int dj,int dk) {
 574         int l,dijk=di+nx*(dj+oy*dk),odijk,ima=step_div(dj-ey,ny);
 575         int qua=di+step_int(-ima*bxy*xsp),quadiv=step_div(qua,nx);
 576         int fi=qua-quadiv*nx,fijk=fi+nx*(dj-ima*ny+oy*dk);
 577         double dis=ima*bxy+quadiv*bx,switchx=di*boxx-ima*bxy-quadiv*bx,adis;
 578
 579         // Left image computation
 580         if((img[dijk]&1)==0) {
 581                 if(di>0) {
 582                         odijk=dijk-1;adis=dis;
 583                 } else {
 584                         odijk=dijk+nx-1;adis=dis+bx;
 585                 }
 586                 img[odijk]|=2;
 587                 for(l=0;l<co[fijk];l++) {
 588                         if(p[fijk][ps*l]>switchx) put_image(dijk,fijk,l,dis,by*ima,0);
 589                         else put_image(odijk,fijk,l,adis,by*ima,0);
 590                 }
 591         }
 592
 593         // Right image computation
 594         if((img[dijk]&2)==0) {
 595                 if(fi==nx-1) {
 596                         fijk+=1-nx;switchx+=(1-nx)*boxx;dis+=bx;
 597                 } else {
 598                         fijk++;switchx+=boxx;
 599                 }
 600                 if(di==nx-1) {
 601                         odijk=dijk-nx+1;adis=dis-bx;
 602                 } else {
 603                         odijk=dijk+1;adis=dis;
 604                 }
 605                 img[odijk]|=1;
 606                 for(l=0;l<co[fijk];l++) {
 607                         if(p[fijk][ps*l]<switchx) put_image(dijk,fijk,l,dis,by*ima,0);
 608                         else put_image(odijk,fijk,l,adis,by*ima,0);
 609                 }
 610         }
 611
 612         // All contributions to the block now added, so set both two bits of
 613         // the image information
 614         img[dijk]=3;
 615 }
 616
 617 /** Creates particles within an image block that is not aligned with the
 618  * primary domain in the z axis. In this case, the image block may be comprised
 619  * of particles from four primary blocks. The routine considers these four
 620  * primary blocks, and adds the needed particles to the image. The remaining
 621  * particles from the primary blocks are also filled into the neighboring
 622  * images.
 623  * \param[in] (di,dj,dk) the index of the block to consider. The z index must
 624  *                       satisfy dk<ez or dk>=wz. */
 625 void container_periodic_base::create_vertical_image(int di,int dj,int dk) {
 626         int l,dijk=di+nx*(dj+oy*dk),dijkl,dijkr,ima=step_div(dk-ez,nz);
 627         int qj=dj+step_int(-ima*byz*ysp),qjdiv=step_div(qj-ey,ny);
 628         int qi=di+step_int((-ima*bxz-qjdiv*bxy)*xsp),qidiv=step_div(qi,nx);
 629         int fi=qi-qidiv*nx,fj=qj-qjdiv*ny,fijk=fi+nx*(fj+oy*(dk-ima*nz)),fijk2;
 630         double disy=ima*byz+qjdiv*by,switchy=(dj-ey)*boxy-ima*byz-qjdiv*by;
 631         double disx=ima*bxz+qjdiv*bxy+qidiv*bx,switchx=di*boxx-ima*bxz-qjdiv*bxy-qidiv*bx;
 632         double switchx2,disxl,disxr,disx2,disxr2;
 633
 634         if(di==0) {dijkl=dijk+nx-1;disxl=disx+bx;}
 635         else {dijkl=dijk-1;disxl=disx;}
 636
 637         if(di==nx-1) {dijkr=dijk-nx+1;disxr=disx-bx;}
 638         else {dijkr=dijk+1;disxr=disx;}
 639
 640         // Down-left image computation
 641         bool y_exist=dj!=0;
 642         if((img[dijk]&1)==0) {
 643                 img[dijkl]|=2;
 644                 if(y_exist) {
 645                         img[dijkl-nx]|=8;
 646                         img[dijk-nx]|=4;
 647                 }
 648                 for(l=0;l<co[fijk];l++) {
 649                         if(p[fijk][ps*l+1]>switchy) {
 650                                 if(p[fijk][ps*l]>switchx) put_image(dijk,fijk,l,disx,disy,bz*ima);
 651                                 else put_image(dijkl,fijk,l,disxl,disy,bz*ima);
 652                         } else {
 653                                 if(!y_exist) continue;
 654                                 if(p[fijk][ps*l]>switchx) put_image(dijk-nx,fijk,l,disx,disy,bz*ima);
 655                                 else put_image(dijkl-nx,fijk,l,disxl,disy,bz*ima);
 656                         }
 657                 }
 658         }
 659
 660         // Down-right image computation
 661         if((img[dijk]&2)==0) {
 662                 if(fi==nx-1) {
 663                         fijk2=fijk+1-nx;switchx2=switchx+(1-nx)*boxx;disx2=disx+bx;disxr2=disxr+bx;
 664                 } else {
 665                         fijk2=fijk+1;switchx2=switchx+boxx;disx2=disx;disxr2=disxr;
 666                 }
 667                 img[dijkr]|=1;
 668                 if(y_exist) {
 669                         img[dijkr-nx]|=4;
 670                         img[dijk-nx]|=8;
 671                 }
 672                 for(l=0;l<co[fijk2];l++) {
 673                         if(p[fijk2][ps*l+1]>switchy) {
 674                                 if(p[fijk2][ps*l]>switchx2) put_image(dijkr,fijk2,l,disxr2,disy,bz*ima);
 675                                 else put_image(dijk,fijk2,l,disx2,disy,bz*ima);
 676                         } else {
 677                                 if(!y_exist) continue;
 678                                 if(p[fijk2][ps*l]>switchx2) put_image(dijkr-nx,fijk2,l,disxr2,disy,bz*ima);
 679                                 else put_image(dijk-nx,fijk2,l,disx2,disy,bz*ima);
 680                         }
 681                 }
 682         }
 683
 684         // Recomputation of some intermediate quantities for boundary cases
 685         if(fj==wy-1) {
 686                 fijk+=nx*(1-ny)-fi;
 687                 switchy+=(1-ny)*boxy;
 688                 disy+=by;
 689                 qi=di+step_int(-(ima*bxz+(qjdiv+1)*bxy)*xsp);
 690                 int dqidiv=step_div(qi,nx)-qidiv;qidiv+=dqidiv;
 691                 fi=qi-qidiv*nx;
 692                 fijk+=fi;
 693                 disx+=bxy+bx*dqidiv;
 694                 disxl+=bxy+bx*dqidiv;
 695                 disxr+=bxy+bx*dqidiv;
 696                 switchx-=bxy+bx*dqidiv;
 697         } else {
 698                 fijk+=nx;switchy+=boxy;
 699         }
 700
 701         // Up-left image computation
 702         y_exist=dj!=oy-1;
 703         if((img[dijk]&4)==0) {
 704                 img[dijkl]|=8;
 705                 if(y_exist) {
 706                         img[dijkl+nx]|=2;
 707                         img[dijk+nx]|=1;
 708                 }
 709                 for(l=0;l<co[fijk];l++) {
 710                         if(p[fijk][ps*l+1]>switchy) {
 711                                 if(!y_exist) continue;
 712                                 if(p[fijk][ps*l]>switchx) put_image(dijk+nx,fijk,l,disx,disy,bz*ima);
 713                                 else put_image(dijkl+nx,fijk,l,disxl,disy,bz*ima);
 714                         } else {
 715                                 if(p[fijk][ps*l]>switchx) put_image(dijk,fijk,l,disx,disy,bz*ima);
 716                                 else put_image(dijkl,fijk,l,disxl,disy,bz*ima);
 717                         }
 718                 }
 719         }
 720
 721         // Up-right image computation
 722         if((img[dijk]&8)==0) {
 723                 if(fi==nx-1) {
 724                         fijk2=fijk+1-nx;switchx2=switchx+(1-nx)*boxx;disx2=disx+bx;disxr2=disxr+bx;
 725                 } else {
 726                         fijk2=fijk+1;switchx2=switchx+boxx;disx2=disx;disxr2=disxr;
 727                 }
 728                 img[dijkr]|=4;
 729                 if(y_exist) {
 730                         img[dijkr+nx]|=1;
 731                         img[dijk+nx]|=2;
 732                 }
 733                 for(l=0;l<co[fijk2];l++) {
 734                         if(p[fijk2][ps*l+1]>switchy) {
 735                                 if(!y_exist) continue;
 736                                 if(p[fijk2][ps*l]>switchx2) put_image(dijkr+nx,fijk2,l,disxr2,disy,bz*ima);
 737                                 else put_image(dijk+nx,fijk2,l,disx2,disy,bz*ima);
 738                         } else {
 739                                 if(p[fijk2][ps*l]>switchx2) put_image(dijkr,fijk2,l,disxr2,disy,bz*ima);
 740                                 else put_image(dijk,fijk2,l,disx2,disy,bz*ima);
 741                         }
 742                 }
 743         }
 744
 745         // All contributions to the block now added, so set all four bits of
 746         // the image information
 747         img[dijk]=15;
 748 }
 749
 750 /** Copies a particle position from the primary domain into an image block.
 751  * \param[in] reg the block index within the primary domain that the particle
 752  *                is within.
 753  * \param[in] fijk the index of the image block.
 754  * \param[in] l the index of the particle entry within the primary block.
 755  * \param[in] (dx,dy,dz) the displacement vector to add to the particle. */
 756 void container_periodic_base::put_image(int reg,int fijk,int l,double dx,double dy,double dz) {
 757         if(co[reg]==mem[reg]) add_particle_memory(reg);
 758         double *p1=p[reg]+ps*co[reg],*p2=p[fijk]+ps*l;
 759         *(p1++)=*(p2++)+dx;
 760         *(p1++)=*(p2++)+dy;
 761         *p1=*p2+dz;
 762         if(ps==4) *(++p1)=*(++p2);
 763         id[reg][co[reg]++]=id[fijk][l];
 764 }
 765
 766 }