00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include "cell.hh"
00012 #include "container.hh"
00013
00014
00015
00016
00017
00018
00019 template<class r_option>
00020 container_base<r_option>::container_base(fpoint xa,fpoint xb,fpoint ya,fpoint yb,fpoint za,fpoint zb,
00021 int xn,int yn,int zn,bool xper,bool yper,bool zper,int memi)
00022 : ax(xa),bx(xb),ay(ya),by(yb),az(za),bz(zb),
00023 xsp(xn/(xb-xa)),ysp(yn/(yb-ya)),zsp(zn/(zb-za)),
00024 nx(xn),ny(yn),nz(zn),nxy(xn*yn),nxyz(xn*yn*zn),
00025 hx(xper?2*xn+1:xn),hy(yper?2*yn+1:yn),hz(zper?2*zn+1:zn),hxy(hx*hy),hxyz(hx*hy*hz),
00026 xperiodic(xper),yperiodic(yper),zperiodic(zper),
00027 mv(0),wall_number(0),current_wall_size(init_wall_size),
00028 radius(this),sz(radius.mem_size) {
00029 int l;
00030 co=new int[nxyz];
00031 for(l=0;l<nxyz;l++) co[l]=0;
00032 mem=new int[nxyz];
00033 for(l=0;l<nxyz;l++) mem[l]=memi;
00034 id=new int*[nxyz];
00035 for(l=0;l<nxyz;l++) id[l]=new int[memi];
00036 p=new fpoint*[nxyz];
00037 for(l=0;l<nxyz;l++) p[l]=new fpoint[sz*memi];
00038 mask=new unsigned int[hxyz];
00039 for(l=0;l<hxyz;l++) mask[l]=0;
00040 s_size=3*(hxy+hz*(hx+hy));if(s_size<18) s_size=18;
00041 sl=new int[s_size];
00042 walls=new wall*[current_wall_size];
00043 mrad=new fpoint[hgridsq*seq_length];
00044 initialize_radii();
00045 }
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056 template<class r_option>
00057 void container_base<r_option>::initialize_radii() {
00058 const unsigned int b1=1<<21,b2=1<<22,b3=1<<24,b4=1<<25,b5=1<<27,b6=1<<28;
00059 const fpoint xstep=(bx-ax)/nx/fgrid;
00060 const fpoint ystep=(by-ay)/ny/fgrid;
00061 const fpoint zstep=(bz-az)/nz/fgrid;
00062 int i,j,k,lx,ly,lz,l=0,q;
00063 unsigned int *e=const_cast<unsigned int*> (wl);
00064 fpoint xlo,ylo,zlo,xhi,yhi,zhi,minr;fpoint *radp=mrad;
00065 unsigned int f;
00066 for(zlo=0,zhi=zstep,lz=0;lz<hgrid;zlo=zhi,zhi+=zstep,lz++) {
00067 for(ylo=0,yhi=ystep,ly=0;ly<hgrid;ylo=yhi,yhi+=ystep,ly++) {
00068 for(xlo=0,xhi=xstep,lx=0;lx<hgrid;xlo=xhi,xhi+=xstep,l++,lx++) {
00069 minr=large_number;
00070 for(q=e[0]+1;q<seq_length;q++) {
00071 f=e[q];
00072 i=(f&127)-64;
00073 j=(f>>7&127)-64;
00074 k=(f>>14&127)-64;
00075 if((f&b2)==b2) {
00076 compute_minimum(minr,xlo,xhi,ylo,yhi,zlo,zhi,i-1,j,k);
00077 if((f&b1)==0) compute_minimum(minr,xlo,xhi,ylo,yhi,zlo,zhi,i+1,j,k);
00078 } else if((f&b1)==b1) compute_minimum(minr,xlo,xhi,ylo,yhi,zlo,zhi,i+1,j,k);
00079 if((f&b4)==b4) {
00080 compute_minimum(minr,xlo,xhi,ylo,yhi,zlo,zhi,i,j-1,k);
00081 if((f&b3)==0) compute_minimum(minr,xlo,xhi,ylo,yhi,zlo,zhi,i,j+1,k);
00082 } else if((f&b3)==b3) compute_minimum(minr,xlo,xhi,ylo,yhi,zlo,zhi,i,j+1,k);
00083 if((f&b6)==b6) {
00084 compute_minimum(minr,xlo,xhi,ylo,yhi,zlo,zhi,i,j,k-1);
00085 if((f&b5)==0) compute_minimum(minr,xlo,xhi,ylo,yhi,zlo,zhi,i,j,k+1);
00086 } else if((f&b5)==b5) compute_minimum(minr,xlo,xhi,ylo,yhi,zlo,zhi,i,j,k+1);
00087 }
00088 q--;
00089 while(q>0) {
00090 radp[q]=minr;
00091 f=e[q];
00092 i=(f&127)-64;
00093 j=(f>>7&127)-64;
00094 k=(f>>14&127)-64;
00095 compute_minimum(minr,xlo,xhi,ylo,yhi,zlo,zhi,i,j,k);
00096 q--;
00097 }
00098 radp[0]=minr;
00099 e+=seq_length;
00100 radp+=seq_length;
00101 }
00102 }
00103 }
00104 }
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116 template<class r_option>
00117 inline void container_base<r_option>::compute_minimum(fpoint &minr,fpoint &xlo,fpoint &xhi,fpoint &ylo,fpoint &yhi,fpoint &zlo,fpoint &zhi,int ti,int tj,int tk) {
00118 const fpoint boxx=(bx-ax)/nx,boxy=(by-ay)/ny,boxz=(bz-az)/nz;
00119 fpoint radsq,temp;
00120 if(ti>0) {temp=boxx*ti-xhi;radsq=temp*temp;}
00121 else if(ti<0) {temp=xlo-boxx*(1+ti);radsq=temp*temp;}
00122 else radsq=0;
00123
00124 if(tj>0) {temp=boxy*tj-yhi;radsq+=temp*temp;}
00125 else if(tj<0) {temp=ylo-boxy*(1+tj);radsq+=temp*temp;}
00126
00127 if(tk>0) {temp=boxz*tk-zhi;radsq+=temp*temp;}
00128 else if(tk<0) {temp=zlo-boxz*(1+tk);radsq+=temp*temp;}
00129
00130 if(radsq<minr) minr=radsq;
00131 }
00132
00133
00134 template<class r_option>
00135 container_base<r_option>::~container_base() {
00136 int l;
00137 for(l=0;l<nxyz;l++) delete [] p[l];
00138 for(l=0;l<nxyz;l++) delete [] id[l];
00139 delete [] p;
00140 delete [] id;
00141 delete [] mem;
00142 delete [] co;
00143 delete [] mask;
00144 delete [] sl;
00145 delete [] walls;
00146 delete [] mrad;
00147 }
00148
00149
00150 template<class r_option>
00151 void container_base<r_option>::draw_particles(ostream &os) {
00152 int c,l,i;
00153 for(l=0;l<nxyz;l++) {
00154 for(c=0;c<co[l];c++) {
00155 os << id[l][c];
00156 for(i=sz*c;i<sz*(c+1);i++) os << " " << p[l][i];
00157 os << "\n";
00158 }
00159 }
00160 }
00161
00162
00163
00164 template<class r_option>
00165 void container_base<r_option>::draw_particles() {
00166 draw_particles(cin);
00167 }
00168
00169
00170
00171
00172 template<class r_option>
00173 void container_base<r_option>::draw_particles(const char *filename) {
00174 ofstream os;
00175 os.open(filename,ofstream::out|ofstream::trunc);
00176 draw_particles(os);
00177 os.close();
00178 }
00179
00180
00181 template<class r_option>
00182 void container_base<r_option>::draw_particles_pov(ostream &os) {
00183 int l,c;
00184 for(l=0;l<nxyz;l++) {
00185 for(c=0;c<co[l];c++) {
00186 os << "// id " << id[l][c] << "\n";
00187 os << "sphere{<" << p[l][sz*c] << "," << p[l][sz*c+1] << ","
00188 << p[l][sz*c+2] << ">,";
00189 radius.rad(os,l,c);
00190 os << "}\n";
00191 }
00192 }
00193 }
00194
00195
00196
00197 template<class r_option>
00198 void container_base<r_option>::draw_particles_pov() {
00199 draw_particles_pov(cin);
00200 }
00201
00202
00203
00204
00205 template<class r_option>
00206 void container_base<r_option>::draw_particles_pov(const char *filename) {
00207 ofstream os;
00208 os.open(filename,ofstream::out|ofstream::trunc);
00209 draw_particles_pov(os);
00210 os.close();
00211 }
00212
00213
00214 template<class r_option>
00215 void container_base<r_option>::put(int n,fpoint x,fpoint y,fpoint z) {
00216 if(x>ax&&y>ay&&z>az) {
00217 int i,j,k;
00218 i=int((x-ax)*xsp);j=int((y-ay)*ysp);k=int((z-az)*zsp);
00219 if(i<nx&&j<ny&&k<nz) {
00220 i+=nx*j+nxy*k;
00221 if(co[i]==mem[i]) add_particle_memory(i);
00222 p[i][sz*co[i]]=x;p[i][sz*co[i]+1]=y;p[i][sz*co[i]+2]=z;
00223 radius.store_radius(i,co[i],0.5);
00224 id[i][co[i]++]=n;
00225 }
00226 }
00227 }
00228
00229
00230
00231
00232
00233 template<class r_option>
00234 void container_base<r_option>::put(int n,fpoint x,fpoint y,fpoint z,fpoint r) {
00235 if(x>ax&&y>ay&&z>az) {
00236 int i,j,k;
00237 i=int((x-ax)*xsp);j=int((y-ay)*ysp);k=int((z-az)*zsp);
00238 if(i<nx&&j<ny&&k<nz) {
00239 i+=nx*j+nxy*k;
00240 if(co[i]==mem[i]) add_particle_memory(i);
00241 p[i][sz*co[i]]=x;p[i][sz*co[i]+1]=y;p[i][sz*co[i]+2]=z;
00242 radius.store_radius(i,co[i],r);
00243 id[i][co[i]++]=n;
00244 }
00245 }
00246 }
00247
00248
00249 template<class r_option>
00250 void container_base<r_option>::add_particle_memory(int i) {
00251 int *idp;fpoint *pp;
00252 int l,nmem=2*mem[i];
00253 #if VOROPP_VERBOSE >=3
00254 cerr << "Particle memory in region " << i << " scaled up to " << nmem << endl;
00255 #endif
00256 if(nmem>max_particle_memory) throw fatal_error("Absolute maximum memory allocation exceeded");
00257 idp=new int[nmem];
00258 for(l=0;l<co[i];l++) idp[l]=id[i][l];
00259 pp=new fpoint[sz*nmem];
00260 for(l=0;l<sz*co[i];l++) pp[l]=p[i][l];
00261 mem[i]=nmem;
00262 delete [] id[i];id[i]=idp;
00263 delete [] p[i];p[i]=pp;
00264 }
00265
00266
00267 template<class r_option>
00268 inline void container_base<r_option>::add_list_memory() {
00269 int i,j=0,*ps;
00270 ps=new int[s_size*2];
00271 #if VOROPP_VERBOSE >=2
00272 cerr << "List memory scaled up to " << s_size*2 << endl;
00273 #endif
00274 if(s_start<=s_end) {
00275 for(i=s_start;i<s_end;i++) ps[j++]=sl[i];
00276 } else {
00277 for(i=s_start;i<s_size;i++) ps[j++]=sl[i];
00278 for(i=0;i<s_end;i++) ps[j++]=sl[i];
00279 }
00280 s_size*=2;
00281 s_start=0;s_end=j;
00282 delete [] sl;sl=ps;
00283 }
00284
00285
00286 template<class r_option>
00287 void container_base<r_option>::import(istream &is) {
00288 radius.import(is);
00289 }
00290
00291
00292
00293 template<class r_option>
00294 inline void container_base<r_option>::import() {
00295 import(cin);
00296 }
00297
00298
00299
00300
00301 template<class r_option>
00302 inline void container_base<r_option>::import(const char *filename) {
00303 ifstream is;
00304 is.open(filename,ifstream::in);
00305 import(is);
00306 is.close();
00307 }
00308
00309
00310 template<class r_option>
00311 void container_base<r_option>::region_count() {
00312 int i,j,k,ijk=0;
00313 for(k=0;k<nz;k++) {
00314 for(j=0;j<ny;j++) {
00315 for(i=0;i<nx;i++) cout << "Region (" << i << "," << j << "," << k << "): " << co[ijk++] << " particles" << endl;
00316 }
00317 }
00318 }
00319
00320
00321 template<class r_option>
00322 void container_base<r_option>::clear() {
00323 for(int ijk=0;ijk<nxyz;ijk++) co[ijk]=0;
00324 radius.clear_max();
00325 }
00326
00327
00328
00329
00330 template<class r_option>
00331 void container_base<r_option>::draw_cells_gnuplot(const char *filename,fpoint xmin,fpoint xmax,fpoint ymin,fpoint ymax,fpoint zmin,fpoint zmax) {
00332 fpoint x,y,z,px,py,pz;
00333 voropp_loop l1(this);
00334 int q,s;
00335 voronoicell c;
00336 ofstream os;
00337 os.open(filename,ofstream::out|ofstream::trunc);
00338 s=l1.init(xmin,xmax,ymin,ymax,zmin,zmax,px,py,pz);
00339 do {
00340 for(q=0;q<co[s];q++) {
00341 x=p[s][sz*q]+px;y=p[s][sz*q+1]+py;z=p[s][sz*q+2]+pz;
00342 if(x>xmin&&x<xmax&&y>ymin&&y<ymax&&z>zmin&&z<zmax) {
00343 if (compute_cell(c,l1.ip,l1.jp,l1.kp,s,q,x,y,z)) c.draw_gnuplot(os,x,y,z);
00344 }
00345 }
00346 } while((s=l1.inc(px,py,pz))!=-1);
00347 os.close();
00348 }
00349
00350
00351
00352 template<class r_option>
00353 void container_base<r_option>::draw_cells_gnuplot(const char *filename) {
00354 draw_cells_gnuplot(filename,ax,bx,ay,by,az,bz);
00355 }
00356
00357
00358
00359
00360 template<class r_option>
00361 void container_base<r_option>::draw_cells_pov(const char *filename,fpoint xmin,fpoint xmax,fpoint ymin,fpoint ymax,fpoint zmin,fpoint zmax) {
00362 fpoint x,y,z,px,py,pz;
00363 voropp_loop l1(this);
00364 int q,s;
00365 voronoicell c;
00366 ofstream os;
00367 os.open(filename,ofstream::out|ofstream::trunc);
00368 s=l1.init(xmin,xmax,ymin,ymax,zmin,zmax,px,py,pz);
00369 do {
00370 for(q=0;q<co[s];q++) {
00371 os << "// cell " << id[s][q] << "\n";
00372 x=p[s][sz*q]+px;y=p[s][sz*q+1]+py;z=p[s][sz*q+2]+pz;
00373 if(x>xmin&&x<xmax&&y>ymin&&y<ymax&&z>zmin&&z<zmax) {
00374 if (compute_cell(c,l1.ip,l1.jp,l1.kp,s,q,x,y,z)) c.draw_pov(os,x,y,z);
00375 }
00376 }
00377 } while((s=l1.inc(px,py,pz))!=-1);
00378 os.close();
00379 }
00380
00381
00382
00383 template<class r_option>
00384 void container_base<r_option>::draw_cells_pov(const char *filename) {
00385 draw_cells_pov(filename,ax,bx,ay,by,az,bz);
00386 }
00387
00388
00389
00390
00391
00392 template<class r_option>
00393 void container_base<r_option>::compute_all_cells() {
00394 voronoicell c;
00395 int i,j,k,ijk=0,q;
00396 for(k=0;k<nz;k++) for(j=0;j<ny;j++) for(i=0;i<nx;i++,ijk++) {
00397 for(q=0;q<co[ijk];q++) compute_cell(c,i,j,k,ijk,q);
00398 }
00399 }
00400
00401
00402
00403
00404 template<class r_option>
00405 void container_base<r_option>::store_cell_volumes(fpoint *bb) {
00406 voronoicell c;
00407 int i,j,k,ijk=0,q;
00408 for(k=0;k<nz;k++) for(j=0;j<ny;j++) for(i=0;i<nx;i++,ijk++) {
00409 for(q=0;q<co[ijk];q++) {
00410 if (compute_cell(c,i,j,k,ijk,q)) {
00411 bb[id[ijk][q]]=c.volume();
00412 } else {
00413 bb[id[ijk][q]]=0;
00414 }
00415 }
00416 }
00417 }
00418
00419
00420
00421
00422
00423
00424
00425
00426 template<class r_option>
00427 fpoint container_base<r_option>::packing_fraction(fpoint *bb,fpoint cx,fpoint cy,fpoint cz,fpoint r) {
00428 voropp_loop l1(this);
00429 fpoint px,py,pz,x,y,z,rsq=r*r,pvol=0,vvol=0;
00430 int q,s;
00431 s=l1.init(cx,cy,cz,r,px,py,pz);
00432 do {
00433 for(q=0;q<co[s];q++) {
00434 x=p[s][sz*q]+px-cx;
00435 y=p[s][sz*q+1]+py-cy;
00436 z=p[s][sz*q+2]+pz-cz;
00437 if (x*x+y*y+z*z<rsq) {
00438 pvol+=radius.volume(s,q);
00439 vvol+=bb[id[s][q]];
00440 }
00441 }
00442 } while((s=l1.inc(px,py,pz))!=-1);
00443 return vvol>tolerance?pvol/vvol*4.1887902047863909846168578443726:0;
00444 }
00445
00446
00447
00448
00449
00450
00451
00452
00453 template<class r_option>
00454 fpoint container_base<r_option>::packing_fraction(fpoint *bb,fpoint xmin,fpoint xmax,fpoint ymin,fpoint ymax,fpoint zmin,fpoint zmax) {
00455 voropp_loop l1(this);
00456 fpoint x,y,z,px,py,pz,pvol=0,vvol=0;
00457 int q,s;
00458 s=l1.init(xmin,xmax,ymin,ymax,zmin,zmax,px,py,pz);
00459 do {
00460 for(q=0;q<co[s];q++) {
00461 x=p[s][sz*q]+px;
00462 y=p[s][sz*q+1]+py;
00463 z=p[s][sz*q+2]+pz;
00464 if(x>xmin&&x<xmax&&y>ymin&&y<ymax&&z>zmin&&z<zmax) {
00465 pvol+=radius.volume(s,q);
00466 vvol+=bb[id[s][q]];
00467 }
00468 }
00469 } while((s=l1.inc(px,py,pz))!=-1);
00470 return vvol>tolerance?pvol/vvol*4.1887902047863909846168578443726:0;
00471 }
00472
00473
00474
00475 template<class r_option>
00476 fpoint container_base<r_option>::sum_cell_volumes() {
00477 voronoicell c;
00478 int i,j,k,ijk=0,q;fpoint vol=0;
00479 for(k=0;k<nz;k++) for(j=0;j<ny;j++) for(i=0;i<nx;i++,ijk++) {
00480 for(q=0;q<co[ijk];q++) if (compute_cell(c,i,j,k,ijk,q)) vol+=c.volume();
00481 }
00482 return vol;
00483 }
00484
00485
00486
00487 template<class r_option>
00488 inline void container_base<r_option>::count_all_faces(ostream &os) {
00489 voronoicell c;
00490 fpoint x,y,z;
00491 int i,j,k,ijk=0,q;
00492 for(k=0;k<nz;k++) for(j=0;j<ny;j++) for(i=0;i<nx;i++,ijk++) {
00493 for(q=0;q<co[ijk];q++) {
00494 x=p[ijk][sz*q];y=p[ijk][sz*q+1];z=p[ijk][sz*q+2];
00495 os << id[ijk][q] << " " << x << " " << y << " " << z;
00496 radius.print(os,ijk,q);
00497 if (compute_cell(c,i,j,k,ijk,q,x,y,z)) {
00498 os << " " << c.number_of_faces() << "\n";
00499 } else os << " 0\n";
00500 }
00501 }
00502 }
00503
00504
00505
00506 template<class r_option>
00507 void container_base<r_option>::count_all_faces() {
00508 count_all_faces(cout);
00509 }
00510
00511
00512
00513
00514 template<class r_option>
00515 void container_base<r_option>::count_all_faces(const char* filename) {
00516 ofstream os;
00517 os.open(filename,ofstream::out|ofstream::trunc);
00518 count_all_faces(os);
00519 os.close();
00520 }
00521
00522
00523
00524 template<class r_option>
00525 template<class n_option>
00526 inline void container_base<r_option>::print_all(ostream &os,voronoicell_base<n_option> &c) {
00527 fpoint x,y,z;
00528 int i,j,k,ijk=0,q;
00529 for(k=0;k<nz;k++) for(j=0;j<ny;j++) for(i=0;i<nx;i++,ijk++) {
00530 for(q=0;q<co[ijk];q++) {
00531 x=p[ijk][sz*q];y=p[ijk][sz*q+1];z=p[ijk][sz*q+2];
00532 os << id[ijk][q] << " " << x << " " << y << " " << z;
00533 radius.print(os,ijk,q);
00534 if (compute_cell(c,i,j,k,ijk,q,x,y,z)) {
00535 os << " " << c.volume();
00536 c.neighbors(os);
00537 os << "\n";
00538 } else os << " 0\n";
00539 }
00540 }
00541 }
00542
00543
00544
00545 template<class r_option>
00546 void container_base<r_option>::print_all(ostream &os) {
00547 voronoicell c;
00548 print_all(os,c);
00549 }
00550
00551
00552 template<class r_option>
00553 void container_base<r_option>::print_all() {
00554 voronoicell c;
00555 print_all(cout,c);
00556 }
00557
00558
00559
00560
00561 template<class r_option>
00562 inline void container_base<r_option>::print_all(const char* filename) {
00563 voronoicell c;
00564 ofstream os;
00565 os.open(filename,ofstream::out|ofstream::trunc);
00566 print_all(os,c);
00567 os.close();
00568 }
00569
00570
00571
00572
00573 template<class r_option>
00574 void container_base<r_option>::print_all_neighbor(ostream &os) {
00575 voronoicell_neighbor c;
00576 print_all(os,c);
00577 }
00578
00579
00580
00581 template<class r_option>
00582 void container_base<r_option>::print_all_neighbor() {
00583 voronoicell_neighbor c;
00584 print_all(cout,c);
00585 }
00586
00587
00588
00589
00590 template<class r_option>
00591 inline void container_base<r_option>::print_all_neighbor(const char* filename) {
00592 voronoicell_neighbor c;
00593 ofstream os;
00594 os.open(filename,ofstream::out|ofstream::trunc);
00595 print_all(os,c);
00596 os.close();
00597 }
00598
00599
00600
00601
00602
00603
00604 template<class r_option>
00605 template<class n_option>
00606 inline bool container_base<r_option>::initialize_voronoicell(voronoicell_base<n_option> &c,fpoint x,fpoint y,fpoint z) {
00607 fpoint x1,x2,y1,y2,z1,z2;
00608 if(xperiodic) x1=-(x2=0.5*(bx-ax));else {x1=ax-x;x2=bx-x;}
00609 if(yperiodic) y1=-(y2=0.5*(by-ay));else {y1=ay-y;y2=by-y;}
00610 if(zperiodic) z1=-(z2=0.5*(bz-az));else {z1=az-z;z2=bz-z;}
00611 c.init(x1,x2,y1,y2,z1,z2);
00612 for(int j=0;j<wall_number;j++) {
00613 if (!(walls[j]->cut_cell(c,x,y,z))) return false;
00614 }
00615 return true;
00616 }
00617
00618
00619
00620
00621
00622
00623 template<class r_option>
00624 bool container_base<r_option>::point_inside(fpoint x,fpoint y,fpoint z) {
00625 if(x<ax||x>bx||y<ay||y>by||z<az||z>bz) return false;
00626 return point_inside_walls(x,y,z);
00627 }
00628
00629
00630
00631
00632
00633
00634
00635 template<class r_option>
00636 bool container_base<r_option>::point_inside_walls(fpoint x,fpoint y,fpoint z) {
00637 for(int j=0;j<wall_number;j++) if(!walls[j]->point_inside(x,y,z)) return false;
00638 return true;
00639 }
00640
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653
00654
00655 template<class r_option>
00656 template<class n_option>
00657 bool container_base<r_option>::compute_cell_sphere(voronoicell_base<n_option> &c,int i,int j,int k,int ijk,int s,fpoint x,fpoint y,fpoint z) {
00658
00659
00660
00661 const fpoint length_scale=1;
00662 fpoint x1,y1,z1,qx,qy,qz,lr=0,lrs=0,ur,urs,rs;
00663 int q,t;
00664 voropp_loop l(this);
00665 if (!initialize_voronoicell(c,x,y,z)) return false;
00666
00667
00668
00669
00670 radius.init(ijk,s);
00671 while(radius.cutoff(lrs)<c.maxradsq()) {
00672 ur=lr+0.5*length_scale;urs=ur*ur;
00673 t=l.init(x,y,z,ur,qx,qy,qz);
00674 do {
00675 for(q=0;q<co[t];q++) {
00676 x1=p[t][sz*q]+qx-x;y1=p[t][sz*q+1]+qy-y;z1=p[t][sz*q+2]+qz-z;
00677 rs=x1*x1+y1*y1+z1*z1;
00678 if(lrs-tolerance<rs&&rs<urs&&(q!=s||ijk!=t)) {
00679 if (!c.nplane(x1,y1,z1,radius.scale(rs,t,q),id[t][q])) return false;
00680 }
00681 }
00682 } while((t=l.inc(qx,qy,qz))!=-1);
00683 lr=ur;lrs=urs;
00684 }
00685 return true;
00686 }
00687
00688
00689
00690
00691
00692
00693
00694
00695
00696 template<class r_option>
00697 template<class n_option>
00698 inline bool container_base<r_option>::compute_cell_sphere(voronoicell_base<n_option> &c,int i,int j,int k,int ijk,int s) {
00699 fpoint x=p[ijk][sz*s],y=p[ijk][sz*s+1],z=p[ijk][sz*s+2];
00700 return compute_cell_sphere(c,i,j,k,ijk,s,x,y,z);
00701 }
00702
00703
00704
00705
00706
00707
00708
00709
00710
00711
00712 template<class r_option>
00713 template<class n_option>
00714 inline bool container_base<r_option>::compute_cell(voronoicell_base<n_option> &c,int i,int j,int k,int ijk,int s) {
00715 fpoint x=p[ijk][sz*s],y=p[ijk][sz*s+1],z=p[ijk][sz*s+2];
00716 return compute_cell(c,i,j,k,ijk,s,x,y,z);
00717 }
00718
00719
00720
00721
00722
00723
00724
00725
00726
00727
00728
00729
00730
00731
00732
00733
00734
00735
00736
00737
00738
00739
00740
00741
00742
00743 template<class r_option>
00744 template<class n_option>
00745 bool container_base<r_option>::compute_cell(voronoicell_base<n_option> &c,int i,int j,int k,int ijk,int s,fpoint x,fpoint y,fpoint z) {
00746 const fpoint boxx=(bx-ax)/nx,boxy=(by-ay)/ny,boxz=(bz-az)/nz;
00747 fpoint x1,y1,z1,qx=0,qy=0,qz=0;
00748 fpoint xlo,ylo,zlo,xhi,yhi,zhi,rs;
00749 int ci,cj,ck,di,dj,dk,dijk,ei,ej,ek,eijk,si,sj,sk,sijk;
00750 fpoint gxs,gys,gzs,*radp;
00751 int f,g,l;unsigned int q,*e;
00752 const unsigned int b1=1<<21,b2=1<<22,b3=1<<24,b4=1<<25,b5=1<<27,b6=1<<28;
00753
00754 radius.init(ijk,s);
00755
00756
00757 if (!initialize_voronoicell(c,x,y,z)) return false;
00758 fpoint crs,mrs;
00759
00760 int next_count=3,list_index=0,list_size=8;
00761 int count_list[]={7,11,15,19,26,35,45,59};
00762
00763
00764 for(l=0;l<s;l++) {
00765 x1=p[ijk][sz*l]-x;
00766 y1=p[ijk][sz*l+1]-y;
00767 z1=p[ijk][sz*l+2]-z;
00768 rs=radius.scale(x1*x1+y1*y1+z1*z1,ijk,l);
00769 if (!c.nplane(x1,y1,z1,rs,id[ijk][l])) return false;
00770 }
00771 l++;
00772 while(l<co[ijk]) {
00773 x1=p[ijk][sz*l]-x;
00774 y1=p[ijk][sz*l+1]-y;
00775 z1=p[ijk][sz*l+2]-z;
00776 rs=radius.scale(x1*x1+y1*y1+z1*z1,ijk,l);
00777 if (!c.nplane(x1,y1,z1,rs,id[ijk][l])) return false;
00778 l++;
00779 }
00780
00781
00782
00783
00784 mrs=c.maxradsq();
00785
00786
00787
00788
00789
00790 unsigned int m1,m2;
00791 fpoint fx=x-ax-boxx*i,fy=y-ay-boxy*j,fz=z-az-boxz*k;
00792 si=int(fx*xsp*fgrid);sj=int(fy*ysp*fgrid);sk=int(fz*zsp*fgrid);
00793
00794
00795
00796
00797
00798
00799
00800
00801 if(si>=hgrid) {
00802 gxs=fx;
00803 m1=127+(3<<21);si=fgrid-1-si;m2=1+(1<<21);
00804 } else {m1=m2=0;gxs=boxx-fx;}
00805 if(sj>=hgrid) {
00806 gys=fy;
00807 m1|=(127<<7)+(3<<24);sj=fgrid-1-sj;m2|=(1<<7)+(1<<24);
00808 } else gys=boxy-fy;
00809 if(sk>=hgrid) {
00810 gzs=fz;
00811 m1|=(127<<14)+(3<<27);sk=fgrid-1-sk;m2|=(1<<14)+(1<<27);
00812 } else gzs=boxz-fz;
00813 gxs*=gxs;gys*=gys;gzs*=gzs;
00814
00815
00816
00817
00818 if(si<0) si=0;if(sj<0) sj=0;if(sk<0) sk=0;
00819
00820
00821
00822 sijk=si+hgrid*(sj+hgrid*sk);
00823 radp=mrad+sijk*seq_length;
00824 e=(const_cast<unsigned int*> (wl))+sijk*seq_length;
00825
00826
00827
00828 f=e[0];g=0;
00829 do {
00830
00831
00832
00833 if(g==next_count) {
00834 mrs=c.maxradsq();
00835 if(list_index!=list_size) next_count=count_list[list_index++];
00836 }
00837
00838
00839
00840 if(mrs<radius.cutoff(radp[g])) return true;
00841 g++;
00842
00843
00844
00845
00846 q=e[g];q^=m1;q+=m2;
00847 di=q&127;di-=64;
00848 dj=(q>>7)&127;dj-=64;
00849 dk=(q>>14)&127;dk-=64;
00850
00851
00852 if(xperiodic) {if(di<-nx) continue;else if(di>nx) continue;}
00853 else {if(di<-i) continue;else if(di>=nx-i) continue;}
00854 if(yperiodic) {if(dj<-ny) continue;else if(dj>ny) continue;}
00855 else {if(dj<-j) continue;else if(dj>=ny-j) continue;}
00856 if(zperiodic) {if(dk<-nz) continue;else if(dk>nz) continue;}
00857 else {if(dk<-k) continue;else if(dk>=nz-k) continue;}
00858
00859
00860
00861
00862
00863
00864 if(compute_min_max_radius(di,dj,dk,fx,fy,fz,gxs,gys,gzs,crs,mrs)) continue;
00865
00866
00867
00868 di+=i;dj+=j;dk+=k;
00869 if(xperiodic) {if(di<0) {qx=ax-bx;di+=nx;} else if(di>=nx) {qx=bx-ax;di-=nx;} else qx=0;}
00870 if(yperiodic) {if(dj<0) {qy=ay-by;dj+=ny;} else if(dj>=ny) {qy=by-ay;dj-=ny;} else qy=0;}
00871 if(zperiodic) {if(dk<0) {qz=az-bz;dk+=nz;} else if(dk>=nz) {qz=bz-az;dk-=nz;} else qz=0;}
00872 dijk=di+nx*(dj+ny*dk);
00873
00874
00875
00876
00877
00878 if(mrs>radius.cutoff(crs)) {
00879 for(l=0;l<co[dijk];l++) {
00880 x1=p[dijk][sz*l]+qx-x;
00881 y1=p[dijk][sz*l+1]+qy-y;
00882 z1=p[dijk][sz*l+2]+qz-z;
00883 rs=radius.scale(x1*x1+y1*y1+z1*z1,dijk,l);
00884 if (!c.nplane(x1,y1,z1,rs,id[dijk][l])) return false;
00885 }
00886 } else {
00887 for(l=0;l<co[dijk];l++) {
00888 x1=p[dijk][sz*l]+qx-x;
00889 y1=p[dijk][sz*l+1]+qy-y;
00890 z1=p[dijk][sz*l+2]+qz-z;
00891 rs=radius.scale(x1*x1+y1*y1+z1*z1,dijk,l);
00892 if(rs<mrs) {
00893 if (!c.nplane(x1,y1,z1,rs,id[dijk][l])) return false;
00894 }
00895 }
00896 }
00897 } while(g<f);
00898
00899
00900
00901
00902
00903
00904
00905
00906
00907 ci=xperiodic?nx:i;
00908 cj=yperiodic?ny:j;
00909 ck=zperiodic?nz:k;
00910
00911
00912
00913 mv++;
00914 if(mv==0) {
00915 for(l=0;l<hxyz;l++) mask[l]=0;
00916 mv=1;
00917 }
00918
00919
00920 s_start=s_end=0;
00921
00922 while(g<seq_length-1) {
00923
00924
00925
00926 if(g==next_count) {
00927 mrs=c.maxradsq();
00928 if(list_index!=list_size) next_count=count_list[list_index++];
00929 }
00930
00931
00932
00933 if(mrs<radius.cutoff(radp[g])) return true;
00934 g++;
00935
00936
00937
00938
00939 q=e[g];q^=m1;q+=m2;
00940 di=q&127;di-=64;
00941 dj=(q>>7)&127;dj-=64;
00942 dk=(q>>14)&127;dk-=64;
00943
00944
00945
00946
00947 ei=ci+di;
00948 ej=cj+dj;
00949 ek=ck+dk;
00950 if(ei<0) continue;else if(ei>=hx) continue;
00951 if(ej<0) continue;else if(ej>=hy) continue;
00952 if(ek<0) continue;else if(ek>=hz) continue;
00953 eijk=ei+hx*(ej+hy*ek);
00954 mask[eijk]=mv;
00955
00956
00957
00958
00959
00960
00961 if(compute_min_max_radius(di,dj,dk,fx,fy,fz,gxs,gys,gzs,crs,mrs)) continue;
00962
00963
00964
00965 di+=i;dj+=j;dk+=k;
00966 if(xperiodic) {if(di<0) {qx=ax-bx;di+=nx;} else if(di>=nx) {qx=bx-ax;di-=nx;} else qx=0;}
00967 if(yperiodic) {if(dj<0) {qy=ay-by;dj+=ny;} else if(dj>=ny) {qy=by-ay;dj-=ny;} else qy=0;}
00968 if(zperiodic) {if(dk<0) {qz=az-bz;dk+=nz;} else if(dk>=nz) {qz=bz-az;dk-=nz;} else qz=0;}
00969 dijk=di+nx*(dj+ny*dk);
00970
00971
00972
00973
00974
00975 if(mrs>radius.cutoff(crs)) {
00976 for(l=0;l<co[dijk];l++) {
00977 x1=p[dijk][sz*l]+qx-x;
00978 y1=p[dijk][sz*l+1]+qy-y;
00979 z1=p[dijk][sz*l+2]+qz-z;
00980 rs=radius.scale(x1*x1+y1*y1+z1*z1,dijk,l);
00981 if (!c.nplane(x1,y1,z1,rs,id[dijk][l])) return false;
00982 }
00983 } else {
00984 for(l=0;l<co[dijk];l++) {
00985 x1=p[dijk][sz*l]+qx-x;
00986 y1=p[dijk][sz*l+1]+qy-y;
00987 z1=p[dijk][sz*l+2]+qz-z;
00988 rs=radius.scale(x1*x1+y1*y1+z1*z1,dijk,l);
00989 if(rs<mrs) {
00990 if(!c.nplane(x1,y1,z1,rs,id[dijk][l])) return false;
00991 }
00992 }
00993 }
00994
00995
00996
00997 if(s_end+18>s_size) add_list_memory();
00998
00999
01000
01001
01002 if((q&b2)==b2) {
01003 if(ei>0) if(mask[eijk-1]!=mv) {mask[eijk-1]=mv;sl[s_end++]=ei-1;sl[s_end++]=ej;sl[s_end++]=ek;}
01004 if((q&b1)==0) if(ei<hx-1) if(mask[eijk+1]!=mv) {mask[eijk+1]=mv;sl[s_end++]=ei+1;sl[s_end++]=ej;sl[s_end++]=ek;}
01005 } else if((q&b1)==b1) {if(ei<hx-1) if(mask[eijk+1]!=mv) {mask[eijk+1]=mv;sl[s_end++]=ei+1;sl[s_end++]=ej;sl[s_end++]=ek;}}
01006 if((q&b4)==b4) {if(ej>0) if(mask[eijk-hx]!=mv) {mask[eijk-hx]=mv;sl[s_end++]=ei;sl[s_end++]=ej-1;sl[s_end++]=ek;}
01007 if((q&b3)==0) if(ej<hy-1) if(mask[eijk+hx]!=mv) {mask[eijk+hx]=mv;sl[s_end++]=ei;sl[s_end++]=ej+1;sl[s_end++]=ek;}
01008 } else if((q&b3)==b3) {if(ej<hy-1) if(mask[eijk+hx]!=mv) {mask[eijk+hx]=mv;sl[s_end++]=ei;sl[s_end++]=ej+1;sl[s_end++]=ek;}}
01009 if((q&b6)==b6) {if(ek>0) if(mask[eijk-hxy]!=mv) {mask[eijk-hxy]=mv;sl[s_end++]=ei;sl[s_end++]=ej;sl[s_end++]=ek-1;}
01010 if((q&b5)==0) if(ek<hz-1) if(mask[eijk+hxy]!=mv) {mask[eijk+hxy]=mv;sl[s_end++]=ei;sl[s_end++]=ej;sl[s_end++]=ek+1;}
01011 } else if((q&b5)==b5) if(ek<hz-1) if(mask[eijk+hxy]!=mv) {mask[eijk+hxy]=mv;sl[s_end++]=ei;sl[s_end++]=ej;sl[s_end++]=ek+1;}
01012 }
01013
01014
01015 if(mrs<radius.cutoff(radp[g])) return true;
01016
01017
01018
01019 fx+=boxx*ci;fy+=boxy*cj;fz+=boxz*ck;
01020
01021
01022
01023
01024 while(s_start!=s_end) {
01025
01026
01027 if(s_start==s_size) s_start=0;
01028
01029
01030
01031 di=sl[s_start++];dj=sl[s_start++];dk=sl[s_start++];
01032 xlo=di*boxx-fx;xhi=xlo+boxx;
01033 ylo=dj*boxy-fy;yhi=ylo+boxy;
01034 zlo=dk*boxz-fz;zhi=zlo+boxz;
01035
01036
01037
01038 if(di>ci) {
01039 if(dj>cj) {
01040 if(dk>ck) {
01041 if(corner_test(c,xlo,ylo,zlo,xhi,yhi,zhi)) continue;
01042 } else if(dk<ck) {
01043 if(corner_test(c,xlo,ylo,zhi,xhi,yhi,zlo)) continue;
01044 } else {
01045 if(edge_z_test(c,xlo,ylo,zlo,xhi,yhi,zhi)) continue;
01046 }
01047 } else if(dj<cj) {
01048 if(dk>ck) {
01049 if(corner_test(c,xlo,yhi,zlo,xhi,ylo,zhi)) continue;
01050 } else if(dk<ck) {
01051 if(corner_test(c,xlo,yhi,zhi,xhi,ylo,zlo)) continue;
01052 } else {
01053 if(edge_z_test(c,xlo,yhi,zlo,xhi,ylo,zhi)) continue;
01054 }
01055 } else {
01056 if(dk>ck) {
01057 if(edge_y_test(c,xlo,ylo,zlo,xhi,yhi,zhi)) continue;
01058 } else if(dk<ck) {
01059 if(edge_y_test(c,xlo,ylo,zhi,xhi,yhi,zlo)) continue;
01060 } else {
01061 if(face_x_test(c,xlo,ylo,zlo,yhi,zhi)) continue;
01062 }
01063 }
01064 } else if(di<ci) {
01065 if(dj>cj) {
01066 if(dk>ck) {
01067 if(corner_test(c,xhi,ylo,zlo,xlo,yhi,zhi)) continue;
01068 } else if(dk<ck) {
01069 if(corner_test(c,xhi,ylo,zhi,xlo,yhi,zlo)) continue;
01070 } else {
01071 if(edge_z_test(c,xhi,ylo,zlo,xlo,yhi,zhi)) continue;
01072 }
01073 } else if(dj<cj) {
01074 if(dk>ck) {
01075 if(corner_test(c,xhi,yhi,zlo,xlo,ylo,zhi)) continue;
01076 } else if(dk<ck) {
01077 if(corner_test(c,xhi,yhi,zhi,xlo,ylo,zlo)) continue;
01078 } else {
01079 if(edge_z_test(c,xhi,yhi,zlo,xlo,ylo,zhi)) continue;
01080 }
01081 } else {
01082 if(dk>ck) {
01083 if(edge_y_test(c,xhi,ylo,zlo,xlo,yhi,zhi)) continue;
01084 } else if(dk<ck) {
01085 if(edge_y_test(c,xhi,ylo,zhi,xlo,yhi,zlo)) continue;
01086 } else {
01087 if(face_x_test(c,xhi,ylo,zlo,yhi,zhi)) continue;
01088 }
01089 }
01090 } else {
01091 if(dj>cj) {
01092 if(dk>ck) {
01093 if(edge_x_test(c,xlo,ylo,zlo,xhi,yhi,zhi)) continue;
01094 } else if(dk<ck) {
01095 if(edge_x_test(c,xlo,ylo,zhi,xhi,yhi,zlo)) continue;
01096 } else {
01097 if(face_y_test(c,xlo,ylo,zlo,xhi,zhi)) continue;
01098 }
01099 } else if(dj<cj) {
01100 if(dk>ck) {
01101 if(edge_x_test(c,xlo,yhi,zlo,xhi,ylo,zhi)) continue;
01102 } else if(dk<ck) {
01103 if(edge_x_test(c,xlo,yhi,zhi,xhi,ylo,zlo)) continue;
01104 } else {
01105 if(face_y_test(c,xlo,yhi,zlo,xhi,zhi)) continue;
01106 }
01107 } else {
01108 if(dk>ck) {
01109 if(face_z_test(c,xlo,ylo,zlo,xhi,yhi)) continue;
01110 } else if(dk<ck) {
01111 if(face_z_test(c,xlo,ylo,zhi,xhi,yhi)) continue;
01112 } else {
01113 cout << "Can't happen ... happened\n";
01114 }
01115 }
01116 }
01117
01118
01119
01120 if(xperiodic) {ei=i+di-nx;if(ei<0) {qx=ax-bx;ei+=nx;} else if(ei>=nx) {qx=bx-ax;ei-=nx;} else qx=0;} else ei=di;
01121 if(yperiodic) {ej=j+dj-ny;if(ej<0) {qy=ay-by;ej+=ny;} else if(ej>=ny) {qy=by-ay;ej-=ny;} else qy=0;} else ej=dj;
01122 if(zperiodic) {ek=k+dk-nz;if(ek<0) {qz=az-bz;ek+=nz;} else if(ek>=nz) {qz=bz-az;ek-=nz;} else qz=0;} else ek=dk;
01123 eijk=ei+nx*(ej+ny*ek);
01124
01125
01126
01127
01128 for(l=0;l<co[eijk];l++) {
01129 x1=p[eijk][sz*l]+qx-x;
01130 y1=p[eijk][sz*l+1]+qy-y;
01131 z1=p[eijk][sz*l+2]+qz-z;
01132 rs=radius.scale(x1*x1+y1*y1+z1*z1,eijk,l);
01133 if (!c.nplane(x1,y1,z1,rs,id[eijk][l])) return false;
01134 }
01135
01136
01137 if((s_start<=s_end?s_size-s_end+s_start:s_end-s_start)<18) add_list_memory();
01138
01139
01140
01141 dijk=di+hx*(dj+hy*dk);
01142 if(di>0) if(mask[dijk-1]!=mv) {if(s_end==s_size) s_end=0;mask[dijk-1]=mv;sl[s_end++]=di-1;sl[s_end++]=dj;sl[s_end++]=dk;}
01143 if(dj>0) if(mask[dijk-hx]!=mv) {if(s_end==s_size) s_end=0;mask[dijk-hx]=mv;sl[s_end++]=di;sl[s_end++]=dj-1;sl[s_end++]=dk;}
01144 if(dk>0) if(mask[dijk-hxy]!=mv) {if(s_end==s_size) s_end=0;mask[dijk-hxy]=mv;sl[s_end++]=di;sl[s_end++]=dj;sl[s_end++]=dk-1;}
01145 if(di<hx-1) if(mask[dijk+1]!=mv) {if(s_end==s_size) s_end=0;mask[dijk+1]=mv;sl[s_end++]=di+1;sl[s_end++]=dj;sl[s_end++]=dk;}
01146 if(dj<hy-1) if(mask[dijk+hx]!=mv) {if(s_end==s_size) s_end=0;mask[dijk+hx]=mv;sl[s_end++]=di;sl[s_end++]=dj+1;sl[s_end++]=dk;}
01147 if(dk<hz-1) if(mask[dijk+hxy]!=mv) {if(s_end==s_size) s_end=0;mask[dijk+hxy]=mv;sl[s_end++]=di;sl[s_end++]=dj;sl[s_end++]=dk+1;}
01148 }
01149
01150 return true;
01151 }
01152
01153
01154
01155
01156
01157
01158
01159
01160
01161
01162 template<class r_option>
01163 template<class n_option>
01164 inline bool container_base<r_option>::corner_test(voronoicell_base<n_option> &c,fpoint xl,fpoint yl,fpoint zl,fpoint xh,fpoint yh,fpoint zh) {
01165 if(c.plane_intersects_guess(xh,yl,zl,radius.cutoff(xl*xh+yl*yl+zl*zl))) return false;
01166 if(c.plane_intersects(xh,yh,zl,radius.cutoff(xl*xh+yl*yh+zl*zl))) return false;
01167 if(c.plane_intersects(xl,yh,zl,radius.cutoff(xl*xl+yl*yh+zl*zl))) return false;
01168 if(c.plane_intersects(xl,yh,zh,radius.cutoff(xl*xl+yl*yh+zl*zh))) return false;
01169 if(c.plane_intersects(xl,yl,zh,radius.cutoff(xl*xl+yl*yl+zl*zh))) return false;
01170 if(c.plane_intersects(xh,yl,zh,radius.cutoff(xl*xh+yl*yl+zl*zh))) return false;
01171 return true;
01172 }
01173
01174
01175
01176
01177
01178
01179
01180
01181
01182
01183
01184
01185 template<class r_option>
01186 template<class n_option>
01187 inline bool container_base<r_option>::edge_x_test(voronoicell_base<n_option> &c,fpoint x0,fpoint yl,fpoint zl,fpoint x1,fpoint yh,fpoint zh) {
01188 if(c.plane_intersects_guess(x0,yl,zh,radius.cutoff(yl*yl+zl*zh))) return false;
01189 if(c.plane_intersects(x1,yl,zh,radius.cutoff(yl*yl+zl*zh))) return false;
01190 if(c.plane_intersects(x1,yl,zl,radius.cutoff(yl*yl+zl*zl))) return false;
01191 if(c.plane_intersects(x0,yl,zl,radius.cutoff(yl*yl+zl*zl))) return false;
01192 if(c.plane_intersects(x0,yh,zl,radius.cutoff(yl*yh+zl*zl))) return false;
01193 if(c.plane_intersects(x1,yh,zl,radius.cutoff(yl*yh+zl*zl))) return false;
01194 return true;
01195 }
01196
01197
01198
01199
01200
01201
01202
01203
01204
01205
01206
01207
01208 template<class r_option>
01209 template<class n_option>
01210 inline bool container_base<r_option>::edge_y_test(voronoicell_base<n_option> &c,fpoint xl,fpoint y0,fpoint zl,fpoint xh,fpoint y1,fpoint zh) {
01211 if(c.plane_intersects_guess(xl,y0,zh,radius.cutoff(xl*xl+zl*zh))) return false;
01212 if(c.plane_intersects(xl,y1,zh,radius.cutoff(xl*xl+zl*zh))) return false;
01213 if(c.plane_intersects(xl,y1,zl,radius.cutoff(xl*xl+zl*zl))) return false;
01214 if(c.plane_intersects(xl,y0,zl,radius.cutoff(xl*xl+zl*zl))) return false;
01215 if(c.plane_intersects(xh,y0,zl,radius.cutoff(xl*xh+zl*zl))) return false;
01216 if(c.plane_intersects(xh,y1,zl,radius.cutoff(xl*xh+zl*zl))) return false;
01217 return true;
01218 }
01219
01220
01221
01222
01223
01224
01225
01226
01227
01228
01229
01230
01231 template<class r_option>
01232 template<class n_option>
01233 inline bool container_base<r_option>::edge_z_test(voronoicell_base<n_option> &c,fpoint xl,fpoint yl,fpoint z0,fpoint xh,fpoint yh,fpoint z1) {
01234 if(c.plane_intersects_guess(xl,yh,z0,radius.cutoff(xl*xl+yl*yh))) return false;
01235 if(c.plane_intersects(xl,yh,z1,radius.cutoff(xl*xl+yl*yh))) return false;
01236 if(c.plane_intersects(xl,yl,z1,radius.cutoff(xl*xl+yl*yl))) return false;
01237 if(c.plane_intersects(xl,yl,z0,radius.cutoff(xl*xl+yl*yl))) return false;
01238 if(c.plane_intersects(xh,yl,z0,radius.cutoff(xl*xh+yl*yl))) return false;
01239 if(c.plane_intersects(xh,yl,z1,radius.cutoff(xl*xh+yl*yl))) return false;
01240 return true;
01241 }
01242
01243
01244
01245
01246
01247
01248
01249
01250
01251 template<class r_option>
01252 template<class n_option>
01253 inline bool container_base<r_option>::face_x_test(voronoicell_base<n_option> &c,fpoint xl,fpoint y0,fpoint z0,fpoint y1,fpoint z1) {
01254 if(c.plane_intersects_guess(xl,y0,z0,radius.cutoff(xl*xl))) return false;
01255 if(c.plane_intersects(xl,y0,z1,radius.cutoff(xl*xl))) return false;
01256 if(c.plane_intersects(xl,y1,z1,radius.cutoff(xl*xl))) return false;
01257 if(c.plane_intersects(xl,y1,z0,radius.cutoff(xl*xl))) return false;
01258 return true;
01259 }
01260
01261
01262
01263
01264
01265
01266
01267
01268
01269 template<class r_option>
01270 template<class n_option>
01271 inline bool container_base<r_option>::face_y_test(voronoicell_base<n_option> &c,fpoint x0,fpoint yl,fpoint z0,fpoint x1,fpoint z1) {
01272 if(c.plane_intersects_guess(x0,yl,z0,radius.cutoff(yl*yl))) return false;
01273 if(c.plane_intersects(x0,yl,z1,radius.cutoff(yl*yl))) return false;
01274 if(c.plane_intersects(x1,yl,z1,radius.cutoff(yl*yl))) return false;
01275 if(c.plane_intersects(x1,yl,z0,radius.cutoff(yl*yl))) return false;
01276 return true;
01277 }
01278
01279
01280
01281
01282
01283
01284
01285
01286
01287 template<class r_option>
01288 template<class n_option>
01289 inline bool container_base<r_option>::face_z_test(voronoicell_base<n_option> &c,fpoint x0,fpoint y0,fpoint zl,fpoint x1,fpoint y1) {
01290 if(c.plane_intersects_guess(x0,y0,zl,radius.cutoff(zl*zl))) return false;
01291 if(c.plane_intersects(x0,y1,zl,radius.cutoff(zl*zl))) return false;
01292 if(c.plane_intersects(x1,y1,zl,radius.cutoff(zl*zl))) return false;
01293 if(c.plane_intersects(x1,y0,zl,radius.cutoff(zl*zl))) return false;
01294 return true;
01295 }
01296
01297
01298
01299 template<class r_option>
01300 voropp_loop::voropp_loop(container_base<r_option> *q) : sx(q->bx-q->ax), sy(q->by-q->ay), sz(q->bz-q->az),
01301 xsp(q->xsp),ysp(q->ysp),zsp(q->zsp),
01302 ax(q->ax),ay(q->ay),az(q->az),
01303 nx(q->nx),ny(q->ny),nz(q->nz),nxy(q->nxy),nxyz(q->nxyz),
01304 xperiodic(q->xperiodic),yperiodic(q->yperiodic),zperiodic(q->zperiodic) {}
01305
01306
01307
01308
01309 inline int voropp_loop::init(fpoint vx,fpoint vy,fpoint vz,fpoint r,fpoint &px,fpoint &py,fpoint &pz) {
01310 ai=step_int((vx-ax-r)*xsp);
01311 bi=step_int((vx-ax+r)*xsp);
01312 if(!xperiodic) {
01313 if(ai<0) {ai=0;if(bi<0) bi=0;}
01314 if(bi>=nx) {bi=nx-1;if(ai>=nx) ai=nx-1;}
01315 }
01316 aj=step_int((vy-ay-r)*ysp);
01317 bj=step_int((vy-ay+r)*ysp);
01318 if(!yperiodic) {
01319 if(aj<0) {aj=0;if(bj<0) bj=0;}
01320 if(bj>=ny) {bj=ny-1;if(aj>=ny) aj=ny-1;}
01321 }
01322 ak=step_int((vz-az-r)*zsp);
01323 bk=step_int((vz-az+r)*zsp);
01324 if(!zperiodic) {
01325 if(ak<0) {ak=0;if(bk<0) bk=0;}
01326 if(bk>=nz) {bk=nz-1;if(ak>=nz) ak=nz-1;}
01327 }
01328 i=ai;j=aj;k=ak;
01329 aip=ip=step_mod(i,nx);apx=px=step_div(i,nx)*sx;
01330 ajp=jp=step_mod(j,ny);apy=py=step_div(j,ny)*sy;
01331 akp=kp=step_mod(k,nz);apz=pz=step_div(k,nz)*sz;
01332 inc1=aip-step_mod(bi,nx);
01333 inc2=nx*(ny+ajp-step_mod(bj,ny))+inc1;
01334 inc1+=nx;
01335 s=aip+nx*(ajp+ny*akp);
01336 return s;
01337 }
01338
01339
01340
01341
01342
01343 inline int voropp_loop::init(fpoint xmin,fpoint xmax,fpoint ymin,fpoint ymax,fpoint zmin,fpoint zmax,fpoint &px,fpoint &py,fpoint &pz) {
01344 ai=step_int((xmin-ax)*xsp);
01345 bi=step_int((xmax-ax)*xsp);
01346 if(!xperiodic) {
01347 if(ai<0) {ai=0;if(bi<0) bi=0;}
01348 if(bi>=nx) {bi=nx-1;if(ai>=nx) ai=nx-1;}
01349 }
01350 aj=step_int((ymin-ay)*ysp);
01351 bj=step_int((ymax-ay)*ysp);
01352 if(!yperiodic) {
01353 if(aj<0) {aj=0;if(bj<0) bj=0;}
01354 if(bj>=ny) {bj=ny-1;if(aj>=ny) aj=ny-1;}
01355 }
01356 ak=step_int((zmin-az)*zsp);
01357 bk=step_int((zmax-az)*zsp);
01358 if(!zperiodic) {
01359 if(ak<0) {ak=0;if(bk<0) bk=0;}
01360 if(bk>=nz) {bk=nz-1;if(ak>=nz) ak=nz-1;}
01361 }
01362 i=ai;j=aj;k=ak;
01363 aip=ip=step_mod(i,nx);apx=px=step_div(i,nx)*sx;
01364 ajp=jp=step_mod(j,ny);apy=py=step_div(j,ny)*sy;
01365 akp=kp=step_mod(k,nz);apz=pz=step_div(k,nz)*sz;
01366 inc1=aip-step_mod(bi,nx);
01367 inc2=nx*(ny+ajp-step_mod(bj,ny))+inc1;
01368 inc1+=nx;
01369 s=aip+nx*(ajp+ny*akp);
01370 return s;
01371 }
01372
01373
01374
01375 inline int voropp_loop::inc(fpoint &px,fpoint &py,fpoint &pz) {
01376 if(i<bi) {
01377 i++;
01378 if(ip<nx-1) {ip++;s++;} else {ip=0;s+=1-nx;px+=sx;}
01379 return s;
01380 } else if(j<bj) {
01381 i=ai;ip=aip;px=apx;j++;
01382 if(jp<ny-1) {jp++;s+=inc1;} else {jp=0;s+=inc1-nxy;py+=sy;}
01383 return s;
01384 } else if(k<bk) {
01385 i=ai;ip=aip;j=aj;jp=ajp;px=apx;py=apy;k++;
01386 if(kp<nz-1) {kp++;s+=inc2;} else {kp=0;s+=inc2-nxyz;pz+=sz;}
01387 return s;
01388 } else return -1;
01389 }
01390
01391
01392
01393
01394 inline int voropp_loop::step_int(fpoint a) {
01395 return a<0?int(a)-1:int(a);
01396 }
01397
01398
01399 inline int voropp_loop::step_mod(int a,int b) {
01400 return a>=0?a%b:b-1-(b-1-a)%b;
01401 }
01402
01403
01404 inline int voropp_loop::step_div(int a,int b) {
01405 return a>=0?a/b:-1+(a+1)/b;
01406 }
01407
01408
01409
01410 template<class r_option>
01411 void container_base<r_option>::add_wall(wall& w) {
01412 if(wall_number==current_wall_size) {
01413 current_wall_size*=2;
01414 if(current_wall_size>max_wall_size) throw fatal_error("Wall memory allocation exceeded absolute maximum");
01415 wall **pwall;
01416 pwall=new wall*[current_wall_size];
01417 for(int i=0;i<wall_number;i++) pwall[i]=walls[i];
01418 delete [] walls;walls=pwall;
01419 }
01420 walls[wall_number++]=&w;
01421 }
01422
01423
01424
01425
01426
01427 inline void radius_poly::store_radius(int i,int j,fpoint r) {
01428 cc->p[i][4*j+3]=r;
01429 if(r>max_radius) max_radius=r;
01430 }
01431
01432
01433 inline void radius_poly::clear_max() {
01434 max_radius=0;
01435 }
01436
01437
01438
01439
01440 inline void radius_mono::import(istream &is) {
01441 int n;fpoint x,y,z;
01442 is >> n >> x >> y >> z;
01443 while(!is.eof()) {
01444 cc->put(n,x,y,z);
01445 is >> n >> x >> y >> z;
01446 }
01447 }
01448
01449
01450
01451
01452 inline void radius_poly::import(istream &is) {
01453 int n;fpoint x,y,z,r;
01454 is >> n >> x >> y >> z >> r;
01455 while(!is.eof()) {
01456 cc->put(n,x,y,z,r);
01457 is >> n >> x >> y >> z >> r;
01458 }
01459 }
01460
01461
01462
01463
01464
01465
01466 inline void radius_poly::init(int ijk,int s) {
01467 crad=cc->p[ijk][4*s+3];
01468 mul=1+(crad*crad-max_radius*max_radius)/((max_radius+crad)*(max_radius+crad));
01469 crad*=crad;
01470 }
01471
01472
01473
01474
01475
01476
01477
01478 inline fpoint radius_poly::cutoff(fpoint lrs) {
01479 return mul*lrs;
01480 }
01481
01482
01483
01484
01485
01486
01487 inline fpoint radius_mono::cutoff(fpoint lrs) {
01488 return lrs;
01489 }
01490
01491
01492
01493
01494
01495 inline void radius_mono::rad(ostream &os,int l,int c) {
01496 os << "s";
01497 }
01498
01499
01500
01501
01502
01503 inline void radius_poly::rad(ostream &os,int l,int c) {
01504 os << cc->p[l][4*c+3];
01505 }
01506
01507
01508
01509
01510
01511
01512
01513
01514 inline fpoint radius_mono::volume(int ijk,int s) {
01515 return 0.125;
01516 }
01517
01518
01519
01520
01521
01522 inline fpoint radius_poly::volume(int ijk,int s) {
01523 fpoint a=cc->p[ijk][4*s+3];
01524 return a*a*a;
01525 }
01526
01527
01528
01529
01530
01531
01532
01533 inline fpoint radius_poly::scale(fpoint rs,int t,int q) {
01534 return rs+crad-cc->p[t][4*q+3]*cc->p[t][4*q+3];
01535 }
01536
01537
01538
01539
01540
01541
01542 inline fpoint radius_mono::scale(fpoint rs,int t,int q) {
01543 return rs;
01544 }
01545
01546
01547
01548
01549
01550 inline void radius_poly::print(ostream &os,int ijk,int q) {
01551 os << " " << cc->p[ijk][4*q+3];
01552 }
01553
01554
01555
01556
01557
01558
01559
01560
01561
01562
01563
01564
01565
01566
01567 template<class r_option>
01568 inline bool container_base<r_option>::compute_min_max_radius(int di,int dj,int dk,fpoint fx,fpoint fy,fpoint fz,fpoint gxs,fpoint gys,fpoint gzs,fpoint &crs,fpoint mrs) {
01569 fpoint xlo,ylo,zlo;
01570 const fpoint boxx=(bx-ax)/nx,boxy=(by-ay)/ny,boxz=(bz-az)/nz;
01571 const fpoint bxsq=boxx*boxx+boxy*boxy+boxz*boxz;
01572 if(di>0) {
01573 xlo=di*boxx-fx;
01574 crs=xlo*xlo;
01575 if(dj>0) {
01576 ylo=dj*boxy-fy;
01577 crs+=ylo*ylo;
01578 if(dk>0) {
01579 zlo=dk*boxz-fz;
01580 crs+=zlo*zlo;if(radius.cutoff(crs)>mrs) return true;
01581 crs+=bxsq+2*(boxx*xlo+boxy*ylo+boxz*zlo);
01582 } else if(dk<0) {
01583 zlo=(dk+1)*boxz-fz;
01584 crs+=zlo*zlo;if(radius.cutoff(crs)>mrs) return true;
01585 crs+=bxsq+2*(boxx*xlo+boxy*ylo-boxz*zlo);
01586 } else {
01587 if(radius.cutoff(crs)>mrs) return true;
01588 crs+=boxx*(2*xlo+boxx)+boxy*(2*ylo+boxy)+gzs;
01589 }
01590 } else if(dj<0) {
01591 ylo=(dj+1)*boxy-fy;
01592 crs+=ylo*ylo;
01593 if(dk>0) {
01594 zlo=dk*boxz-fz;
01595 crs+=zlo*zlo;if(radius.cutoff(crs)>mrs) return true;
01596 crs+=bxsq+2*(boxx*xlo-boxy*ylo+boxz*zlo);
01597 } else if(dk<0) {
01598 zlo=(dk+1)*boxz-fz;
01599 crs+=zlo*zlo;if(radius.cutoff(crs)>mrs) return true;
01600 crs+=bxsq+2*(boxx*xlo-boxy*ylo-boxz*zlo);
01601 } else {
01602 if(radius.cutoff(crs)>mrs) return true;
01603 crs+=boxx*(2*xlo+boxx)+boxy*(-2*ylo+boxy)+gzs;
01604 }
01605 } else {
01606 if(dk>0) {
01607 zlo=dk*boxz-fz;
01608 crs+=zlo*zlo;if(radius.cutoff(crs)>mrs) return true;
01609 crs+=boxz*(2*zlo+boxz);
01610 } else if(dk<0) {
01611 zlo=(dk+1)*boxz-fz;
01612 crs+=zlo*zlo;if(radius.cutoff(crs)>mrs) return true;
01613 crs+=boxz*(-2*zlo+boxz);
01614 } else {
01615 if(radius.cutoff(crs)>mrs) return true;
01616 crs+=gzs;
01617 }
01618 crs+=gys+boxx*(2*xlo+boxx);
01619 }
01620 } else if(di<0) {
01621 xlo=(di+1)*boxx-fx;
01622 crs=xlo*xlo;
01623 if(dj>0) {
01624 ylo=dj*boxy-fy;
01625 crs+=ylo*ylo;
01626 if(dk>0) {
01627 zlo=dk*boxz-fz;
01628 crs+=zlo*zlo;if(radius.cutoff(crs)>mrs) return true;
01629 crs+=bxsq+2*(-boxx*xlo+boxy*ylo+boxz*zlo);
01630 } else if(dk<0) {
01631 zlo=(dk+1)*boxz-fz;
01632 crs+=zlo*zlo;if(radius.cutoff(crs)>mrs) return true;
01633 crs+=bxsq+2*(-boxx*xlo+boxy*ylo-boxz*zlo);
01634 } else {
01635 if(radius.cutoff(crs)>mrs) return true;
01636 crs+=boxx*(-2*xlo+boxx)+boxy*(2*ylo+boxy)+gzs;
01637 }
01638 } else if(dj<0) {
01639 ylo=(dj+1)*boxy-fy;
01640 crs+=ylo*ylo;
01641 if(dk>0) {
01642 zlo=dk*boxz-fz;
01643 crs+=zlo*zlo;if(radius.cutoff(crs)>mrs) return true;
01644 crs+=bxsq+2*(-boxx*xlo-boxy*ylo+boxz*zlo);
01645 } else if(dk<0) {
01646 zlo=(dk+1)*boxz-fz;
01647 crs+=zlo*zlo;if(radius.cutoff(crs)>mrs) return true;
01648 crs+=bxsq+2*(-boxx*xlo-boxy*ylo-boxz*zlo);
01649 } else {
01650 if(radius.cutoff(crs)>mrs) return true;
01651 crs+=boxx*(-2*xlo+boxx)+boxy*(-2*ylo+boxy)+gzs;
01652 }
01653 } else {
01654 if(dk>0) {
01655 zlo=dk*boxz-fz;
01656 crs+=zlo*zlo;if(radius.cutoff(crs)>mrs) return true;
01657 crs+=boxz*(2*zlo+boxz);
01658 } else if(dk<0) {
01659 zlo=(dk+1)*boxz-fz;
01660 crs+=zlo*zlo;if(radius.cutoff(crs)>mrs) return true;
01661 crs+=boxz*(-2*zlo+boxz);
01662 } else {
01663 if(radius.cutoff(crs)>mrs) return true;
01664 crs+=gzs;
01665 }
01666 crs+=gys+boxx*(-2*xlo+boxx);
01667 }
01668 } else {
01669 if(dj>0) {
01670 ylo=dj*boxy-fy;
01671 crs=ylo*ylo;
01672 if(dk>0) {
01673 zlo=dk*boxz-fz;
01674 crs+=zlo*zlo;if(radius.cutoff(crs)>mrs) return true;
01675 crs+=boxz*(2*zlo+boxz);
01676 } else if(dk<0) {
01677 zlo=(dk+1)*boxz-fz;
01678 crs+=zlo*zlo;if(radius.cutoff(crs)>mrs) return true;
01679 crs+=boxz*(-2*zlo+boxz);
01680 } else {
01681 if(radius.cutoff(crs)>mrs) return true;
01682 crs+=gzs;
01683 }
01684 crs+=boxy*(2*ylo+boxy);
01685 } else if(dj<0) {
01686 ylo=(dj+1)*boxy-fy;
01687 crs=ylo*ylo;
01688 if(dk>0) {
01689 zlo=dk*boxz-fz;
01690 crs+=zlo*zlo;if(radius.cutoff(crs)>mrs) return true;
01691 crs+=boxz*(2*zlo+boxz);
01692 } else if(dk<0) {
01693 zlo=(dk+1)*boxz-fz;
01694 crs+=zlo*zlo;if(radius.cutoff(crs)>mrs) return true;
01695 crs+=boxz*(-2*zlo+boxz);
01696 } else {
01697 if(radius.cutoff(crs)>mrs) return true;
01698 crs+=gzs;
01699 }
01700 crs+=boxy*(-2*ylo+boxy);
01701 } else {
01702 if(dk>0) {
01703 zlo=dk*boxz-fz;crs=zlo*zlo;if(radius.cutoff(crs)>mrs) return true;
01704 crs+=boxz*(2*zlo+boxz);
01705 } else if(dk<0) {
01706 zlo=(dk+1)*boxz-fz;crs=zlo*zlo;if(radius.cutoff(crs)>mrs) return true;
01707 crs+=boxz*(-2*zlo+boxz);
01708 } else {
01709 crs=0;
01710 cout << "Can't happen ... happened\n";
01711 }
01712 crs+=gys;
01713 }
01714 crs+=gxs;
01715 }
01716 return false;
01717 }
01718
01719 #include "worklist.cc"