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