trunk/src/v_base.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_base.cc
   8  * \brief Function implementations for the base Voronoi container class. */
   9
  10 #include "v_base.hh"
  11 #include "config.hh"
  12
  13 namespace voro {
  14
  15 /** This function is called during container construction. The routine scans
  16  * all of the worklists in the wl[] array. For a given worklist of blocks
  17  * labeled \f$w_1\f$ to \f$w_n\f$, it computes a sequence \f$r_0\f$ to
  18  * \f$r_n\f$ so that $r_i$ is the minimum distance to all the blocks
  19  * \f$w_{j}\f$ where \f$j>i\f$ and all blocks outside the worklist. The values
  20  * of \f$r_n\f$ is calculated first, as the minimum distance to any block in
  21  * the shell surrounding the worklist. The \f$r_i\f$ are then computed in
  22  * reverse order by considering the distance to \f$w_{i+1}\f$. */
  23 voro_base::voro_base(int nx_,int ny_,int nz_,double boxx_,double boxy_,double boxz_) :
  24         nx(nx_), ny(ny_), nz(nz_), nxy(nx_*ny_), nxyz(nxy*nz_), boxx(boxx_), boxy(boxy_), boxz(boxz_),
  25         xsp(1/boxx_), ysp(1/boxy_), zsp(1/boxz_), mrad(new double[wl_hgridcu*wl_seq_length]) {
  26         const unsigned int b1=1<<21,b2=1<<22,b3=1<<24,b4=1<<25,b5=1<<27,b6=1<<28;
  27         const double xstep=boxx/wl_fgrid,ystep=boxy/wl_fgrid,zstep=boxz/wl_fgrid;
  28         int i,j,k,lx,ly,lz,q;
  29         unsigned int f,*e=const_cast<unsigned int*> (wl);
  30         double xlo,ylo,zlo,xhi,yhi,zhi,minr,*radp=mrad;
  31         for(zlo=0,zhi=zstep,lz=0;lz<wl_hgrid;zlo=zhi,zhi+=zstep,lz++) {
  32                 for(ylo=0,yhi=ystep,ly=0;ly<wl_hgrid;ylo=yhi,yhi+=ystep,ly++) {
  33                         for(xlo=0,xhi=xstep,lx=0;lx<wl_hgrid;xlo=xhi,xhi+=xstep,lx++) {
  34                                 minr=large_number;
  35                                 for(q=e[0]+1;q<wl_seq_length;q++) {
  36                                         f=e[q];
  37                                         i=(f&127)-64;
  38                                         j=(f>>7&127)-64;
  39                                         k=(f>>14&127)-64;
  40                                         if((f&b2)==b2) {
  41                                                 compute_minimum(minr,xlo,xhi,ylo,yhi,zlo,zhi,i-1,j,k);
  42                                                 if((f&b1)==0) compute_minimum(minr,xlo,xhi,ylo,yhi,zlo,zhi,i+1,j,k);
  43                                         } else if((f&b1)==b1) compute_minimum(minr,xlo,xhi,ylo,yhi,zlo,zhi,i+1,j,k);
  44                                         if((f&b4)==b4) {
  45                                                 compute_minimum(minr,xlo,xhi,ylo,yhi,zlo,zhi,i,j-1,k);
  46                                                 if((f&b3)==0) compute_minimum(minr,xlo,xhi,ylo,yhi,zlo,zhi,i,j+1,k);
  47                                         } else if((f&b3)==b3) compute_minimum(minr,xlo,xhi,ylo,yhi,zlo,zhi,i,j+1,k);
  48                                         if((f&b6)==b6) {
  49                                                 compute_minimum(minr,xlo,xhi,ylo,yhi,zlo,zhi,i,j,k-1);
  50                                                 if((f&b5)==0) compute_minimum(minr,xlo,xhi,ylo,yhi,zlo,zhi,i,j,k+1);
  51                                         } else if((f&b5)==b5) compute_minimum(minr,xlo,xhi,ylo,yhi,zlo,zhi,i,j,k+1);
  52                                 }
  53                                 q--;
  54                                 while(q>0) {
  55                                         radp[q]=minr;
  56                                         f=e[q];
  57                                         i=(f&127)-64;
  58                                         j=(f>>7&127)-64;
  59                                         k=(f>>14&127)-64;
  60                                         compute_minimum(minr,xlo,xhi,ylo,yhi,zlo,zhi,i,j,k);
  61                                         q--;
  62                                 }
  63                                 *radp=minr;
  64                                 e+=wl_seq_length;
  65                                 radp+=wl_seq_length;
  66                         }
  67                 }
  68         }
  69 }
  70
  71 /** Computes the minimum distance from a subregion to a given block. If this distance
  72  * is smaller than the value of minr, then it passes
  73  * \param[in,out] minr a pointer to the current minimum distance. If the distance
  74  *                     computed in this function is smaller, then this distance is
  75  *                     set to the new one.
  76  * \param[out] (xlo,ylo,zlo) the lower coordinates of the subregion being
  77  *                           considered.
  78  * \param[out] (xhi,yhi,zhi) the upper coordinates of the subregion being
  79  *                           considered.
  80  * \param[in] (ti,tj,tk) the coordinates of the block. */
  81 void voro_base::compute_minimum(double &minr,double &xlo,double &xhi,double &ylo,double &yhi,double &zlo,double &zhi,int ti,int tj,int tk) {
  82         double radsq,temp;
  83         if(ti>0) {temp=boxx*ti-xhi;radsq=temp*temp;}
  84         else if(ti<0) {temp=xlo-boxx*(1+ti);radsq=temp*temp;}
  85         else radsq=0;
  86
  87         if(tj>0) {temp=boxy*tj-yhi;radsq+=temp*temp;}
  88         else if(tj<0) {temp=ylo-boxy*(1+tj);radsq+=temp*temp;}
  89
  90         if(tk>0) {temp=boxz*tk-zhi;radsq+=temp*temp;}
  91         else if(tk<0) {temp=zlo-boxz*(1+tk);radsq+=temp*temp;}
  92
  93         if(radsq<minr) minr=radsq;
  94 }
  95
  96 /** Checks to see whether "%n" appears in a format sequence to determine
  97  * whether neighbor information is required or not.
  98  * \param[in] format the format string to check.
  99  * \return True if a "%n" is found, false otherwise. */
 100 bool voro_base::contains_neighbor(const char *format) {
 101         char *fmp=(const_cast<char*>(format));
 102
 103         // Check to see if "%n" appears in the format sequence
 104         while(*fmp!=0) {
 105                 if(*fmp=='%') {
 106                         fmp++;
 107                         if(*fmp=='n') return true;
 108                         else if(*fmp==0) return false;
 109                 }
 110                 fmp++;
 111         }
 112
 113         return false;
 114 }
 115
 116 #include "v_base_wl.cc"
 117
 118 }