branches/2d/src/ctr_boundary_2d.cc

   1 // Voro++, a cell-based Voronoi library
   2 //
   3 // Authors  : Chris H. Rycroft (LBL / UC Berkeley)
   4 //            Cody Robert Dance (UC Berkeley)
   5 // Email    : chr@alum.mit.edu
   6 // Date     : August 30th 2011
   7
   8 /** \file ctr_boundary_2d.cc
   9  * \brief Function implementations for the ctr_boundary_2d and related classes. */
  10
  11 #include "ctr_boundary_2d.hh"
  12
  13 namespace voro {
  14
  15 /** The class constructor sets up the geometry of container, initializing the
  16  * minimum and maximum coordinates in each direction, and setting whether each
  17  * direction is periodic or not. It divides the container into a rectangular
  18  * grid of blocks, and allocates memory for each of these for storing particle
  19  * positions and IDs.
  20  * \param[in] (ax_,bx_) the minimum and maximum x coordinates.
  21  * \param[in] (ay_,by_) the minimum and maximum y coordinates.
  22  * \param[in] (nx_,ny_) the number of grid blocks in each of the three
  23  *                      coordinate directions.
  24  * \param[in] (xperiodic_,yperiodic_) flags setting whether the container is
  25  *                                    periodic in each coordinate direction.
  26  * \param[in] init_mem the initial memory allocation for each block.
  27  * \param[in] ps_ the number of floating point entries to store for each
  28  *                particle. */
  29 container_boundary_2d::container_boundary_2d(double ax_,double bx_,double ay_,double by_,
  30                 int nx_,int ny_,bool xperiodic_,bool yperiodic_,int init_mem)
  31         : voro_base_2d(nx_,ny_,(bx_-ax_)/nx_,(by_-ay_)/ny_),
  32         ax(ax_), bx(bx_), ay(ay_), by(by_), xperiodic(xperiodic_), yperiodic(yperiodic_),
  33         id(new int*[nxy]), p(new double*[nxy]), co(new int[nxy]), mem(new int[nxy]),
  34         wid(new int*[nxy]), nlab(new int*[nxy]), plab(new int**[nxy]), bndpts(new int*[nxy]),
  35         boundary_track(-1), edbc(0), edbm(init_boundary_size),
  36         edb(new int[2*edbm]), bnds(new double[2*edbm]), ps(2), soi(NULL),
  37         vc(*this,xperiodic_?2*nx_+1:nx_,yperiodic_?2*ny_+1:ny_) {
  38         int l;
  39 //        totpar=0;
  40         for(l=0;l<nxy;l++) co[l]=0;
  41         for(l=0;l<nxy;l++) mem[l]=init_mem;
  42         for(l=0;l<nxy;l++) id[l]=new int[init_mem];
  43         for(l=0;l<nxy;l++) p[l]=new double[ps*init_mem];
  44         for(l=0;l<nxy;l++) nlab[l]=new int[init_mem];
  45         for(l=0;l<nxy;l++) plab[l]=new int*[init_mem];
  46         for(l=0;l<nxy;l++) bndpts[l]=new int[init_mem];
  47
  48         for(l=0;l<nxy;l++) {wid[l]=new int[init_wall_tag_size+2];*(wid[l])=0;wid[l][1]=init_wall_tag_size;}
  49 }
  50
  51 /** The container destructor frees the dynamically allocated memory. */
  52 container_boundary_2d::~container_boundary_2d() {
  53         int l;
  54
  55         // Clear "sphere of influence" array if it has been allocated
  56         if(soi!=NULL) delete [] soi;
  57
  58         // Deallocate the block-level arrays
  59         for(l=nxy-1;l>=0;l--) delete [] wid[l];
  60         for(l=nxy-1;l>=0;l--) delete [] bndpts[l];
  61         for(l=nxy-1;l>=0;l--) delete [] plab[l];
  62         for(l=nxy-1;l>=0;l--) delete [] nlab[l];
  63         for(l=nxy-1;l>=0;l--) delete [] p[l];
  64         for(l=nxy-1;l>=0;l--) delete [] id[l];
  65
  66         // Delete the two-dimensional arrays
  67         delete [] id;
  68         delete [] p;
  69         delete [] co;
  70         delete [] mem;
  71 }
  72
  73
  74 /** Put a particle into the correct region of the container.
  75  * \param[in] n the numerical ID of the inserted particle.
  76  * \param[in] (x,y) the position vector of the inserted particle. */
  77 void container_boundary_2d::put(int n,double x,double y) {
  78         int ij;
  79         if(put_locate_block(ij,x,y)) {
  80                 //totpar++;
  81                 id[ij][co[ij]]=n;
  82                 if(boundary_track!=-1) {
  83                         bndpts[ij][co[ij]]=edbc;
  84                         register_boundary(x,y);
  85                 } else bndpts[ij][co[ij]]=-1;
  86                 double *pp=p[ij]+2*co[ij]++;
  87                 *(pp++)=x;*pp=y;
  88         }
  89 }
  90
  91 /** Put a particle into the correct region of the container, also recording
  92  * into which region it was stored.
  93  * \param[in] vo the ordering class in which to record the region.
  94  * \param[in] n the numerical ID of the inserted particle.
  95  * \param[in] (x,y) the position vector of the inserted particle. */
  96 void container_boundary_2d::put(particle_order &vo,int n,double x,double y) {
  97         int ij;
  98         if(put_locate_block(ij,x,y)) {
  99                 //totpar++;
 100                 id[ij][co[ij]]=n;
 101                 if(boundary_track!=-1) {
 102                         bndpts[ij][co[ij]]=edbc;
 103                         register_boundary(x,y);
 104                 } else bndpts[ij][co[ij]]=-1;
 105                 vo.add(ij,co[ij]);
 106                 double *pp=p[ij]+2*co[ij]++;
 107                 *(pp++)=x;*pp=y;
 108         }
 109 }
 110
 111 /** This routine takes a particle position vector, tries to remap it into the
 112  * primary domain. If successful, it computes the region into which it can be
 113  * stored and checks that there is enough memory within this region to store
 114  * it.
 115  * \param[out] ij the region index.
 116  * \param[in,out] (x,y) the particle position, remapped into the primary
 117  *                        domain if necessary.
 118  * \return True if the particle can be successfully placed into the container,
 119  * false otherwise. */
 120 inline bool container_boundary_2d::put_locate_block(int &ij,double &x,double &y) {
 121         if(put_remap(ij,x,y)) {
 122                 if(co[ij]==mem[ij]) add_particle_memory(ij);
 123                 return true;
 124         }
 125 #if VOROPP_REPORT_OUT_OF_BOUNDS ==1
 126         fprintf(stderr,"Out of bounds: (x,y)=(%g,%g)\n",x,y);
 127 #endif
 128         return false;
 129 }
 130
 131 /** Takes a particle position vector and computes the region index into which
 132  * it should be stored. If the container is periodic, then the routine also
 133  * maps the particle position to ensure it is in the primary domain. If the
 134  * container is not periodic, the routine bails out.
 135  * \param[out] ij the region index.
 136  * \param[in,out] (x,y) the particle position, remapped into the primary domain
 137  *                      if necessary.
 138  * \return True if the particle can be successfully placed into the container,
 139  * false otherwise. */
 140 inline bool container_boundary_2d::put_remap(int &ij,double &x,double &y) {
 141         int l;
 142
 143         ij=step_int((x-ax)*xsp);
 144         if(xperiodic) {l=step_mod(ij,nx);x+=boxx*(l-ij);ij=l;}
 145         else if(ij<0||ij>=nx) return false;
 146
 147         int j=step_int((y-ay)*ysp);
 148         if(yperiodic) {l=step_mod(j,ny);y+=boxy*(l-j);j=l;}
 149         else if(j<0||j>=ny) return false;
 150
 151         ij+=nx*j;
 152         return true;
 153 }
 154
 155 /** Increase memory for a particular region.
 156  * \param[in] i the index of the region to reallocate. */
 157 void container_boundary_2d::add_particle_memory(int i) {
 158         int l,nmem=mem[i]<<1;
 159
 160         // Carry out a check on the memory allocation size, and
 161         // print a status message if requested
 162         if(nmem>max_particle_memory_2d)
 163                 voro_fatal_error("Absolute maximum memory allocation exceeded",VOROPP_MEMORY_ERROR);
 164 #if VOROPP_VERBOSE >=3
 165         fprintf(stderr,"Particle memory in region %d scaled up to %d\n",i,nmem);
 166 #endif
 167
 168         // Allocate new memory and copy in the contents of the old arrays
 169         int *idp=new int[nmem];
 170         for(l=0;l<co[i];l++) idp[l]=id[i][l];
 171         double *pp=new double[ps*nmem];
 172         for(l=0;l<ps*co[i];l++) pp[l]=p[i][l];
 173         int *nlabp=new int[nmem];
 174         for(l=0;l<co[i];l++) nlabp[l]=nlab[i][l];
 175         int **plabp=new int*[nmem];
 176         for(l=0;l<co[i];l++) plabp[l]=plab[i][l];
 177         int *bndptsp=new int[nmem];
 178         for(l=0;l<co[i];l++) bndptsp[l]=bndpts[i][l];
 179
 180         // Update pointers and delete old arrays
 181         mem[i]=nmem;
 182         delete [] id[i];id[i]=idp;
 183         delete [] p[i];p[i]=pp;
 184         delete [] nlab[i];nlab[i]=nlabp;
 185         delete [] plab[i];plab[i]=plabp;
 186         delete [] bndpts[i];bndpts[i]=bndptsp;
 187 }
 188
 189 /** Outputs the a list of all the container regions along with the number of
 190  * particles stored within each. */
 191 void container_boundary_2d::region_count() {
 192         int i,j,*cop=co;
 193         for(j=0;j<ny;j++) for(i=0;i<nx;i++)
 194                 printf("Region (%d,%d): %d particles\n",i,j,*(cop++));
 195 }
 196
 197 /** Clears a container of particles. */
 198 void container_boundary_2d::clear() {
 199         for(int *cop=co;cop<co+nxy;cop++) *cop=0;
 200 }
 201
 202 /** Computes all the Voronoi cells and saves customized information about them.
 203  * \param[in] format the custom output string to use.
 204  * \param[in] fp a file handle to write to. */
 205 void container_boundary_2d::print_custom(const char *format,FILE *fp) {
 206         c_loop_all_2d vl(*this);
 207         print_custom(vl,format,fp);
 208 }
 209
 210 /** Computes all the Voronoi cells and saves customized information about them.
 211  * \param[in] format the custom output string to use.
 212  * \param[in] filename the name of the file to write to. */
 213 void container_boundary_2d::print_custom(const char *format,const char *filename) {
 214         FILE *fp=safe_fopen(filename,"w");
 215         print_custom(format,fp);
 216         fclose(fp);
 217 }
 218
 219 /** Computes all of the Voronoi cells in the container, but does nothing
 220  * with the output. It is useful for measuring the pure computation time
 221  * of the Voronoi algorithm, without any additional calculations such as
 222  * volume evaluation or cell output. */
 223 void container_boundary_2d::compute_all_cells() {
 224         voronoicell_nonconvex_2d c;
 225         c_loop_all_2d vl(*this);
 226         if(vl.start()) do compute_cell(c,vl);
 227         while(vl.inc());
 228 }
 229
 230 /** Calculates all of the Voronoi cells and sums their volumes. In most cases
 231  * without walls, the sum of the Voronoi cell volumes should equal the volume
 232  * of the container to numerical precision.
 233  * \return The sum of all of the computed Voronoi volumes. */
 234 double container_boundary_2d::sum_cell_areas() {
 235         voronoicell_nonconvex_2d c;
 236         double area=0;
 237         c_loop_all_2d vl(*this);
 238         if(vl.start()) do if(compute_cell(c,vl)) area+=c.area();while(vl.inc());
 239         return area;
 240 }
 241
 242 /** Draws an outline of the domain in gnuplot format.
 243  * \param[in] fp the file handle to write to. */
 244 void container_boundary_2d::draw_domain_gnuplot(FILE *fp) {
 245         fprintf(fp,"%g %g\n%g %g\n%g %g\n%g %g\n%g %g\n",ax,ay,bx,ay,bx,by,ax,by,ax,ay);
 246 }
 247
 248 /** Draws an outline of the domain in POV-Ray format.
 249  * \param[in] fp the file handle to write to. */
 250 void container_boundary_2d::draw_domain_pov(FILE *fp) {
 251         fprintf(fp,"cylinder{<%g,%g,0>,<%g,%g,0>,rr}\n"
 252                    "cylinder{<%g,%g,0>,<%g,%g,0>,rr}\n",ax,ay,bx,ay,ax,by,bx,by);
 253         fprintf(fp,"cylinder{<%g,%g,0>,<%g,%g,0>,rr}\n"
 254                    "cylinder{<%g,%g,0>,<%g,%g,0>,rr}\n",ax,ay,ax,by,bx,ay,bx,by);
 255         fprintf(fp,"sphere{<%g,%g,0>,rr}\nsphere{<%g,%g,0>,rr}\n"
 256                    "sphere{<%g,%g,0>,rr}\nsphere{<%g,%g,0>,rr}\n",ax,ay,bx,ay,ax,by,bx,by);
 257 }
 258
 259 /** This does the additional set-up for non-convex containers. We assume that
 260  * **p, **id, *co, *mem, *bnds, and edbc have already been setup. We then
 261  * proceed to setup **wid, *soi, and THE PROBLEM POINTS BOOLEAN ARRAY.
 262  * This algorithm keeps the importing seperate from the set-up */
 263 void container_boundary_2d::setup(){
 264         double lx,ly,cx,cy,nx,ny;//last (x,y),current (x,y),next (x,y)
 265         int widl=1,maxwid=1,fwid=1,nwid,lwid;
 266         bool first=true;
 267
 268         tmp=tmpp=new int[3*init_temp_label_size];
 269         tmpe=tmp+3*init_temp_label_size;
 270
 271         while(widl!=edbc){
 272                 cx=bnds[2*widl];cy=bnds[2*widl+1];
 273                 nwid=edb[2*widl];lwid=edb[2*widl+1];
 274                 lx=bnds[lwid*2];ly=bnds[lwid*2+1];
 275                 nx=bnds[2*nwid];ny=bnds[2*nwid+1];
 276
 277                 tag_walls(cx,cy,nx,ny,widl);
 278                 semi_circle_labeling(cx,cy,nx,ny,widl);
 279
 280                 //make sure that the cos(angle)>1 and the angle points inward
 281                 //probpts=(lx-cx)*(nx-cx)+(ly-cy)*(ny-cy)>tolerance &&
 282                 //      cross_product(lx-cx,ly-cy,nx-cx,ny-cy);
 283
 284                 widl=edb[2*widl];
 285                 if(widl>maxwid) maxwid=widl;
 286                 if(widl==fwid){
 287                         widl=maxwid+1;
 288                         fwid=widl;
 289                         maxwid++;
 290                         first=false;
 291                 }
 292         }
 293
 294         // The temporary array can now be used to set up the label table
 295         create_label_table();
 296
 297         // Remove temporary array
 298         delete [] tmp;
 299 }
 300
 301 /** Given two points, tags all the computational boxes that the line segment
 302  * specified by the two points
 303  * goes through. param[in] (x1,y1) this is one point
 304  * \param[in] (x2,y2) this is the other point.
 305  * \param[in] wid this is the wall id bnds[2*wid] is the x index of the first
 306  *                vertex in the c-c direction. */
 307 void container_boundary_2d::tag_walls(double x1,double y1,double x2,double y2,int wid_) {
 308
 309         // Find which boxes these points are within
 310         int i1=int((x1-ax)*xsp),j1=int((y1-ay)*ysp);
 311         int i2=int((x2-ax)*xsp),j2=int((y2-ay)*ysp),k,ij;
 312
 313         // Swap to ensure that i1 is smaller than i2
 314         double q,yfac;
 315         if(i2<i1) {
 316                 q=x1;x1=x2;x2=q;q=y1;y1=y2;y2=q;
 317                 k=i1;i1=i2;i2=k;k=j1;j1=j2;j2=k;
 318         }
 319
 320         ij=i1+j1*nx;
 321         if(j1<j2) {
 322                 yfac=1/(y2-y1);
 323                 do {
 324                         j1++;
 325                         q=ay+j1*boxy;
 326                         k=int((((q-y1)*x2+x1*(y2-q))*yfac-ax)*xsp);
 327                         if(k>=nx) k=nx-1;
 328                         tag_line(ij,(j1-1)*nx+k,wid_);
 329                         ij+=nx;
 330                 } while(j1<j2);
 331         } else if(j1>j2) {
 332                 yfac=1/(y2-y1);
 333                 do {
 334                         q=ay+j1*boxy;
 335                         k=int((((q-y1)*x2+x1*(y2-q))*yfac-ax)*xsp);
 336                         if(k>=nx) k=nx-1;
 337                         tag_line(ij,j1*nx+k,wid_);
 338                         ij-=nx;
 339                         j1--;
 340                 } while(j1>j2);
 341         }
 342         tag_line(ij,i2+j2*nx,wid_);
 343 }
 344
 345 void container_boundary_2d::tag_line(int &ij,int ije,int wid_) {
 346         tag(ij,wid_);
 347         while(ij<ije) {
 348                 ij++;
 349                 tag(ij,wid_);
 350         }
 351 }
 352
 353 inline void container_boundary_2d::tag(int ij,int wid_) {
 354         int *&wp(wid[ij]);
 355         if(*wp==wp[1]) {
 356                 int nws=wp[1]<<1;
 357                 if(nws>max_wall_tag_size) voro_fatal_error("Maximum wall tag memory exceeded",VOROPP_MEMORY_ERROR);
 358                 int *np=new int[nws+2];
 359                 *np=*wp;np[1]=nws;
 360                 for(int i=2;i<*wp+2;i++) np[i]=wp[i];
 361                 delete [] wp;
 362                 wp=np;
 363         }
 364         wp[2+(*wp)++]=wid_;
 365 }
 366
 367 /* Tags particles that are within a semicircle (on the appropriate side) of a
 368  * boundary.
 369  * \param[in] (x1,y1) the start point of the wall segment, arranged so that it
 370  *                    is the first point reached in the counter-clockwise
 371  *                    direction.
 372  * \param[in] (x2,y2) the end points of the wall segment. */
 373 void container_boundary_2d::semi_circle_labeling(double x1,double y1,double x2,double y2,int bid) {
 374
 375         double radius=sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2))*0.5,
 376                midx=(x1+x2)*0.5,midy=(y1+y2)*0.5,cpx,cpy;
 377         int ai=int((midx-radius-ax)*xsp),
 378             bi=int((midx+radius-ax)*xsp),
 379             aj=int((midy-radius-ay)*ysp),
 380             bj=int((midy+radius-ay)*ysp),i,j,ij,k;
 381         if(ai<0) ai=0;if(ai>=nx) ai=nx-1;
 382         if(bi<0) bi=0;if(bi>=nx) bi=nx-1;
 383         if(aj<0) aj=0;if(aj>=ny) aj=ny-1;
 384         if(bj<0) bj=0;if(bj>=ny) bj=ny-1;
 385
 386         // Now loop through all the particles in the boxes we found, tagging
 387         // the ones that are within radius of (midx,midy) and are on the
 388         // appropriate side of the wall
 389         for(j=aj;j<=bj;j++) for(i=ai;i<=bi;i++) {
 390                 ij=i+nx*j;
 391                 for(k=0;k<co[ij];k++) {
 392                         cpx=p[ij][2*k];
 393                         cpy=p[ij][2*k+1];
 394                         if((midx-cpx)*(midx-cpx)+(midy-cpy)*(midy-cpy)<=radius*radius&&
 395                         cross_product((x1-x2),(y1-y2),(cpx-x2),(cpy-y2))&&
 396                         (cpx!=x1||cpy==y1)&&(cpx!=x2||cpy!=y2)) {
 397
 398                                 if(tmpp==tmpe) add_temporary_label_memory();
 399                                 *(tmpp++)=ij;
 400                                 *(tmpp++)=k;
 401                                 *(tmpp++)=bid;
 402                         }
 403                 }
 404         }
 405 }
 406
 407 void container_boundary_2d::create_label_table() {
 408         int ij,q,*pp,tlab=0;
 409
 410         // Clear label counters
 411         for(ij=0;ij<nxy;ij++) for(q=0;q<co[ij];q++) nlab[ij][q]=0;
 412
 413         // Increment label counters
 414         for(pp=tmp;pp<tmpp;pp+=3) {nlab[*pp][pp[1]]++;tlab++;}
 415
 416         // Check for case of no labels at all (which may be common)
 417         if(tlab==0) {
 418 #if VOROPP_VERBOSE >=2
 419                 fputs("No labels needed\n",stderr);
 420 #endif
 421                 return;
 422         }
 423
 424         // If there was already a table from a previous call, remove it
 425         if(soi!=NULL) delete [] soi;
 426
 427         // Allocate the label array, and set up pointers from each particle
 428         // to the corresponding location
 429         pp=soi=new int[tlab];
 430         for(ij=0;ij<nxy;ij++) for(q=0;q<co[ij];pp+=nlab[ij][q++]) plab[ij][q]=pp;
 431
 432         // Fill in the label entries
 433         for(pp=tmp;pp<tmpp;pp+=3) *(plab[*pp][pp[1]]++)=pp[2];
 434
 435         // Reset the label pointers
 436         pp=soi;
 437         for(ij=0;ij<nxy;ij++) for(q=0;q<co[ij];pp+=nlab[ij][q++]) plab[ij][q]=pp;
 438 }
 439
 440 /** Draws the boundaries. (Note: this currently assumes that each boundary loop
 441  * is a continuous block in the bnds array, which will be true for the import
 442  * function. However, it may not be true in other cases, in which case this
 443  * routine would have to be extended.) */
 444 void container_boundary_2d::draw_boundary_gnuplot(FILE *fp) {
 445         int i;
 446
 447         for(i=0;i<edbc;i++) {
 448                 fprintf(fp,"%g %g\n",bnds[2*i],bnds[2*i+1]);
 449
 450                 // If a loop is detected, than complete the loop in the output file
 451                 // and insert a newline
 452                 if(edb[2*i]<i) fprintf(fp,"%g %g\n\n",bnds[2*edb[2*i]],bnds[2*edb[2*i]+1]);
 453         }
 454 }
 455
 456 bool container_boundary_2d::point_inside(double x,double y) {
 457         int i=0,j,k=0;
 458         bool sleft,left,nleft;
 459
 460         while(i<edbc) {
 461                 sleft=left=bnds[2*i]<x;
 462                 do {
 463                         j=edb[2*i];
 464                         nleft=j<i?sleft:bnds[2*j]<x;
 465                         if(nleft!=left) {
 466                                 if(left) {
 467                                         if(bnds[2*j+1]*(x-bnds[2*i])+bnds[2*i+1]*(bnds[2*j]-x)<y*(bnds[2*j]-bnds[2*i])) k++;
 468                                 } else {
 469                                         if(bnds[2*j+1]*(x-bnds[2*i])+bnds[2*i+1]*(bnds[2*j]-x)>y*(bnds[2*j]-bnds[2*i])) k--;
 470                                 }
 471                         }
 472                         left=nleft;
 473                         i++;
 474                 } while(j==i);
 475         }
 476
 477 #if VOROPP_VERBOSE >=2
 478         if(k<0) fprintf(stderr,"Negative winding number of %d for (%g,%g)\n",j,x,y);
 479         else if(k>1) fprintf(stderr,"Winding number of %d for (%g,%g)\n",j,x,y);
 480 #endif
 481         return k>0;
 482 }
 483
 484 template<class v_cell_2d>
 485 bool container_boundary_2d::boundary_cuts(v_cell_2d &c,int ij,double x,double y) {
 486         int i,j,k;
 487         double lx,ly,dx,dy,dr;
 488         for(i=2;i<*(wid[ij])+2;i++) {
 489                 j=2*wid[ij][i];k=2*edb[j];
 490                 dx=bnds[k]-bnds[j];dy=bnds[k+1]-bnds[j+1];
 491                 dr=dy*(bnds[j]-x)-dx*(bnds[j+1]-y);
 492                 if(dr<tolerance) continue;
 493                 lx=bnds[j]+bnds[k]-2*x;
 494                 ly=bnds[j+1]+bnds[k+1]-2*y;
 495                 if(lx*lx+ly*ly>dx*dx+dy*dy) continue;
 496                 if(!c.plane(dy,-dx,2*dr)) return false;
 497         }
 498         return true;
 499 }
 500
 501 bool container_boundary_2d::skip(int ij,int q,double x,double y) {
 502         int j;
 503         double wx1,wy1,wx2,wy2,dx,dy,lx,ly;
 504         double cx=p[ij][ps*q],cy=p[ij][ps*q+1];
 505
 506         for(int i=0;i<nlab[ij][q];i++) {
 507                 j=2*plab[ij][q][i];
 508                 wx1=bnds[j];
 509                 wy1=bnds[j+1];
 510                 j=2*edb[j];
 511                 wx2=bnds[j];
 512                 wy2=bnds[j+1];
 513                 dx=wx1-wx2;
 514                 dy=wy1-wy2;
 515                 wx1-=x;wy1-=y;
 516                 if(dx*wy1-dy*wx1>tolerance) {
 517                         lx=cx-x;ly=cy-y;
 518                         if(wx1*ly-wy1*lx>tolerance&&lx*(wy2-y)-ly*(wx2-x)>tolerance) return true;
 519                 }
 520         }
 521         return false;
 522
 523 }
 524
 525 /** Imports a list of particles from an input stream.
 526  * \param[in] fp a file handle to read from. */
 527 void container_boundary_2d::import(FILE *fp) {
 528         int i;
 529         double x,y;
 530         char *buf(new char[512]);
 531
 532         while(fgets(buf,512,fp)!=NULL) {
 533                 if(strcmp(buf,"#Start\n")==0||strcmp(buf,"# Start\n")==0) {
 534
 535                         // Check that two consecutive start tokens haven't been
 536                         // encountered
 537                         if(boundary_track!=-1) voro_fatal_error("File import error - two consecutive start tokens found",VOROPP_FILE_ERROR);
 538                         start_boundary();
 539
 540                 } else if(strcmp(buf,"#End\n")==0||strcmp(buf,"# End\n")==0||
 541                           strcmp(buf,"#End")==0||strcmp(buf,"# End")==0) {
 542
 543                         // Check that two consecutive end tokens haven't been
 544                         // encountered
 545                         if(boundary_track==-1) voro_fatal_error("File import error - found end token without start token",VOROPP_FILE_ERROR);
 546                         end_boundary();
 547                 } else {
 548
 549                         // Try and read three entries from the line
 550                         if(sscanf(buf,"%d %lg %lg",&i,&x,&y)!=3) voro_fatal_error("File import error - can't parse particle information",VOROPP_FILE_ERROR);
 551                         put(i,x,y);
 552                 }
 553         }
 554         if(boundary_track!=-1) voro_fatal_error("File import error - end of file reached without finding end token",VOROPP_FILE_ERROR);
 555
 556         if(!feof(fp)) voro_fatal_error("File import error - error reading string from file",VOROPP_FILE_ERROR);
 557         delete [] buf;
 558 }
 559
 560 void container_boundary_2d::end_boundary() {
 561         if(boundary_track!=edbc) {
 562                 edb[2*boundary_track+1]=edbc-1;
 563                 edb[2*(edbc-1)]=boundary_track;
 564         }
 565         boundary_track=-1;
 566 }
 567
 568 void container_boundary_2d::register_boundary(double x,double y) {
 569         if(edbc==edbm) add_boundary_memory();
 570         if(edbc!=boundary_track) {
 571                 edb[2*edbc-2]=edbc;
 572                 edb[2*edbc+1]=edbc-1;
 573         }
 574         bnds[2*edbc]=x;
 575         bnds[2*(edbc++)+1]=y;
 576 }
 577
 578 /** Increases the size of the temporary label memory. */
 579 void container_boundary_2d::add_temporary_label_memory() {
 580         int size(tmpe-tmp);
 581         size<<=1;
 582         if(size>3*max_temp_label_size)
 583                 voro_fatal_error("Absolute temporary label memory allocation exceeded",VOROPP_MEMORY_ERROR);
 584 #if VOROPP_VERBOSE >=3
 585         fprintf(stderr,"Temporary label memory in region scaled up to %d\n",size);
 586 #endif
 587         int *ntmp(new int[size]),*tp(tmp);tmpp=ntmp;
 588         while(tp<tmpe) *(tmpp++)=*(tp++);
 589         delete [] tmp;
 590         tmp=ntmp;tmpe=tmp+size;
 591 }
 592
 593 /** Increases the memory allocation for the boundary points. */
 594 void container_boundary_2d::add_boundary_memory() {
 595         int i;
 596         edbm<<=1;
 597         if(edbm>max_boundary_size)
 598                 voro_fatal_error("Absolute boundary memory allocation exceeded",VOROPP_MEMORY_ERROR);
 599 #if VOROPP_VERBOSE >=3
 600         fprintf(stderr,"Boundary memory scaled up to %d\n",size);
 601 #endif
 602
 603         // Reallocate the boundary vertex information
 604         double *nbnds(new double[2*edbm]);
 605         for(i=0;i<2*edbc;i++) nbnds[i]=bnds[i];
 606         delete [] nbnds;bnds=nbnds;
 607
 608         // Reallocate the edge information
 609         int *nedb(new int[2*edbm]);
 610         for(i=0;i<2*edbc;i++) nedb[i]=edb[i];
 611         delete [] edb;edb=nedb;
 612 }
 613
 614 // Explicit instantiation
 615 template bool container_boundary_2d::boundary_cuts(voronoicell_nonconvex_2d&,int,double,double);
 616 template bool container_boundary_2d::boundary_cuts(voronoicell_nonconvex_neighbor_2d&,int,double,double);
 617
 618 }