branches/2d/src/v_connect/v_base_2d.cc

   1 // Voro++, a 2D and 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_2d.cc
   8  * \brief Function implementations for the base 2D Voronoi container class. */
   9 //#include <stdio.h>
  10
  11 using namespace std;
  12 #include "v_base_2d.hh"
  13 #include "config.hh"
  14
  15 namespace voro {
  16
  17 /** This function is called during container construction. The routine scans
  18  * all of the worklists in the wl[] array. For a given worklist of blocks
  19  * labeled \f$w_1\f$ to \f$w_n\f$, it computes a sequence \f$r_0\f$ to
  20  * \f$r_n\f$ so that $r_i$ is the minimum distance to all the blocks
  21  * \f$w_{j}\f$ where \f$j>i\f$ and all blocks outside the worklist. The values
  22  * of \f$r_n\f$ is calculated first, as the minimum distance to any block in
  23  * the shell surrounding the worklist. The \f$r_i\f$ are then computed in
  24  * reverse order by considering the distance to \f$w_{i+1}\f$. */
  25 voro_base_2d::voro_base_2d(int nx_,int ny_,double boxx_,double boxy_) :
  26         nx(nx_), ny(ny_), nxy(nx_*ny_), boxx(boxx_), boxy(boxy_),
  27         xsp(1/boxx_), ysp(1/boxy_), mrad(new double[wl_hgridsq_2d*wl_seq_length_2d]) {
  28         full_connect=false;
  29         const unsigned int b1=1<<21,b2=1<<22,b3=1<<24,b4=1<<25;
  30         const double xstep=boxx/wl_fgrid_2d,ystep=boxy/wl_fgrid_2d;
  31         int i,j,lx,ly,q;
  32         unsigned int f,*e=const_cast<unsigned int*> (wl);
  33         double xlo,ylo,xhi,yhi,minr,*radp=mrad;
  34         for(ylo=0,yhi=ystep,ly=0;ly<wl_hgrid_2d;ylo=yhi,yhi+=ystep,ly++) {
  35                 for(xlo=0,xhi=xstep,lx=0;lx<wl_hgrid_2d;xlo=xhi,xhi+=xstep,lx++) {
  36                         minr=large_number;
  37                         for(q=e[0]+1;q<wl_seq_length_2d;q++) {
  38                                 f=e[q];
  39                                 i=(f&127)-64;
  40                                 j=(f>>7&127)-64;
  41                                 if((f&b2)==b2) {
  42                                         compute_minimum(minr,xlo,xhi,ylo,yhi,i-1,j);
  43                                         if((f&b1)==0) compute_minimum(minr,xlo,xhi,ylo,yhi,i+1,j);
  44                                 } else if((f&b1)==b1) compute_minimum(minr,xlo,xhi,ylo,yhi,i+1,j);
  45                                 if((f&b4)==b4) {
  46                                         compute_minimum(minr,xlo,xhi,ylo,yhi,i,j-1);
  47                                         if((f&b3)==0) compute_minimum(minr,xlo,xhi,ylo,yhi,i,j+1);
  48                                 } else if((f&b3)==b3) compute_minimum(minr,xlo,xhi,ylo,yhi,i,j+1);
  49                         }
  50                         q--;
  51                         while(q>0) {
  52                                 radp[q]=minr;
  53                                 f=e[q];
  54                                 i=(f&127)-64;
  55                                 j=(f>>7&127)-64;
  56                                 compute_minimum(minr,xlo,xhi,ylo,yhi,i,j);
  57                                 q--;
  58                         }
  59                         *radp=minr;
  60                         e+=wl_seq_length_2d;
  61                         radp+=wl_seq_length_2d;
  62                 }
  63         }
  64 }
  65
  66 /** Computes the minimum distance from a subregion to a given block. If this distance
  67  * is smaller than the value of minr, then it passes
  68  * \param[in,out] minr a pointer to the current minimum distance. If the distance
  69  *                     computed in this function is smaller, then this distance is
  70  *                     set to the new one.
  71  * \param[out] (xlo,ylo) the lower coordinates of the subregion being
  72  *                       considered.
  73  * \param[out] (xhi,yhi) the upper coordinates of the subregion being
  74  *                       considered.
  75  * \param[in] (ti,tj) the coordinates of the block. */
  76 void voro_base_2d::compute_minimum(double &minr,double &xlo,double &xhi,double &ylo,double &yhi,int ti,int tj) {
  77         double radsq,temp;
  78         if(ti>0) {temp=boxx*ti-xhi;radsq=temp*temp;}
  79         else if(ti<0) {temp=xlo-boxx*(1+ti);radsq=temp*temp;}
  80         else radsq=0;
  81
  82         if(tj>0) {temp=boxy*tj-yhi;radsq+=temp*temp;}
  83         else if(tj<0) {temp=ylo-boxy*(1+tj);radsq+=temp*temp;}
  84
  85         if(radsq<minr) minr=radsq;
  86 }
  87
  88 /** Checks to see whether "%n" appears in a format sequence to determine
  89  * whether neighbor information is required or not.
  90  * \param[in] format the format string to check.
  91  * \return True if a "%n" is found, false otherwise. */
  92 bool voro_base_2d::contains_neighbor(const char *format) {
  93         char *fmp=(const_cast<char*>(format));
  94
  95         // Check to see if "%n" appears in the format sequence
  96         while(*fmp!=0) {
  97                 if(*fmp=='%') {
  98                         fmp++;
  99                         if(*fmp=='n') return true;
 100                         else if(*fmp==0) return false;
 101                 }
 102                 fmp++;
 103         }
 104
 105         return false;
 106 }
 107
 108 /*bool voro_base_2d::contains_neighbor_global(const char *format) {
 109                 char *fmp=(const_cast<char*>(format));
 110
 111                         // Check to see if "%N" appears in the format sequence
 112                         while(*fmp!=0) {
 113                                 if(*fmp=='%') {
 114                                         fmp++;
 115                                         if(*fmp=='N') return true;
 116                                         else if(*fmp==0) return false;
 117                                 }
 118                                 fmp++;
 119                         }
 120
 121                         return false;
 122 }
 123 void voro_base_2d::add_globne_info(int pid, int *nel, int length){
 124         for(int i=0;i<length;i++){
 125                 if(nel[i]>=0){
 126                         globne[((totpar*pid)+nel[i])/32] |= (unsigned int)(1 << (((totpar*pid)+nel[i])%32));
 127                 }
 128         }
 129
 130 }
 131 void voro_base_2d::print_globne(FILE *fp){
 132         int j=0;
 133         fprintf(fp, "Global neighbor info: Format \n [Particle-ID] \t [Neighbors] \n [Particle-ID] \t [Neighbors]");
 134         for(int i=0; i<totpar; i++){
 135                 fprintf(fp,"\n %d \t",i);
 136                 for(j=0;j<totpar;j++){
 137                         if((globne[(((i*totpar)+j)/32)] & (unsigned int)(1 << (((i*totpar)+j)%32))) != 0){
 138                                 fprintf(fp,"%d \t",j);
 139                         }
 140                 }
 141         }
 142 }
 143 //returns true if a and be are connected(share an edge,) note that add_globne_info must be called with all cells first.
 144 bool voro_base_2d::connected(int a, int b){
 145         return (((globne[(((a*totpar)+b)/32)]) & (unsigned int)(1 << (((a*totpar)+b)%32))) != 0);
 146 }*/
 147 #include "v_base_wl_2d.cc"
 148 }