trunk/src/v_compute.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 v_compute.cc
   8  * \brief Function implementantions for the voro_compute template. */
   9
  10 #include "worklist.hh"
  11 #include "v_compute.hh"
  12 #include "rad_option.hh"
  13 #include "container.hh"
  14 #include "container_prd.hh"
  15
  16 namespace voro {
  17
  18 /** The class constructor initializes constants from the container class, and
  19  * sets up the mask and queue used for Voronoi computations.
  20  * \param[in] con_ a reference to the container class to use.
  21  * \param[in] (hx_,hy_,hz_) the size of the mask to use. */
  22 template<class c_class>
  23 voro_compute<c_class>::voro_compute(c_class &con_,int hx_,int hy_,int hz_) :
  24         con(con_), boxx(con_.boxx), boxy(con_.boxy), boxz(con_.boxz),
  25         xsp(con_.xsp), ysp(con_.ysp), zsp(con_.zsp),
  26         hx(hx_), hy(hy_), hz(hz_), hxy(hx_*hy_), hxyz(hxy*hz_), ps(con_.ps),
  27         id(con_.id), p(con_.p), co(con_.co), bxsq(boxx*boxx+boxy*boxy+boxz*boxz),
  28         mv(0), qu_size(3*(3+hxy+hz*(hx+hy))), wl(con_.wl), mrad(con_.mrad),
  29         mask(new unsigned int[hxyz]), qu(new int[qu_size]), qu_l(qu+qu_size) {
  30         reset_mask();
  31 }
  32
  33 /** Scans all of the particles within a block to see if any of them have a
  34  * smaller distance to the given test vector. If one is found, the routine
  35  * updates the minimum distance and store information about this particle.
  36  * \param[in] ijk the index of the block.
  37  * \param[in] (x,y,z) the test vector to consider (which may have already had a
  38  *                    periodic displacement applied to it).
  39  * \param[in] (di,dj,dk) the coordinates of the current block, to store if the
  40  *                       particle record is updated.
  41  * \param[in,out] w a reference to a particle record in which to store
  42  *                  information about the particle whose Voronoi cell the
  43  *                  vector is within.
  44  * \param[in,out] mrs the current minimum distance, that may be updated if a
  45  *                    closer particle is found. */
  46 template<class c_class>
  47 inline void voro_compute<c_class>::scan_all(int ijk,double x,double y,double z,int di,int dj,int dk,particle_record &w,double &mrs) {
  48         double x1,y1,z1,rs;bool in_block=false;
  49         for(int l=0;l<co[ijk];l++) {
  50                 x1=p[ijk][ps*l]-x;
  51                 y1=p[ijk][ps*l+1]-y;
  52                 z1=p[ijk][ps*l+2]-z;
  53                 rs=con.r_current_sub(x1*x1+y1*y1+z1*z1,ijk,l);
  54                 if(rs<mrs) {mrs=rs;w.l=l;in_block=true;}
  55         }
  56         if(in_block) {w.ijk=ijk;w.di=di;w.dj=dj,w.dk=dk;}
  57 }
  58
  59 /** Finds the Voronoi cell that given vector is within. For containers that are
  60  * not radially dependent, this corresponds to findig the particle that is
  61  * closest to the vector; for the radical tessellation containers, this
  62  * corresponds to a finding the minimum weighted distance.
  63  * \param[in] (x,y,z) the vector to consider.
  64  * \param[in] (ci,cj,ck) the coordinates of the block that the test particle is
  65  *                       in relative to the container data structure.
  66  * \param[in] ijk the index of the block that the test particle is in.
  67  * \param[out] w a reference to a particle record in which to store information
  68  *               about the particle whose Voronoi cell the vector is within.
  69  * \param[out] mrs the minimum computed distance. */
  70 template<class c_class>
  71 void voro_compute<c_class>::find_voronoi_cell(double x,double y,double z,int ci,int cj,int ck,int ijk,particle_record &w,double &mrs) {
  72         double qx=0,qy=0,qz=0,rs;
  73         int i,j,k,di,dj,dk,ei,ej,ek,f,g,disp;
  74         double fx,fy,fz,mxs,mys,mzs,*radp;
  75         unsigned int q,*e,*mijk;
  76
  77         // Init setup for parameters to return
  78         w.ijk=-1;mrs=large_number;
  79
  80         con.initialize_search(ci,cj,ck,ijk,i,j,k,disp);
  81
  82         // Test all particles in the particle's local region first
  83         scan_all(ijk,x,y,z,0,0,0,w,mrs);
  84
  85         // Now compute the fractional position of the particle within its
  86         // region and store it in (fx,fy,fz). We use this to compute an index
  87         // (di,dj,dk) of which subregion the particle is within.
  88         unsigned int m1,m2;
  89         con.frac_pos(x,y,z,ci,cj,ck,fx,fy,fz);
  90         di=int(fx*xsp*wl_fgrid);dj=int(fy*ysp*wl_fgrid);dk=int(fz*zsp*wl_fgrid);
  91
  92         // The indices (di,dj,dk) tell us which worklist to use, to test the
  93         // blocks in the optimal order. But we only store worklists for the
  94         // eighth of the region where di, dj, and dk are all less than half the
  95         // full grid. The rest of the cases are handled by symmetry. In this
  96         // section, we detect for these cases, by reflecting high values of di,
  97         // dj, and dk. For these cases, a mask is constructed in m1 and m2
  98         // which is used to flip the worklist information when it is loaded.
  99         if(di>=wl_hgrid) {
 100                 mxs=boxx-fx;
 101                 m1=127+(3<<21);m2=1+(1<<21);di=wl_fgrid-1-di;if(di<0) di=0;
 102         } else {m1=m2=0;mxs=fx;}
 103         if(dj>=wl_hgrid) {
 104                 mys=boxy-fy;
 105                 m1|=(127<<7)+(3<<24);m2|=(1<<7)+(1<<24);dj=wl_fgrid-1-dj;if(dj<0) dj=0;
 106         } else mys=fy;
 107         if(dk>=wl_hgrid) {
 108                 mzs=boxz-fz;
 109                 m1|=(127<<14)+(3<<27);m2|=(1<<14)+(1<<27);dk=wl_fgrid-1-dk;if(dk<0) dk=0;
 110         } else mzs=fz;
 111
 112         // Do a quick test to account for the case when the minimum radius is
 113         // small enought that no other blocks need to be considered
 114         rs=con.r_max_add(mrs);
 115         if(mxs*mxs>rs&&mys*mys>rs&&mzs*mzs>rs) return;
 116
 117         // Now compute which worklist we are going to use, and set radp and e to
 118         // point at the right offsets
 119         ijk=di+wl_hgrid*(dj+wl_hgrid*dk);
 120         radp=mrad+ijk*wl_seq_length;
 121         e=(const_cast<unsigned int*> (wl))+ijk*wl_seq_length;
 122
 123         // Read in how many items in the worklist can be tested without having to
 124         // worry about writing to the mask
 125         f=e[0];g=0;
 126         do {
 127
 128                 // If mrs is less than the minimum distance to any untested
 129                 // block, then we are done
 130                 if(con.r_max_add(mrs)<radp[g]) return;
 131                 g++;
 132
 133                 // Load in a block off the worklist, permute it with the
 134                 // symmetry mask, and decode its position. These are all
 135                 // integer bit operations so they should run very fast.
 136                 q=e[g];q^=m1;q+=m2;
 137                 di=q&127;di-=64;
 138                 dj=(q>>7)&127;dj-=64;
 139                 dk=(q>>14)&127;dk-=64;
 140
 141                 // Check that the worklist position is in range
 142                 ei=di+i;if(ei<0||ei>=hx) continue;
 143                 ej=dj+j;if(ej<0||ej>=hy) continue;
 144                 ek=dk+k;if(ek<0||ek>=hz) continue;
 145
 146                 // Call the compute_min_max_radius() function. This returns
 147                 // true if the minimum distance to the block is bigger than the
 148                 // current mrs, in which case we skip this block and move on.
 149                 // Otherwise, it computes the maximum distance to the block and
 150                 // returns it in crs.
 151                 if(compute_min_radius(di,dj,dk,fx,fy,fz,mrs)) continue;
 152
 153                 // Now compute which region we are going to loop over, adding a
 154                 // displacement for the periodic cases
 155                 ijk=con.region_index(ci,cj,ck,ei,ej,ek,qx,qy,qz,disp);
 156
 157                 // If mrs is bigger than the maximum distance to the block,
 158                 // then we have to test all particles in the block for
 159                 // intersections. Otherwise, we do additional checks and skip
 160                 // those particles which can't possibly intersect the block.
 161                 scan_all(ijk,x-qx,y-qy,z-qz,di,dj,dk,w,mrs);
 162         } while(g<f);
 163
 164         // Update mask value and initialize queue
 165         mv++;
 166         if(mv==0) {reset_mask();mv=1;}
 167         int *qu_s=qu,*qu_e=qu;
 168
 169         while(g<wl_seq_length-1) {
 170
 171                 // If mrs is less than the minimum distance to any untested
 172                 // block, then we are done
 173                 if(con.r_max_add(mrs)<radp[g]) return;
 174                 g++;
 175
 176                 // Load in a block off the worklist, permute it with the
 177                 // symmetry mask, and decode its position. These are all
 178                 // integer bit operations so they should run very fast.
 179                 q=e[g];q^=m1;q+=m2;
 180                 di=q&127;di-=64;
 181                 dj=(q>>7)&127;dj-=64;
 182                 dk=(q>>14)&127;dk-=64;
 183
 184                 // Compute the position in the mask of the current block. If
 185                 // this lies outside the mask, then skip it. Otherwise, mark
 186                 // it.
 187                 ei=di+i;if(ei<0||ei>=hx) continue;
 188                 ej=dj+j;if(ej<0||ej>=hy) continue;
 189                 ek=dk+k;if(ek<0||ek>=hz) continue;
 190                 mijk=mask+ei+hx*(ej+hy*ek);
 191                 *mijk=mv;
 192
 193                 // Skip this block if it is further away than the current
 194                 // minimum radius
 195                 if(compute_min_radius(di,dj,dk,fx,fy,fz,mrs)) continue;
 196
 197                 // Now compute which region we are going to loop over, adding a
 198                 // displacement for the periodic cases
 199                 ijk=con.region_index(ci,cj,ck,ei,ej,ek,qx,qy,qz,disp);
 200                 scan_all(ijk,x-qx,y-qy,z-qz,di,dj,dk,w,mrs);
 201
 202                 if(qu_e>qu_l-18) add_list_memory(qu_s,qu_e);
 203                 scan_bits_mask_add(q,mijk,ei,ej,ek,qu_e);
 204         }
 205
 206         // Do a check to see if we've reached the radius cutoff
 207         if(con.r_max_add(mrs)<radp[g]) return;
 208
 209         // We were unable to completely compute the cell based on the blocks in
 210         // the worklist, so now we have to go block by block, reading in items
 211         // off the list
 212         while(qu_s!=qu_e) {
 213
 214                 // Read the next entry of the queue
 215                 if(qu_s==qu_l) qu_s=qu;
 216                 ei=*(qu_s++);ej=*(qu_s++);ek=*(qu_s++);
 217                 di=ei-i;dj=ej-j;dk=ek-k;
 218                 if(compute_min_radius(di,dj,dk,fx,fy,fz,mrs)) continue;
 219
 220                 ijk=con.region_index(ci,cj,ck,ei,ej,ek,qx,qy,qz,disp);
 221                 scan_all(ijk,x-qx,y-qy,z-qz,di,dj,dk,w,mrs);
 222
 223                 // Test the neighbors of the current block, and add them to the
 224                 // block list if they haven't already been tested
 225                 if((qu_s<=qu_e?(qu_l-qu_e)+(qu_s-qu):qu_s-qu_e)<18) add_list_memory(qu_s,qu_e);
 226                 add_to_mask(ei,ej,ek,qu_e);
 227         }
 228 }
 229
 230 /** Scans the six orthogonal neighbors of a given block and adds them to the
 231  * queue if they haven't been considered already. It assumes that the queue
 232  * will definitely have enough memory to add six entries at the end.
 233  * \param[in] (ei,ej,ek) the block to consider.
 234  * \param[in,out] qu_e a pointer to the end of the queue. */
 235 template<class c_class>
 236 inline void voro_compute<c_class>::add_to_mask(int ei,int ej,int ek,int *&qu_e) {
 237         unsigned int *mijk=mask+ei+hx*(ej+hy*ek);
 238         if(ek>0) if(*(mijk-hxy)!=mv) {if(qu_e==qu_l) qu_e=qu;*(mijk-hxy)=mv;*(qu_e++)=ei;*(qu_e++)=ej;*(qu_e++)=ek-1;}
 239         if(ej>0) if(*(mijk-hx)!=mv) {if(qu_e==qu_l) qu_e=qu;*(mijk-hx)=mv;*(qu_e++)=ei;*(qu_e++)=ej-1;*(qu_e++)=ek;}
 240         if(ei>0) if(*(mijk-1)!=mv) {if(qu_e==qu_l) qu_e=qu;*(mijk-1)=mv;*(qu_e++)=ei-1;*(qu_e++)=ej;*(qu_e++)=ek;}
 241         if(ei<hx-1) if(*(mijk+1)!=mv) {if(qu_e==qu_l) qu_e=qu;*(mijk+1)=mv;*(qu_e++)=ei+1;*(qu_e++)=ej;*(qu_e++)=ek;}
 242         if(ej<hy-1) if(*(mijk+hx)!=mv) {if(qu_e==qu_l) qu_e=qu;*(mijk+hx)=mv;*(qu_e++)=ei;*(qu_e++)=ej+1;*(qu_e++)=ek;}
 243         if(ek<hz-1) if(*(mijk+hxy)!=mv) {if(qu_e==qu_l) qu_e=qu;*(mijk+hxy)=mv;*(qu_e++)=ei;*(qu_e++)=ej;*(qu_e++)=ek+1;}
 244 }
 245
 246 /** Scans a worklist entry and adds any blocks to the queue
 247  * \param[in] (ei,ej,ek) the block to consider.
 248  * \param[in,out] qu_e a pointer to the end of the queue. */
 249 template<class c_class>
 250 inline void voro_compute<c_class>::scan_bits_mask_add(unsigned int q,unsigned int *mijk,int ei,int ej,int ek,int *&qu_e) {
 251         const unsigned int b1=1<<21,b2=1<<22,b3=1<<24,b4=1<<25,b5=1<<27,b6=1<<28;
 252         if((q&b2)==b2) {
 253                 if(ei>0) {*(mijk-1)=mv;*(qu_e++)=ei-1;*(qu_e++)=ej;*(qu_e++)=ek;}
 254                 if((q&b1)==0&&ei<hx-1) {*(mijk+1)=mv;*(qu_e++)=ei+1;*(qu_e++)=ej;*(qu_e++)=ek;}
 255         } else if((q&b1)==b1&&ei<hx-1) {*(mijk+1)=mv;*(qu_e++)=ei+1;*(qu_e++)=ej;*(qu_e++)=ek;}
 256         if((q&b4)==b4) {
 257                 if(ej>0) {*(mijk-hx)=mv;*(qu_e++)=ei;*(qu_e++)=ej-1;*(qu_e++)=ek;}
 258                 if((q&b3)==0&&ej<hy-1) {*(mijk+hx)=mv;*(qu_e++)=ei;*(qu_e++)=ej+1;*(qu_e++)=ek;}
 259         } else if((q&b3)==b3&&ej<hy-1) {*(mijk+hx)=mv;*(qu_e++)=ei;*(qu_e++)=ej+1;*(qu_e++)=ek;}
 260         if((q&b6)==b6) {
 261                 if(ek>0) {*(mijk-hxy)=mv;*(qu_e++)=ei;*(qu_e++)=ej;*(qu_e++)=ek-1;}
 262                 if((q&b5)==0&&ek<hz-1) {*(mijk+hxy)=mv;*(qu_e++)=ei;*(qu_e++)=ej;*(qu_e++)=ek+1;}
 263         } else if((q&b5)==b5&&ek<hz-1) {*(mijk+hxy)=mv;*(qu_e++)=ei;*(qu_e++)=ej;*(qu_e++)=ek+1;}
 264 }
 265
 266 /** This routine computes a Voronoi cell for a single particle in the
 267  * container. It can be called by the user, but is also forms the core part of
 268  * several of the main functions, such as store_cell_volumes(), print_all(),
 269  * and the drawing routines. The algorithm constructs the cell by testing over
 270  * the neighbors of the particle, working outwards until it reaches those
 271  * particles which could not possibly intersect the cell. For maximum
 272  * efficiency, this algorithm is divided into three parts. In the first
 273  * section, the algorithm tests over the blocks which are in the immediate
 274  * vicinity of the particle, by making use of one of the precomputed worklists.
 275  * The code then continues to test blocks on the worklist, but also begins to
 276  * construct a list of neighboring blocks outside the worklist which may need
 277  * to be test. In the third section, the routine starts testing these
 278  * neighboring blocks, evaluating whether or not a particle in them could
 279  * possibly intersect the cell. For blocks that intersect the cell, it tests
 280  * the particles in that block, and then adds the block neighbors to the list
 281  * of potential places to consider.
 282  * \param[in,out] c a reference to a voronoicell object.
 283  * \param[in] ijk the index of the block that the test particle is in.
 284  * \param[in] s the index of the particle within the test block.
 285  * \param[in] (ci,cj,ck) the coordinates of the block that the test particle is
 286  *                       in relative to the container data structure.
 287  * \return False if the Voronoi cell was completely removed during the
 288  *         computation and has zero volume, true otherwise. */
 289 template<class c_class>
 290 template<class v_cell>
 291 bool voro_compute<c_class>::compute_cell(v_cell &c,int ijk,int s,int ci,int cj,int ck) {
 292         static const int count_list[8]={7,11,15,19,26,35,45,59},*count_e=count_list+8;
 293         double x,y,z,x1,y1,z1,qx=0,qy=0,qz=0;
 294         double xlo,ylo,zlo,xhi,yhi,zhi,x2,y2,z2,rs;
 295         int i,j,k,di,dj,dk,ei,ej,ek,f,g,l,disp;
 296         double fx,fy,fz,gxs,gys,gzs,*radp;
 297         unsigned int q,*e,*mijk;
 298
 299         if(!con.initialize_voronoicell(c,ijk,s,ci,cj,ck,i,j,k,x,y,z,disp)) return false;
 300         con.r_init(ijk,s);
 301
 302         // Initialize the Voronoi cell to fill the entire container
 303         double crs,mrs;
 304
 305         int next_count=3,*count_p=(const_cast<int*> (count_list));
 306
 307         // Test all particles in the particle's local region first
 308         for(l=0;l<s;l++) {
 309                 x1=p[ijk][ps*l]-x;
 310                 y1=p[ijk][ps*l+1]-y;
 311                 z1=p[ijk][ps*l+2]-z;
 312                 rs=con.r_scale(x1*x1+y1*y1+z1*z1,ijk,l);
 313                 if(!c.nplane(x1,y1,z1,rs,id[ijk][l])) return false;
 314         }
 315         l++;
 316         while(l<co[ijk]) {
 317                 x1=p[ijk][ps*l]-x;
 318                 y1=p[ijk][ps*l+1]-y;
 319                 z1=p[ijk][ps*l+2]-z;
 320                 rs=con.r_scale(x1*x1+y1*y1+z1*z1,ijk,l);
 321                 if(!c.nplane(x1,y1,z1,rs,id[ijk][l])) return false;
 322                 l++;
 323         }
 324
 325         // Now compute the maximum distance squared from the cell center to a
 326         // vertex. This is used to cut off the calculation since we only need
 327         // to test out to twice this range.
 328         mrs=c.max_radius_squared();
 329
 330         // Now compute the fractional position of the particle within its
 331         // region and store it in (fx,fy,fz). We use this to compute an index
 332         // (di,dj,dk) of which subregion the particle is within.
 333         unsigned int m1,m2;
 334         con.frac_pos(x,y,z,ci,cj,ck,fx,fy,fz);
 335         di=int(fx*xsp*wl_fgrid);dj=int(fy*ysp*wl_fgrid);dk=int(fz*zsp*wl_fgrid);
 336
 337         // The indices (di,dj,dk) tell us which worklist to use, to test the
 338         // blocks in the optimal order. But we only store worklists for the
 339         // eighth of the region where di, dj, and dk are all less than half the
 340         // full grid. The rest of the cases are handled by symmetry. In this
 341         // section, we detect for these cases, by reflecting high values of di,
 342         // dj, and dk. For these cases, a mask is constructed in m1 and m2
 343         // which is used to flip the worklist information when it is loaded.
 344         if(di>=wl_hgrid) {
 345                 gxs=fx;
 346                 m1=127+(3<<21);m2=1+(1<<21);di=wl_fgrid-1-di;if(di<0) di=0;
 347         } else {m1=m2=0;gxs=boxx-fx;}
 348         if(dj>=wl_hgrid) {
 349                 gys=fy;
 350                 m1|=(127<<7)+(3<<24);m2|=(1<<7)+(1<<24);dj=wl_fgrid-1-dj;if(dj<0) dj=0;
 351         } else gys=boxy-fy;
 352         if(dk>=wl_hgrid) {
 353                 gzs=fz;
 354                 m1|=(127<<14)+(3<<27);m2|=(1<<14)+(1<<27);dk=wl_fgrid-1-dk;if(dk<0) dk=0;
 355         } else gzs=boxz-fz;
 356         gxs*=gxs;gys*=gys;gzs*=gzs;
 357
 358         // Now compute which worklist we are going to use, and set radp and e to
 359         // point at the right offsets
 360         ijk=di+wl_hgrid*(dj+wl_hgrid*dk);
 361         radp=mrad+ijk*wl_seq_length;
 362         e=(const_cast<unsigned int*> (wl))+ijk*wl_seq_length;
 363
 364         // Read in how many items in the worklist can be tested without having to
 365         // worry about writing to the mask
 366         f=e[0];g=0;
 367         do {
 368
 369                 // At the intervals specified by count_list, we recompute the
 370                 // maximum radius squared
 371                 if(g==next_count) {
 372                         mrs=c.max_radius_squared();
 373                         if(count_p!=count_e) next_count=*(count_p++);
 374                 }
 375
 376                 // If mrs is less than the minimum distance to any untested
 377                 // block, then we are done
 378                 if(con.r_ctest(radp[g],mrs)) return true;
 379                 g++;
 380
 381                 // Load in a block off the worklist, permute it with the
 382                 // symmetry mask, and decode its position. These are all
 383                 // integer bit operations so they should run very fast.
 384                 q=e[g];q^=m1;q+=m2;
 385                 di=q&127;di-=64;
 386                 dj=(q>>7)&127;dj-=64;
 387                 dk=(q>>14)&127;dk-=64;
 388
 389                 // Check that the worklist position is in range
 390                 ei=di+i;if(ei<0||ei>=hx) continue;
 391                 ej=dj+j;if(ej<0||ej>=hy) continue;
 392                 ek=dk+k;if(ek<0||ek>=hz) continue;
 393
 394                 // Call the compute_min_max_radius() function. This returns
 395                 // true if the minimum distance to the block is bigger than the
 396                 // current mrs, in which case we skip this block and move on.
 397                 // Otherwise, it computes the maximum distance to the block and
 398                 // returns it in crs.
 399                 if(compute_min_max_radius(di,dj,dk,fx,fy,fz,gxs,gys,gzs,crs,mrs)) continue;
 400
 401                 // Now compute which region we are going to loop over, adding a
 402                 // displacement for the periodic cases
 403                 ijk=con.region_index(ci,cj,ck,ei,ej,ek,qx,qy,qz,disp);
 404
 405                 // If mrs is bigger than the maximum distance to the block,
 406                 // then we have to test all particles in the block for
 407                 // intersections. Otherwise, we do additional checks and skip
 408                 // those particles which can't possibly intersect the block.
 409                 if(co[ijk]>0) {
 410                         l=0;x2=x-qx;y2=y-qy;z2=z-qz;
 411                         if(!con.r_ctest(crs,mrs)) {
 412                                 do {
 413                                         x1=p[ijk][ps*l]-x2;
 414                                         y1=p[ijk][ps*l+1]-y2;
 415                                         z1=p[ijk][ps*l+2]-z2;
 416                                         rs=con.r_scale(x1*x1+y1*y1+z1*z1,ijk,l);
 417                                         if(!c.nplane(x1,y1,z1,rs,id[ijk][l])) return false;
 418                                         l++;
 419                                 } while (l<co[ijk]);
 420                         } else {
 421                                 do {
 422                                         x1=p[ijk][ps*l]-x2;
 423                                         y1=p[ijk][ps*l+1]-y2;
 424                                         z1=p[ijk][ps*l+2]-z2;
 425                                         rs=x1*x1+y1*y1+z1*z1;
 426                                         if(con.r_scale_check(rs,mrs,ijk,l)&&!c.nplane(x1,y1,z1,rs,id[ijk][l])) return false;
 427                                         l++;
 428                                 } while (l<co[ijk]);
 429                         }
 430                 }
 431         } while(g<f);
 432
 433         // If we reach here, we were unable to compute the entire cell using
 434         // the first part of the worklist. This section of the algorithm
 435         // continues the worklist, but it now starts preparing the mask that we
 436         // need if we end up going block by block. We do the same as before,
 437         // but we put a mark down on the mask for every block that's tested.
 438         // The worklist also contains information about which neighbors of each
 439         // block are not also on the worklist, and we start storing those
 440         // points in a list in case we have to go block by block. Update the
 441         // mask counter, and if it wraps around then reset the whole mask; that
 442         // will only happen once every 2^32 tries.
 443         mv++;
 444         if(mv==0) {reset_mask();mv=1;}
 445
 446         // Set the queue pointers
 447         int *qu_s=qu,*qu_e=qu;
 448
 449         while(g<wl_seq_length-1) {
 450
 451                 // At the intervals specified by count_list, we recompute the
 452                 // maximum radius squared
 453                 if(g==next_count) {
 454                         mrs=c.max_radius_squared();
 455                         if(count_p!=count_e) next_count=*(count_p++);
 456                 }
 457
 458                 // If mrs is less than the minimum distance to any untested
 459                 // block, then we are done
 460                 if(con.r_ctest(radp[g],mrs)) return true;
 461                 g++;
 462
 463                 // Load in a block off the worklist, permute it with the
 464                 // symmetry mask, and decode its position. These are all
 465                 // integer bit operations so they should run very fast.
 466                 q=e[g];q^=m1;q+=m2;
 467                 di=q&127;di-=64;
 468                 dj=(q>>7)&127;dj-=64;
 469                 dk=(q>>14)&127;dk-=64;
 470
 471                 // Compute the position in the mask of the current block. If
 472                 // this lies outside the mask, then skip it. Otherwise, mark
 473                 // it.
 474                 ei=di+i;if(ei<0||ei>=hx) continue;
 475                 ej=dj+j;if(ej<0||ej>=hy) continue;
 476                 ek=dk+k;if(ek<0||ek>=hz) continue;
 477                 mijk=mask+ei+hx*(ej+hy*ek);
 478                 *mijk=mv;
 479
 480                 // Call the compute_min_max_radius() function. This returns
 481                 // true if the minimum distance to the block is bigger than the
 482                 // current mrs, in which case we skip this block and move on.
 483                 // Otherwise, it computes the maximum distance to the block and
 484                 // returns it in crs.
 485                 if(compute_min_max_radius(di,dj,dk,fx,fy,fz,gxs,gys,gzs,crs,mrs)) continue;
 486
 487                 // Now compute which region we are going to loop over, adding a
 488                 // displacement for the periodic cases
 489                 ijk=con.region_index(ci,cj,ck,ei,ej,ek,qx,qy,qz,disp);
 490
 491                 // If mrs is bigger than the maximum distance to the block,
 492                 // then we have to test all particles in the block for
 493                 // intersections. Otherwise, we do additional checks and skip
 494                 // those particles which can't possibly intersect the block.
 495                 if(co[ijk]>0) {
 496                         l=0;x2=x-qx;y2=y-qy;z2=z-qz;
 497                         if(!con.r_ctest(crs,mrs)) {
 498                                 do {
 499                                         x1=p[ijk][ps*l]-x2;
 500                                         y1=p[ijk][ps*l+1]-y2;
 501                                         z1=p[ijk][ps*l+2]-z2;
 502                                         rs=con.r_scale(x1*x1+y1*y1+z1*z1,ijk,l);
 503                                         if(!c.nplane(x1,y1,z1,rs,id[ijk][l])) return false;
 504                                         l++;
 505                                 } while (l<co[ijk]);
 506                         } else {
 507                                 do {
 508                                         x1=p[ijk][ps*l]-x2;
 509                                         y1=p[ijk][ps*l+1]-y2;
 510                                         z1=p[ijk][ps*l+2]-z2;
 511                                         rs=x1*x1+y1*y1+z1*z1;
 512                                         if(con.r_scale_check(rs,mrs,ijk,l)&&!c.nplane(x1,y1,z1,rs,id[ijk][l])) return false;
 513                                         l++;
 514                                 } while (l<co[ijk]);
 515                         }
 516                 }
 517
 518                 // If there might not be enough memory on the list for these
 519                 // additions, then add more
 520                 if(qu_e>qu_l-18) add_list_memory(qu_s,qu_e);
 521
 522                 // Test the parts of the worklist element which tell us what
 523                 // neighbors of this block are not on the worklist. Store them
 524                 // on the block list, and mark the mask.
 525                 scan_bits_mask_add(q,mijk,ei,ej,ek,qu_e);
 526         }
 527
 528         // Do a check to see if we've reached the radius cutoff
 529         if(con.r_ctest(radp[g],mrs)) return true;
 530
 531         // We were unable to completely compute the cell based on the blocks in
 532         // the worklist, so now we have to go block by block, reading in items
 533         // off the list
 534         while(qu_s!=qu_e) {
 535
 536                 // If we reached the end of the list memory loop back to the
 537                 // start
 538                 if(qu_s==qu_l) qu_s=qu;
 539
 540                 // Read in a block off the list, and compute the upper and lower
 541                 // coordinates in each of the three dimensions
 542                 ei=*(qu_s++);ej=*(qu_s++);ek=*(qu_s++);
 543                 xlo=(ei-i)*boxx-fx;xhi=xlo+boxx;
 544                 ylo=(ej-j)*boxy-fy;yhi=ylo+boxy;
 545                 zlo=(ek-k)*boxz-fz;zhi=zlo+boxz;
 546
 547                 // Carry out plane tests to see if any particle in this block
 548                 // could possibly intersect the cell
 549                 if(ei>i) {
 550                         if(ej>j) {
 551                                 if(ek>k) {if(corner_test(c,xlo,ylo,zlo,xhi,yhi,zhi)) continue;}
 552                                 else if(ek<k) {if(corner_test(c,xlo,ylo,zhi,xhi,yhi,zlo)) continue;}
 553                                 else {if(edge_z_test(c,xlo,ylo,zlo,xhi,yhi,zhi)) continue;}
 554                         } else if(ej<j) {
 555                                 if(ek>k) {if(corner_test(c,xlo,yhi,zlo,xhi,ylo,zhi)) continue;}
 556                                 else if(ek<k) {if(corner_test(c,xlo,yhi,zhi,xhi,ylo,zlo)) continue;}
 557                                 else {if(edge_z_test(c,xlo,yhi,zlo,xhi,ylo,zhi)) continue;}
 558                         } else {
 559                                 if(ek>k) {if(edge_y_test(c,xlo,ylo,zlo,xhi,yhi,zhi)) continue;}
 560                                 else if(ek<k) {if(edge_y_test(c,xlo,ylo,zhi,xhi,yhi,zlo)) continue;}
 561                                 else {if(face_x_test(c,xlo,ylo,zlo,yhi,zhi)) continue;}
 562                         }
 563                 } else if(ei<i) {
 564                         if(ej>j) {
 565                                 if(ek>k) {if(corner_test(c,xhi,ylo,zlo,xlo,yhi,zhi)) continue;}
 566                                 else if(ek<k) {if(corner_test(c,xhi,ylo,zhi,xlo,yhi,zlo)) continue;}
 567                                 else {if(edge_z_test(c,xhi,ylo,zlo,xlo,yhi,zhi)) continue;}
 568                         } else if(ej<j) {
 569                                 if(ek>k) {if(corner_test(c,xhi,yhi,zlo,xlo,ylo,zhi)) continue;}
 570                                 else if(ek<k) {if(corner_test(c,xhi,yhi,zhi,xlo,ylo,zlo)) continue;}
 571                                 else {if(edge_z_test(c,xhi,yhi,zlo,xlo,ylo,zhi)) continue;}
 572                         } else {
 573                                 if(ek>k) {if(edge_y_test(c,xhi,ylo,zlo,xlo,yhi,zhi)) continue;}
 574                                 else if(ek<k) {if(edge_y_test(c,xhi,ylo,zhi,xlo,yhi,zlo)) continue;}
 575                                 else {if(face_x_test(c,xhi,ylo,zlo,yhi,zhi)) continue;}
 576                         }
 577                 } else {
 578                         if(ej>j) {
 579                                 if(ek>k) {if(edge_x_test(c,xlo,ylo,zlo,xhi,yhi,zhi)) continue;}
 580                                 else if(ek<k) {if(edge_x_test(c,xlo,ylo,zhi,xhi,yhi,zlo)) continue;}
 581                                 else {if(face_y_test(c,xlo,ylo,zlo,xhi,zhi)) continue;}
 582                         } else if(ej<j) {
 583                                 if(ek>k) {if(edge_x_test(c,xlo,yhi,zlo,xhi,ylo,zhi)) continue;}
 584                                 else if(ek<k) {if(edge_x_test(c,xlo,yhi,zhi,xhi,ylo,zlo)) continue;}
 585                                 else {if(face_y_test(c,xlo,yhi,zlo,xhi,zhi)) continue;}
 586                         } else {
 587                                 if(ek>k) {if(face_z_test(c,xlo,ylo,zlo,xhi,yhi)) continue;}
 588                                 else if(ek<k) {if(face_z_test(c,xlo,ylo,zhi,xhi,yhi)) continue;}
 589                                 else voro_fatal_error("Compute cell routine revisiting central block, which should never\nhappen.",VOROPP_INTERNAL_ERROR);
 590                         }
 591                 }
 592
 593                 // Now compute the region that we are going to test over, and
 594                 // set a displacement vector for the periodic cases
 595                 ijk=con.region_index(ci,cj,ck,ei,ej,ek,qx,qy,qz,disp);
 596
 597                 // Loop over all the elements in the block to test for cuts. It
 598                 // would be possible to exclude some of these cases by testing
 599                 // against mrs, but this will probably not save time.
 600                 if(co[ijk]>0) {
 601                         l=0;x2=x-qx;y2=y-qy;z2=z-qz;
 602                         do {
 603                                 x1=p[ijk][ps*l]-x2;
 604                                 y1=p[ijk][ps*l+1]-y2;
 605                                 z1=p[ijk][ps*l+2]-z2;
 606                                 rs=con.r_scale(x1*x1+y1*y1+z1*z1,ijk,l);
 607                                 if(!c.nplane(x1,y1,z1,rs,id[ijk][l])) return false;
 608                                 l++;
 609                         } while (l<co[ijk]);
 610                 }
 611
 612                 // If there's not much memory on the block list then add more
 613                 if((qu_s<=qu_e?(qu_l-qu_e)+(qu_s-qu):qu_s-qu_e)<18) add_list_memory(qu_s,qu_e);
 614
 615                 // Test the neighbors of the current block, and add them to the
 616                 // block list if they haven't already been tested
 617                 add_to_mask(ei,ej,ek,qu_e);
 618         }
 619
 620         return true;
 621 }
 622
 623 /** This function checks to see whether a particular block can possibly have
 624  * any intersection with a Voronoi cell, for the case when the closest point
 625  * from the cell center to the block is at a corner.
 626  * \param[in,out] c a reference to a Voronoi cell.
 627  * \param[in] (xl,yl,zl) the relative coordinates of the corner of the block
 628  *                       closest to the cell center.
 629  * \param[in] (xh,yh,zh) the relative coordinates of the corner of the block
 630  *                       furthest away from the cell center.
 631  * \return False if the block may intersect, true if does not. */
 632 template<class c_class>
 633 template<class v_cell>
 634 bool voro_compute<c_class>::corner_test(v_cell &c,double xl,double yl,double zl,double xh,double yh,double zh) {
 635         con.r_prime(xl*xl+yl*yl+zl*zl);
 636         if(c.plane_intersects_guess(xh,yl,zl,con.r_cutoff(xl*xh+yl*yl+zl*zl))) return false;
 637         if(c.plane_intersects(xh,yh,zl,con.r_cutoff(xl*xh+yl*yh+zl*zl))) return false;
 638         if(c.plane_intersects(xl,yh,zl,con.r_cutoff(xl*xl+yl*yh+zl*zl))) return false;
 639         if(c.plane_intersects(xl,yh,zh,con.r_cutoff(xl*xl+yl*yh+zl*zh))) return false;
 640         if(c.plane_intersects(xl,yl,zh,con.r_cutoff(xl*xl+yl*yl+zl*zh))) return false;
 641         if(c.plane_intersects(xh,yl,zh,con.r_cutoff(xl*xh+yl*yl+zl*zh))) return false;
 642         return true;
 643 }
 644
 645 /** This function checks to see whether a particular block can possibly have
 646  * any intersection with a Voronoi cell, for the case when the closest point
 647  * from the cell center to the block is on an edge which points along the x
 648  * direction.
 649  * \param[in,out] c a reference to a Voronoi cell.
 650  * \param[in] (x0,x1) the minimum and maximum relative x coordinates of the
 651  *                    block.
 652  * \param[in] (yl,zl) the relative y and z coordinates of the corner of the
 653  *                    block closest to the cell center.
 654  * \param[in] (yh,zh) the relative y and z coordinates of the corner of the
 655  *                    block furthest away from the cell center.
 656  * \return False if the block may intersect, true if does not. */
 657 template<class c_class>
 658 template<class v_cell>
 659 inline bool voro_compute<c_class>::edge_x_test(v_cell &c,double x0,double yl,double zl,double x1,double yh,double zh) {
 660         con.r_prime(yl*yl+zl*zl);
 661         if(c.plane_intersects_guess(x0,yl,zh,con.r_cutoff(yl*yl+zl*zh))) return false;
 662         if(c.plane_intersects(x1,yl,zh,con.r_cutoff(yl*yl+zl*zh))) return false;
 663         if(c.plane_intersects(x1,yl,zl,con.r_cutoff(yl*yl+zl*zl))) return false;
 664         if(c.plane_intersects(x0,yl,zl,con.r_cutoff(yl*yl+zl*zl))) return false;
 665         if(c.plane_intersects(x0,yh,zl,con.r_cutoff(yl*yh+zl*zl))) return false;
 666         if(c.plane_intersects(x1,yh,zl,con.r_cutoff(yl*yh+zl*zl))) return false;
 667         return true;
 668 }
 669
 670 /** This function checks to see whether a particular block can possibly have
 671  * any intersection with a Voronoi cell, for the case when the closest point
 672  * from the cell center to the block is on an edge which points along the y
 673  * direction.
 674  * \param[in,out] c a reference to a Voronoi cell.
 675  * \param[in] (y0,y1) the minimum and maximum relative y coordinates of the
 676  *                    block.
 677  * \param[in] (xl,zl) the relative x and z coordinates of the corner of the
 678  *                    block closest to the cell center.
 679  * \param[in] (xh,zh) the relative x and z coordinates of the corner of the
 680  *                    block furthest away from the cell center.
 681  * \return False if the block may intersect, true if does not. */
 682 template<class c_class>
 683 template<class v_cell>
 684 inline bool voro_compute<c_class>::edge_y_test(v_cell &c,double xl,double y0,double zl,double xh,double y1,double zh) {
 685         con.r_prime(xl*xl+zl*zl);
 686         if(c.plane_intersects_guess(xl,y0,zh,con.r_cutoff(xl*xl+zl*zh))) return false;
 687         if(c.plane_intersects(xl,y1,zh,con.r_cutoff(xl*xl+zl*zh))) return false;
 688         if(c.plane_intersects(xl,y1,zl,con.r_cutoff(xl*xl+zl*zl))) return false;
 689         if(c.plane_intersects(xl,y0,zl,con.r_cutoff(xl*xl+zl*zl))) return false;
 690         if(c.plane_intersects(xh,y0,zl,con.r_cutoff(xl*xh+zl*zl))) return false;
 691         if(c.plane_intersects(xh,y1,zl,con.r_cutoff(xl*xh+zl*zl))) return false;
 692         return true;
 693 }
 694
 695 /** This function checks to see whether a particular block can possibly have
 696  * any intersection with a Voronoi cell, for the case when the closest point
 697  * from the cell center to the block is on an edge which points along the z
 698  * direction.
 699  * \param[in,out] c a reference to a Voronoi cell.
 700  * \param[in] (z0,z1) the minimum and maximum relative z coordinates of the block.
 701  * \param[in] (xl,yl) the relative x and y coordinates of the corner of the
 702  *                    block closest to the cell center.
 703  * \param[in] (xh,yh) the relative x and y coordinates of the corner of the
 704  *                    block furthest away from the cell center.
 705  * \return False if the block may intersect, true if does not. */
 706 template<class c_class>
 707 template<class v_cell>
 708 inline bool voro_compute<c_class>::edge_z_test(v_cell &c,double xl,double yl,double z0,double xh,double yh,double z1) {
 709         con.r_prime(xl*xl+yl*yl);
 710         if(c.plane_intersects_guess(xl,yh,z0,con.r_cutoff(xl*xl+yl*yh))) return false;
 711         if(c.plane_intersects(xl,yh,z1,con.r_cutoff(xl*xl+yl*yh))) return false;
 712         if(c.plane_intersects(xl,yl,z1,con.r_cutoff(xl*xl+yl*yl))) return false;
 713         if(c.plane_intersects(xl,yl,z0,con.r_cutoff(xl*xl+yl*yl))) return false;
 714         if(c.plane_intersects(xh,yl,z0,con.r_cutoff(xl*xh+yl*yl))) return false;
 715         if(c.plane_intersects(xh,yl,z1,con.r_cutoff(xl*xh+yl*yl))) return false;
 716         return true;
 717 }
 718
 719 /** This function checks to see whether a particular block can possibly have
 720  * any intersection with a Voronoi cell, for the case when the closest point
 721  * from the cell center to the block is on a face aligned with the x direction.
 722  * \param[in,out] c a reference to a Voronoi cell.
 723  * \param[in] xl the minimum distance from the cell center to the face.
 724  * \param[in] (y0,y1) the minimum and maximum relative y coordinates of the
 725  *                    block.
 726  * \param[in] (z0,z1) the minimum and maximum relative z coordinates of the
 727  *                    block.
 728  * \return False if the block may intersect, true if does not. */
 729 template<class c_class>
 730 template<class v_cell>
 731 inline bool voro_compute<c_class>::face_x_test(v_cell &c,double xl,double y0,double z0,double y1,double z1) {
 732         con.r_prime(xl*xl);
 733         if(c.plane_intersects_guess(xl,y0,z0,con.r_cutoff(xl*xl))) return false;
 734         if(c.plane_intersects(xl,y0,z1,con.r_cutoff(xl*xl))) return false;
 735         if(c.plane_intersects(xl,y1,z1,con.r_cutoff(xl*xl))) return false;
 736         if(c.plane_intersects(xl,y1,z0,con.r_cutoff(xl*xl))) return false;
 737         return true;
 738 }
 739
 740 /** This function checks to see whether a particular block can possibly have
 741  * any intersection with a Voronoi cell, for the case when the closest point
 742  * from the cell center to the block is on a face aligned with the y direction.
 743  * \param[in,out] c a reference to a Voronoi cell.
 744  * \param[in] yl the minimum distance from the cell center to the face.
 745  * \param[in] (x0,x1) the minimum and maximum relative x coordinates of the
 746  *                    block.
 747  * \param[in] (z0,z1) the minimum and maximum relative z coordinates of the
 748  *                    block.
 749  * \return False if the block may intersect, true if does not. */
 750 template<class c_class>
 751 template<class v_cell>
 752 inline bool voro_compute<c_class>::face_y_test(v_cell &c,double x0,double yl,double z0,double x1,double z1) {
 753         con.r_prime(yl*yl);
 754         if(c.plane_intersects_guess(x0,yl,z0,con.r_cutoff(yl*yl))) return false;
 755         if(c.plane_intersects(x0,yl,z1,con.r_cutoff(yl*yl))) return false;
 756         if(c.plane_intersects(x1,yl,z1,con.r_cutoff(yl*yl))) return false;
 757         if(c.plane_intersects(x1,yl,z0,con.r_cutoff(yl*yl))) return false;
 758         return true;
 759 }
 760
 761 /** This function checks to see whether a particular block can possibly have
 762  * any intersection with a Voronoi cell, for the case when the closest point
 763  * from the cell center to the block is on a face aligned with the z direction.
 764  * \param[in,out] c a reference to a Voronoi cell.
 765  * \param[in] zl the minimum distance from the cell center to the face.
 766  * \param[in] (x0,x1) the minimum and maximum relative x coordinates of the
 767  *                    block.
 768  * \param[in] (y0,y1) the minimum and maximum relative y coordinates of the
 769  *                    block.
 770  * \return False if the block may intersect, true if does not. */
 771 template<class c_class>
 772 template<class v_cell>
 773 inline bool voro_compute<c_class>::face_z_test(v_cell &c,double x0,double y0,double zl,double x1,double y1) {
 774         con.r_prime(zl*zl);
 775         if(c.plane_intersects_guess(x0,y0,zl,con.r_cutoff(zl*zl))) return false;
 776         if(c.plane_intersects(x0,y1,zl,con.r_cutoff(zl*zl))) return false;
 777         if(c.plane_intersects(x1,y1,zl,con.r_cutoff(zl*zl))) return false;
 778         if(c.plane_intersects(x1,y0,zl,con.r_cutoff(zl*zl))) return false;
 779         return true;
 780 }
 781
 782
 783 /** This routine checks to see whether a point is within a particular distance
 784  * of a nearby region. If the point is within the distance of the region, then
 785  * the routine returns true, and computes the maximum distance from the point
 786  * to the region. Otherwise, the routine returns false.
 787  * \param[in] (di,dj,dk) the position of the nearby region to be tested,
 788  *                       relative to the region that the point is in.
 789  * \param[in] (fx,fy,fz) the displacement of the point within its region.
 790  * \param[in] (gxs,gys,gzs) the maximum squared distances from the point to the
 791  *                          sides of its region.
 792  * \param[out] crs a reference in which to return the maximum distance to the
 793  *                 region (only computed if the routine returns false).
 794  * \param[in] mrs the distance to be tested.
 795  * \return True if the region is further away than mrs, false if the region in
 796  *         within mrs. */
 797 template<class c_class>
 798 bool voro_compute<c_class>::compute_min_max_radius(int di,int dj,int dk,double fx,double fy,double fz,double gxs,double gys,double gzs,double &crs,double mrs) {
 799         double xlo,ylo,zlo;
 800         if(di>0) {
 801                 xlo=di*boxx-fx;
 802                 crs=xlo*xlo;
 803                 if(dj>0) {
 804                         ylo=dj*boxy-fy;
 805                         crs+=ylo*ylo;
 806                         if(dk>0) {
 807                                 zlo=dk*boxz-fz;
 808                                 crs+=zlo*zlo;if(con.r_ctest(crs,mrs)) return true;
 809                                 crs+=bxsq+2*(boxx*xlo+boxy*ylo+boxz*zlo);
 810                         } else if(dk<0) {
 811                                 zlo=(dk+1)*boxz-fz;
 812                                 crs+=zlo*zlo;if(con.r_ctest(crs,mrs)) return true;
 813                                 crs+=bxsq+2*(boxx*xlo+boxy*ylo-boxz*zlo);
 814                         } else {
 815                                 if(con.r_ctest(crs,mrs)) return true;
 816                                 crs+=boxx*(2*xlo+boxx)+boxy*(2*ylo+boxy)+gzs;
 817                         }
 818                 } else if(dj<0) {
 819                         ylo=(dj+1)*boxy-fy;
 820                         crs+=ylo*ylo;
 821                         if(dk>0) {
 822                                 zlo=dk*boxz-fz;
 823                                 crs+=zlo*zlo;if(con.r_ctest(crs,mrs)) return true;
 824                                 crs+=bxsq+2*(boxx*xlo-boxy*ylo+boxz*zlo);
 825                         } else if(dk<0) {
 826                                 zlo=(dk+1)*boxz-fz;
 827                                 crs+=zlo*zlo;if(con.r_ctest(crs,mrs)) return true;
 828                                 crs+=bxsq+2*(boxx*xlo-boxy*ylo-boxz*zlo);
 829                         } else {
 830                                 if(con.r_ctest(crs,mrs)) return true;
 831                                 crs+=boxx*(2*xlo+boxx)+boxy*(-2*ylo+boxy)+gzs;
 832                         }
 833                 } else {
 834                         if(dk>0) {
 835                                 zlo=dk*boxz-fz;
 836                                 crs+=zlo*zlo;if(con.r_ctest(crs,mrs)) return true;
 837                                 crs+=boxz*(2*zlo+boxz);
 838                         } else if(dk<0) {
 839                                 zlo=(dk+1)*boxz-fz;
 840                                 crs+=zlo*zlo;if(con.r_ctest(crs,mrs)) return true;
 841                                 crs+=boxz*(-2*zlo+boxz);
 842                         } else {
 843                                 if(con.r_ctest(crs,mrs)) return true;
 844                                 crs+=gzs;
 845                         }
 846                         crs+=gys+boxx*(2*xlo+boxx);
 847                 }
 848         } else if(di<0) {
 849                 xlo=(di+1)*boxx-fx;
 850                 crs=xlo*xlo;
 851                 if(dj>0) {
 852                         ylo=dj*boxy-fy;
 853                         crs+=ylo*ylo;
 854                         if(dk>0) {
 855                                 zlo=dk*boxz-fz;
 856                                 crs+=zlo*zlo;if(con.r_ctest(crs,mrs)) return true;
 857                                 crs+=bxsq+2*(-boxx*xlo+boxy*ylo+boxz*zlo);
 858                         } else if(dk<0) {
 859                                 zlo=(dk+1)*boxz-fz;
 860                                 crs+=zlo*zlo;if(con.r_ctest(crs,mrs)) return true;
 861                                 crs+=bxsq+2*(-boxx*xlo+boxy*ylo-boxz*zlo);
 862                         } else {
 863                                 if(con.r_ctest(crs,mrs)) return true;
 864                                 crs+=boxx*(-2*xlo+boxx)+boxy*(2*ylo+boxy)+gzs;
 865                         }
 866                 } else if(dj<0) {
 867                         ylo=(dj+1)*boxy-fy;
 868                         crs+=ylo*ylo;
 869                         if(dk>0) {
 870                                 zlo=dk*boxz-fz;
 871                                 crs+=zlo*zlo;if(con.r_ctest(crs,mrs)) return true;
 872                                 crs+=bxsq+2*(-boxx*xlo-boxy*ylo+boxz*zlo);
 873                         } else if(dk<0) {
 874                                 zlo=(dk+1)*boxz-fz;
 875                                 crs+=zlo*zlo;if(con.r_ctest(crs,mrs)) return true;
 876                                 crs+=bxsq+2*(-boxx*xlo-boxy*ylo-boxz*zlo);
 877                         } else {
 878                                 if(con.r_ctest(crs,mrs)) return true;
 879                                 crs+=boxx*(-2*xlo+boxx)+boxy*(-2*ylo+boxy)+gzs;
 880                         }
 881                 } else {
 882                         if(dk>0) {
 883                                 zlo=dk*boxz-fz;
 884                                 crs+=zlo*zlo;if(con.r_ctest(crs,mrs)) return true;
 885                                 crs+=boxz*(2*zlo+boxz);
 886                         } else if(dk<0) {
 887                                 zlo=(dk+1)*boxz-fz;
 888                                 crs+=zlo*zlo;if(con.r_ctest(crs,mrs)) return true;
 889                                 crs+=boxz*(-2*zlo+boxz);
 890                         } else {
 891                                 if(con.r_ctest(crs,mrs)) return true;
 892                                 crs+=gzs;
 893                         }
 894                         crs+=gys+boxx*(-2*xlo+boxx);
 895                 }
 896         } else {
 897                 if(dj>0) {
 898                         ylo=dj*boxy-fy;
 899                         crs=ylo*ylo;
 900                         if(dk>0) {
 901                                 zlo=dk*boxz-fz;
 902                                 crs+=zlo*zlo;if(con.r_ctest(crs,mrs)) return true;
 903                                 crs+=boxz*(2*zlo+boxz);
 904                         } else if(dk<0) {
 905                                 zlo=(dk+1)*boxz-fz;
 906                                 crs+=zlo*zlo;if(con.r_ctest(crs,mrs)) return true;
 907                                 crs+=boxz*(-2*zlo+boxz);
 908                         } else {
 909                                 if(con.r_ctest(crs,mrs)) return true;
 910                                 crs+=gzs;
 911                         }
 912                         crs+=boxy*(2*ylo+boxy);
 913                 } else if(dj<0) {
 914                         ylo=(dj+1)*boxy-fy;
 915                         crs=ylo*ylo;
 916                         if(dk>0) {
 917                                 zlo=dk*boxz-fz;
 918                                 crs+=zlo*zlo;if(con.r_ctest(crs,mrs)) return true;
 919                                 crs+=boxz*(2*zlo+boxz);
 920                         } else if(dk<0) {
 921                                 zlo=(dk+1)*boxz-fz;
 922                                 crs+=zlo*zlo;if(con.r_ctest(crs,mrs)) return true;
 923                                 crs+=boxz*(-2*zlo+boxz);
 924                         } else {
 925                                 if(con.r_ctest(crs,mrs)) return true;
 926                                 crs+=gzs;
 927                         }
 928                         crs+=boxy*(-2*ylo+boxy);
 929                 } else {
 930                         if(dk>0) {
 931                                 zlo=dk*boxz-fz;crs=zlo*zlo;if(con.r_ctest(crs,mrs)) return true;
 932                                 crs+=boxz*(2*zlo+boxz);
 933                         } else if(dk<0) {
 934                                 zlo=(dk+1)*boxz-fz;crs=zlo*zlo;if(con.r_ctest(crs,mrs)) return true;
 935                                 crs+=boxz*(-2*zlo+boxz);
 936                         } else {
 937                                 crs=0;
 938                                 voro_fatal_error("Min/max radius function called for central block, which should never\nhappen.",VOROPP_INTERNAL_ERROR);
 939                         }
 940                         crs+=gys;
 941                 }
 942                 crs+=gxs;
 943         }
 944         return false;
 945 }
 946
 947 template<class c_class>
 948 bool voro_compute<c_class>::compute_min_radius(int di,int dj,int dk,double fx,double fy,double fz,double mrs) {
 949         double t,crs;
 950
 951         if(di>0) {t=di*boxx-fx;crs=t*t;}
 952         else if(di<0) {t=(di+1)*boxx-fx;crs=t*t;}
 953         else crs=0;
 954
 955         if(dj>0) {t=dj*boxy-fy;crs+=t*t;}
 956         else if(dj<0) {t=(dj+1)*boxy-fy;crs+=t*t;}
 957
 958         if(dk>0) {t=dk*boxz-fz;crs+=t*t;}
 959         else if(dk<0) {t=(dk+1)*boxz-fz;crs+=t*t;}
 960
 961         return crs>con.r_max_add(mrs);
 962 }
 963
 964 /** Adds memory to the queue.
 965  * \param[in,out] qu_s a reference to the queue start pointer.
 966  * \param[in,out] qu_e a reference to the queue end pointer. */
 967 template<class c_class>
 968 inline void voro_compute<c_class>::add_list_memory(int*& qu_s,int*& qu_e) {
 969         qu_size<<=1;
 970         int *qu_n=new int[qu_size],*qu_c=qu_n;
 971 #if VOROPP_VERBOSE >=2
 972         fprintf(stderr,"List memory scaled up to %d\n",qu_size);
 973 #endif
 974         if(qu_s<=qu_e) {
 975                 while(qu_s<qu_e) *(qu_c++)=*(qu_s++);
 976         } else {
 977                 while(qu_s<qu_l) *(qu_c++)=*(qu_s++);qu_s=qu;
 978                 while(qu_s<qu_e) *(qu_c++)=*(qu_s++);
 979         }
 980         delete [] qu;
 981         qu_s=qu=qu_n;
 982         qu_l=qu+qu_size;
 983         qu_e=qu_c;
 984 }
 985
 986 // Explicit template instantiation
 987 template voro_compute<container>::voro_compute(container&,int,int,int);
 988 template voro_compute<container_poly>::voro_compute(container_poly&,int,int,int);
 989 template bool voro_compute<container>::compute_cell(voronoicell&,int,int,int,int,int);
 990 template bool voro_compute<container>::compute_cell(voronoicell_neighbor&,int,int,int,int,int);
 991 template void voro_compute<container>::find_voronoi_cell(double,double,double,int,int,int,int,particle_record&,double&);
 992 template bool voro_compute<container_poly>::compute_cell(voronoicell&,int,int,int,int,int);
 993 template bool voro_compute<container_poly>::compute_cell(voronoicell_neighbor&,int,int,int,int,int);
 994 template void voro_compute<container_poly>::find_voronoi_cell(double,double,double,int,int,int,int,particle_record&,double&);
 995
 996 // Explicit template instantiation
 997 template voro_compute<container_periodic>::voro_compute(container_periodic&,int,int,int);
 998 template voro_compute<container_periodic_poly>::voro_compute(container_periodic_poly&,int,int,int);
 999 template bool voro_compute<container_periodic>::compute_cell(voronoicell&,int,int,int,int,int);
1000 template bool voro_compute<container_periodic>::compute_cell(voronoicell_neighbor&,int,int,int,int,int);
1001 template void voro_compute<container_periodic>::find_voronoi_cell(double,double,double,int,int,int,int,particle_record&,double&);
1002 template bool voro_compute<container_periodic_poly>::compute_cell(voronoicell&,int,int,int,int,int);
1003 template bool voro_compute<container_periodic_poly>::compute_cell(voronoicell_neighbor&,int,int,int,int,int);
1004 template void voro_compute<container_periodic_poly>::find_voronoi_cell(double,double,double,int,int,int,int,particle_record&,double&);
1005
1006 }