branches/2d/src/c_loops_2d.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 c_loops_2d.cc
   8  * \brief Function implementations for the 2D loop classes. */
   9
  10 #include "c_loops_2d.hh"
  11
  12 namespace voro {
  13
  14 /** Initializes a c_loop_subset_2d object to scan over all particles within a
  15  * given circle.
  16  * \param[in] (vx,vy) the position vector of the center of the circle.
  17  * \param[in] r the radius of the circle.
  18  * \param[in] bounds_test whether to do detailed bounds checking. If this is
  19  *                        false then the class will loop over all particles in
  20  *                        blocks that overlap the given circle. If it is true,
  21  *                        the particle will only loop over the particles which
  22  *                        actually lie within the circle.
  23  * \return True if there is any valid point to loop over, false otherwise. */
  24 void c_loop_subset_2d::setup_circle(double vx,double vy,double r,bool bounds_test) {
  25         if(bounds_test) {mode=circle;v0=vx;v1=vy;v2=r*r;} else mode=no_check;
  26         ai=step_int((vx-ax-r)*xsp);
  27         bi=step_int((vx-ax+r)*xsp);
  28         aj=step_int((vy-ay-r)*ysp);
  29         bj=step_int((vy-ay+r)*ysp);
  30         setup_common();
  31 }
  32
  33 /** Initializes the class to loop over all particles in a rectangular subgrid
  34  * of blocks.
  35  * \param[in] (ai_,bi_) the subgrid range in the x-direction, inclusive of both
  36  *                      ends.
  37  * \param[in] (aj_,bj_) the subgrid range in the y-direction, inclusive of both
  38  *                      ends.
  39  * \return True if there is any valid point to loop over, false otherwise. */
  40 void c_loop_subset_2d::setup_intbox(int ai_,int bi_,int aj_,int bj_) {
  41         ai=ai_;bi=bi_;aj=aj_;bj=bj_;
  42         mode=no_check;
  43         setup_common();
  44 }
  45
  46 /** Sets up all of the common constants used for the loop.
  47  * \return True if there is any valid point to loop over, false otherwise. */
  48 void c_loop_subset_2d::setup_common() {
  49         if(!xperiodic) {
  50                 if(ai<0) {ai=0;if(bi<0) bi=0;}
  51                 if(bi>=nx) {bi=nx-1;if(ai>=nx) ai=nx-1;}
  52         }
  53         if(!yperiodic) {
  54                 if(aj<0) {aj=0;if(bj<0) bj=0;}
  55                 if(bj>=ny) {bj=ny-1;if(aj>=ny) aj=ny-1;}
  56         }
  57         ci=ai;cj=aj;
  58         di=i=step_mod(ci,nx);apx=px=step_div(ci,nx)*sx;
  59         dj=j=step_mod(cj,ny);apy=py=step_div(cj,ny)*sy;
  60         inc1=nx+di-step_mod(bi,nx);
  61         ij=di+nx*dj;
  62         q=0;
  63 }
  64
  65 /** Starts the loop by finding the first particle within the container to
  66  * consider.
  67  * \return True if there is any particle to consider, false otherwise. */
  68 bool c_loop_subset_2d::start() {
  69         while(co[ij]==0) {if(!next_block()) return false;}
  70         while(mode!=no_check&&out_of_bounds()) {
  71                 q++;
  72                 while(q>=co[ij]) {q=0;if(!next_block()) return false;}
  73         }
  74         return true;
  75 }
  76
  77 /** Initializes the class to loop over all particles in a rectangular box.
  78  * \param[in] (xmin,xmax) the minimum and maximum x coordinates of the box.
  79  * \param[in] (ymin,ymax) the minimum and maximum y coordinates of the box.
  80  * \param[in] bounds_test whether to do detailed bounds checking. If this is
  81  *                        false then the class will loop over all particles in
  82  *                        blocks that overlap the given box. If it is true, the
  83  *                        particle will only loop over the particles which
  84  *                        actually lie within the box.
  85  * \return True if there is any valid point to loop over, false otherwise. */
  86 void c_loop_subset_2d::setup_box(double xmin,double xmax,double ymin,double ymax,bool bounds_test) {
  87         if(bounds_test) {mode=rectangle;v0=xmin;v1=xmax;v2=ymin;v3=ymax;} else mode=no_check;
  88         ai=step_int((xmin-ax)*xsp);
  89         bi=step_int((xmax-ax)*xsp);
  90         aj=step_int((ymin-ay)*ysp);
  91         bj=step_int((ymax-ay)*ysp);
  92         setup_common();
  93 }
  94
  95 /** Computes whether the current point is out of bounds, relative to the
  96  * current loop setup.
  97  * \return True if the point is out of bounds, false otherwise. */
  98 bool c_loop_subset_2d::out_of_bounds() {
  99         double *pp(p[ij]+ps*q);
 100         if(mode==circle) {
 101                 double fx(*pp+px-v0),fy(pp[1]+py-v1);
 102                 return fx*fx+fy*fy>v2;
 103         } else {
 104                 double f(*pp+px);if(f<v0||f>v1) return true;
 105                 f=pp[1]+py;return f<v2||f>v3;
 106         }
 107 }
 108
 109 /** Returns the next block to be tested in a loop, and updates the periodicity
 110  * vector if necessary. */
 111 bool c_loop_subset_2d::next_block() {
 112         if(i<bi) {
 113                 i++;
 114                 if(ci<nx-1) {ci++;ij++;} else {ci=0;ij+=1-nx;px+=sx;}
 115                 return true;
 116         } else if(j<bj) {
 117                 i=ai;ci=di;px=apx;j++;
 118                 if(cj<ny-1) {cj++;ij+=inc1;} else {cj=0;ij+=inc1-nxy;py+=sy;}
 119                 return true;
 120         } else return false;
 121 }
 122
 123 /** Extends the memory available for storing the ordering. */
 124 void particle_order::add_ordering_memory() {
 125         int *no=new int[size<<2],*nop=no,*opp=o;
 126         while(opp<op) *(nop++)=*(opp++);
 127         delete [] o;
 128         size<<=1;o=no;op=nop;
 129 }
 130
 131 }