branches/2d/src/container_2d_old.cc

   1 /** \file container_2d.cc
   2  * \brief Function implementations for the container_2d class. */
   3
   4 #include "container_2d.hh"
   5
   6 /** The class constructor sets up the geometry of container, initializing the
   7  * minimum and maximum coordinates in each direction, and setting whether each
   8  * direction is periodic or not. It divides the container into a rectangular
   9  * grid of blocks, and allocates memory for each of these for storing particle
  10  * positions and IDs.
  11  * \param[in] (ax_,bx_) the minimum and maximum x coordinates.
  12  * \param[in] (ay_,by_) the minimum and maximum y coordinates.
  13  * \param[in] (nx_,ny_) the number of grid blocks in each of the three
  14  *                       coordinate directions.
  15  * \param[in] (xperiodic_,yperiodic_) flags setting whether the container is periodic
  16  *                        in each coordinate direction.
  17  * \param[in] init_mem the initial memory allocation for each block. */
  18 container_2d::container_2d(double ax_,double bx_,double ay_,
  19                 double by_,int nx_,int ny_,bool xperiodic_,bool yperiodic_,int init_mem)
  20         : ax(ax_), bx(bx_), ay(ay_), by(by_), boxx((bx_-ax_)/nx_), boxy((by_-ay_)/ny_),
  21         xsp(1/boxx), ysp(1/boxy), nx(nx_), ny(ny_), nxy(nx*ny),
  22         xperiodic(xperiodic_), yperiodic(yperiodic_),
  23         co(new int[nxy]), mem(new int[nxy]), id(new int*[nxy]), p(new double*[nxy]) {
  24         int l;
  25         for(l=0;l<nxy;l++) co[l]=0;
  26         for(l=0;l<nxy;l++) mem[l]=init_mem;
  27         for(l=0;l<nxy;l++) id[l]=new int[init_mem];
  28         for(l=0;l<nxy;l++) p[l]=new double[2*init_mem];
  29 }
  30
  31 /** The container destructor frees the dynamically allocated memory. */
  32 container_2d::~container_2d() {
  33         int l;
  34         for(l=nxy-1;l>=0;l--) delete [] p[l];
  35         for(l=nxy-1;l>=0;l--) delete [] id[l];
  36         delete [] p;
  37         delete [] id;
  38         delete [] mem;
  39         delete [] co;
  40 }
  41
  42 /** Put a particle into the correct region of the container.
  43  * \param[in] n the numerical ID of the inserted particle.
  44  * \param[in] (x,y) the position vector of the inserted particle. */
  45 void container_2d::put(int n,double x,double y) {
  46         int ij;
  47         if(put_locate_block(ij,x,y)) {
  48                 id[ij][co[ij]]=n;
  49                 double *pp(p[ij]+2*co[ij]++);
  50                 *(pp++)=x;*pp=y;
  51         }
  52 }
  53
  54 /** This routine takes a particle position vector, tries to remap it into the
  55  * primary domain. If successful, it computes the region into which it can be
  56  * stored and checks that there is enough memory within this region to store
  57  * it.
  58  * \param[out] ij the region index.
  59  * \param[in,out] (x,y) the particle position, remapped into the primary
  60  *                      domain if necessary.
  61  * \return True if the particle can be successfully placed into the container,
  62  * false otherwise. */
  63 inline bool container_2d::put_locate_block(int &ij,double &x,double &y) {
  64         if(put_remap(ij,x,y)) {
  65                 if(co[ij]==mem[ij]) add_particle_memory(ij);
  66                 return true;
  67         }
  68 #if VOROPP_REPORT_OUT_OF_BOUNDS ==1
  69         fprintf(stderr,"Out of bounds: (x,y)=(%g,%g)\n",x,y);
  70 #endif
  71         return false;
  72 }
  73
  74 /** Takes a particle position vector and computes the region index into which
  75  * it should be stored. If the container is periodic, then the routine also
  76  * maps the particle position to ensure it is in the primary domain. If the
  77  * container is not periodic, the routine bails out.
  78  * \param[out] ij the region index.
  79  * \param[in,out] (x,y) the particle position, remapped into the primary
  80  *                      domain if necessary.
  81  * \return True if the particle can be successfully placed into the container,
  82  * false otherwise. */
  83 inline bool container_2d::put_remap(int &ij,double &x,double &y) {
  84         int l;
  85
  86         ij=step_int((x-ax)*xsp);
  87         if(xperiodic) {l=step_mod(ij,nx);x+=boxx*(l-ij);ij=l;}
  88         else if(ij<0||ij>=nx) return false;
  89
  90         int j(step_int((y-ay)*ysp));
  91         if(yperiodic) {l=step_mod(j,ny);y+=boxy*(l-j);j=l;}
  92         else if(j<0||j>=ny) return false;
  93
  94         ij+=nx*j;
  95         return true;
  96 }
  97
  98 /** Increase memory for a particular region.
  99  * \param[in] i the index of the region to reallocate. */
 100 void container_2d::add_particle_memory(int i) {
 101         int l,*idp;double *pp;
 102         mem[i]<<=1;
 103         if(mem[i]>max_particle_memory)
 104                 voropp_fatal_error("Absolute maximum particle memory allocation exceeded",VOROPP_MEMORY_ERROR);
 105 #if VOROPP_VERBOSE >=3
 106         fprintf(stderr,"Particle memory in region %d scaled up to %d\n",i,mem[i]);
 107 #endif
 108         idp=new int[mem[i]];
 109         for(l=0;l<co[i];l++) idp[l]=id[i][l];
 110         pp=new double[2*mem[i]];
 111         for(l=0;l<2*co[i];l++) pp[l]=p[i][l];
 112         delete [] id[i];id[i]=idp;
 113         delete [] p[i];p[i]=pp;
 114 }
 115
 116 /** Imports a list of particles from an input stream.
 117  * \param[in] fp a file handle to read from. */
 118 void container_2d::import(FILE *fp) {
 119         int i,j;
 120         double x,y;
 121         while((j=fscanf(fp,"%d %lg %lg",&i,&x,&y))==3) put(i,x,y);
 122         if(j!=EOF) voropp_fatal_error("File import error",VOROPP_FILE_ERROR);
 123 }
 124
 125 /** Clears a container of particles. */
 126 void container_2d::clear() {
 127         for(int* cop=co;cop<co+nxy;cop++) *cop=0;
 128 }
 129
 130 /** Dumps all the particle positions and IDs to a file.
 131  * \param[in] fp the file handle to write to. */
 132 void container_2d::draw_particles(FILE *fp) {
 133         int ij,q;
 134         for(ij=0;ij<nxy;ij++) for(q=0;q<co[ij];q++)
 135                 fprintf(fp,"%d %g %g\n",id[ij][q],p[ij][2*q],p[ij][2*q+1]);
 136 }
 137
 138 /** Dumps all the particle positions in POV-Ray format.
 139  * \param[in] fp the file handle to write to. */
 140 void container_2d::draw_particles_pov(FILE *fp) {
 141         int ij,q;
 142         for(ij=0;ij<nxy;ij++) for(q=0;q<co[ij];q++)
 143                 fprintf(fp,"// id %d\nsphere{<%g,%g,0>,s\n",id[ij][q],p[ij][2*q],p[ij][2*q+1]);
 144 }
 145
 146 /** Computes the Voronoi cells for all particles and saves the output in
 147  * gnuplot format.
 148  * \param[in] fp a file handle to write to. */
 149 void container_2d::draw_cells_gnuplot(FILE *fp) {
 150         int i,j,ij=0,q;
 151         double x,y;
 152         voronoicell_2d c;
 153         for(j=0;j<ny;j++) for(i=0;i<nx;i++,ij++) for(q=0;q<co[ij];q++) {
 154                 x=p[ij][2*q];y=p[ij][2*q+1];
 155                 if(compute_cell_sphere(c,i,j,ij,q,x,y)) c.draw_gnuplot(x,y,fp);
 156         }
 157 }
 158
 159 /** Computes the Voronoi cells for all particles within a rectangular box, and
 160  * saves the output in POV-Ray format.
 161  * \param[in] fp a file handle to write to. */
 162 void container_2d::draw_cells_pov(FILE *fp) {
 163         int i,j,ij=0,q;
 164         double x,y;
 165         voronoicell_2d c;
 166         for(j=0;j<ny;j++) for(i=0;i<nx;i++,ij++) for(q=0;q<co[ij];q++) {
 167                 x=p[ij][2*q];y=p[ij][2*q+1];
 168                 if(compute_cell_sphere(c,i,j,ij,q,x,y)) {
 169                         fprintf(fp,"// cell %d\n",id[ij][q]);
 170                         c.draw_pov(x,y,0,fp);
 171                 }
 172         }
 173 }
 174
 175 /** Computes the Voronoi cells for all particles in the container, and for each
 176  * cell, outputs a line containing custom information about the cell structure.
 177  * The output format is specified using an input string with control sequences
 178  * similar to the standard C printf() routine.
 179  * \param[in] format the format of the output lines, using control sequences to
 180  *                   denote the different cell statistics.
 181  * \param[in] fp a file handle to write to. */
 182 void container_2d::print_custom(const char *format,FILE *fp) {
 183         int i,j,ij=0,q;
 184         double x,y;
 185         voronoicell_2d c;
 186         for(j=0;j<ny;j++) for(i=0;i<nx;i++,ij++) for(q=0;q<co[ij];q++) {
 187                 x=p[ij][2*q];y=p[ij][2*q+1];
 188                 if(compute_cell_sphere(c,i,j,ij,q,x,y)) c.output_custom(format,id[ij][q],x,y,default_radius,fp);
 189         }
 190 }
 191
 192 /** Initializes a voronoicell_2d class to fill the entire container.
 193  * \param[in] c a reference to a voronoicell_2d class.
 194  * \param[in] (x,y) the position of the particle that . */
 195 bool container_2d::initialize_voronoicell(voronoicell_2d &c,double x,double y) {
 196         double x1,x2,y1,y2;
 197         if(xperiodic) x1=-(x2=0.5*(bx-ax));else {x1=ax-x;x2=bx-x;}
 198         if(yperiodic) y1=-(y2=0.5*(by-ay));else {y1=ay-y;y2=by-y;}
 199         c.init(x1,x2,y1,y2);
 200         return true;
 201 }
 202
 203 /** Computes all Voronoi cells and sums their areas.
 204  * \return The computed area. */
 205 double container_2d::sum_cell_areas() {
 206         int i,j,ij=0,q;
 207         double x,y,sum=0;
 208         voronoicell_2d c;
 209         for(j=0;j<ny;j++) for(i=0;i<nx;i++,ij++) for(q=0;q<co[ij];q++) {
 210                 x=p[ij][2*q];y=p[ij][2*q+1];
 211                 if(compute_cell_sphere(c,i,j,ij,q,x,y)) sum+=c.area();
 212         }
 213         return sum;
 214 }
 215
 216 /** Computes all of the Voronoi cells in the container, but does nothing
 217  * with the output. It is useful for measuring the pure computation time
 218  * of the Voronoi algorithm, without any additional calculations such as
 219  * volume evaluation or cell output. */
 220 void container_2d::compute_all_cells() {
 221         int i,j,ij=0,q;
 222         voronoicell_2d c;
 223         for(j=0;j<ny;j++) for(i=0;i<nx;i++,ij++) for(q=0;q<co[ij];q++)
 224                 compute_cell_sphere(c,i,j,ij,q);
 225 }
 226
 227 /** This routine computes the Voronoi cell for a give particle, by successively
 228  * testing over particles within larger and larger concentric circles. This
 229  * routine is simple and fast, although it may not work well for anisotropic
 230  * arrangements of particles.
 231  * \param[in,out] c a reference to a voronoicell object.
 232  * \param[in] (i,j) the coordinates of the block that the test particle is
 233  *                  in.
 234  * \param[in] ij the index of the block that the test particle is in, set to
 235  *               i+nx*j.
 236  * \param[in] s the index of the particle within the test block.
 237  * \param[in] (x,y) the coordinates of the particle.
 238  * \return False if the Voronoi cell was completely removed during the
 239  * computation and has zero volume, true otherwise. */
 240 bool container_2d::compute_cell_sphere(voronoicell_2d &c,int i,int j,int ij,int s,double x,double y) {
 241
 242         // This length scale determines how large the spherical shells should
 243         // be, and it should be set to approximately the particle diameter
 244         const double length_scale=0.5*sqrt((bx-ax)*(by-ay)/(nx*ny));
 245
 246         double x1,y1,qx,qy,lr=0,lrs=0,ur,urs,rs;
 247         int q,t;
 248         voropp_loop_2d l(*this);
 249
 250         if(!initialize_voronoicell(c,x,y)) return false;
 251
 252         // Now the cell is cut by testing neighboring particles in concentric
 253         // shells. Once the test shell becomes twice as large as the Voronoi
 254         // cell we can stop testing.
 255         while(lrs<c.max_radius_squared()) {
 256                 ur=lr+0.5*length_scale;urs=ur*ur;
 257                 t=l.init(x,y,ur,qx,qy);
 258                 do {
 259                         for(q=0;q<co[t];q++) {
 260                                 x1=p[t][2*q]+qx-x;y1=p[t][2*q+1]+qy-y;
 261                                 rs=x1*x1+y1*y1;
 262                                 if(lrs-tolerance<rs&&rs<urs&&(q!=s||ij!=t)) {
 263                                         if(!c.plane(x1,y1,rs)) return false;
 264                                 }
 265                         }
 266                 } while((t=l.inc(qx,qy))!=-1);
 267                 lr=ur;lrs=urs;
 268         }
 269         return true;
 270 }
 271
 272 /** Creates a voropp_loop_2d object, by setting the necessary constants about the
 273  * container geometry from a pointer to the current container class.
 274  * \param[in] con a reference to the associated container class. */
 275 voropp_loop_2d::voropp_loop_2d(container_2d &con) : boxx(con.bx-con.ax), boxy(con.by-con.ay),
 276         xsp(con.xsp),ysp(con.ysp),
 277         ax(con.ax),ay(con.ay),nx(con.nx),ny(con.ny),nxy(con.nxy),
 278         xperiodic(con.xperiodic),yperiodic(con.yperiodic) {}
 279
 280 /** Initializes a voropp_loop_2d object, by finding all blocks which are within a
 281  * given sphere. It calculates the index of the first block that needs to be
 282  * tested and sets the periodic displacement vector accordingly.
 283  * \param[in] (vx,vy) the position vector of the center of the sphere.
 284  * \param[in] r the radius of the sphere.
 285  * \param[out] (px,py) the periodic displacement vector for the first block to
 286  *                     be tested.
 287  * \return The index of the first block to be tested. */
 288 int voropp_loop_2d::init(double vx,double vy,double r,double &px,double &py) {
 289         ai=step_int((vx-ax-r)*xsp);
 290         bi=step_int((vx-ax+r)*xsp);
 291         if(!xperiodic) {
 292                 if(ai<0) {ai=0;if(bi<0) bi=0;}
 293                 if(bi>=nx) {bi=nx-1;if(ai>=nx) ai=nx-1;}
 294         }
 295         aj=step_int((vy-ay-r)*ysp);
 296         bj=step_int((vy-ay+r)*ysp);
 297         if(!yperiodic) {
 298                 if(aj<0) {aj=0;if(bj<0) bj=0;}
 299                 if(bj>=ny) {bj=ny-1;if(aj>=ny) aj=ny-1;}
 300         }
 301         i=ai;j=aj;
 302         aip=ip=step_mod(i,nx);apx=px=step_div(i,nx)*boxx;
 303         ajp=jp=step_mod(j,ny);apy=py=step_div(j,ny)*boxy;
 304         inc1=aip-step_mod(bi,nx)+nx;
 305         s=aip+nx*ajp;
 306         return s;
 307 }
 308
 309 /** Initializes a voropp_loop_2d object, by finding all blocks which overlap a given
 310  * rectangular box. It calculates the index of the first block that needs to be
 311  * tested and sets the periodic displacement vector (px,py,pz) accordingly.
 312  * \param[in] (xmin,xmax) the minimum and maximum x coordinates of the box.
 313  * \param[in] (ymin,ymax) the minimum and maximum y coordinates of the box.
 314  * \param[out] (px,py) the periodic displacement vector for the first block
 315  *                     to be tested.
 316  * \return The index of the first block to be tested. */
 317 int voropp_loop_2d::init(double xmin,double xmax,double ymin,double ymax,double &px,double &py) {
 318         ai=step_int((xmin-ax)*xsp);
 319         bi=step_int((xmax-ax)*xsp);
 320         if(!xperiodic) {
 321                 if(ai<0) {ai=0;if(bi<0) bi=0;}
 322                 if(bi>=nx) {bi=nx-1;if(ai>=nx) ai=nx-1;}
 323         }
 324         aj=step_int((ymin-ay)*ysp);
 325         bj=step_int((ymax-ay)*ysp);
 326         if(!yperiodic) {
 327                 if(aj<0) {aj=0;if(bj<0) bj=0;}
 328                 if(bj>=ny) {bj=ny-1;if(aj>=ny) aj=ny-1;}
 329         }
 330         i=ai;j=aj;
 331         aip=ip=step_mod(i,nx);apx=px=step_div(i,nx)*boxx;
 332         ajp=jp=step_mod(j,ny);apy=py=step_div(j,ny)*boxy;
 333         inc1=aip-step_mod(bi,nx)+nx;
 334         s=aip+nx*ajp;
 335         return s;
 336 }
 337
 338 /** Returns the next block to be tested in a loop, and updates the periodicity
 339  * vector if necessary.
 340  * \param[in,out] (px,py) the current block on entering the function, which is
 341  *                        updated to the next block on exiting the function. */
 342 int voropp_loop_2d::inc(double &px,double &py) {
 343         if(i<bi) {
 344                 i++;
 345                 if(ip<nx-1) {ip++;s++;} else {ip=0;s+=1-nx;px+=boxx;}
 346                 return s;
 347         } else if(j<bj) {
 348                 i=ai;ip=aip;px=apx;j++;
 349                 if(jp<ny-1) {jp++;s+=inc1;} else {jp=0;s+=inc1-nxy;py+=boxy;}
 350                 return s;
 351         } else return -1;
 352 }