30 double bxz_,
double byz_,
double bz_,
int nx_,
int ny_,
int nz_,
int init_mem_,
int ps_)
31 :
unitcell(bx_,bxy_,by_,bxz_,byz_,bz_),
voro_base(nx_,ny_,nz_,bx_/nx_,by_/ny_,bz_/nz_),
32 ey(int(max_uv_y*ysp+1)), ez(int(max_uv_z*zsp+1)), wy(ny+ey), wz(nz+ez),
33 oy(ny+2*ey), oz(nz+2*ez), oxyz(nx*oy*oz), id(new int*[oxyz]), p(new double*[oxyz]),
34 co(new int[oxyz]), mem(new int[oxyz]), img(new char[oxyz]), init_mem(init_mem_), ps(ps_) {
38 int *pp=
co;
while(pp<
co+
oxyz) *(pp++)=0;
43 for(k=
ez;k<
wz;k++)
for(j=
ey;j<
wy;j++)
for(i=0;i<
nx;i++) {
53 for(
int l=
oxyz-1;l>=0;l--)
if(
mem[l]>0) {
73 int nx_,
int ny_,
int nz_,
int init_mem_)
74 :
container_periodic_base(bx_,bxy_,by_,bxz_,byz_,bz_,nx_,ny_,nz_,init_mem_,3),
75 vc(*this,2*nx_+1,2*ey+1,2*ez+1) {}
86 int nx_,
int ny_,
int nz_,
int init_mem_)
87 :
container_periodic_base(bx_,bxy_,by_,bxz_,byz_,bz_,nx_,ny_,nz_,init_mem_,4),
88 vc(*this,2*nx_+1,2*ey+1,2*ez+1) {
ppr=
p;}
97 double *pp=
p[ijk]+3*co[ijk]++;
98 *(pp++)=x;*(pp++)=y;*pp=z;
109 double *pp=
p[ijk]+4*co[ijk]++;
110 *(pp++)=x;*(pp++)=y;*(pp++)=z;*pp=r;
124 double *pp=
p[ijk]+3*co[ijk]++;
125 *(pp++)=x;*(pp++)=y;*pp=z;
139 double *pp=
p[ijk]+4*co[ijk]++;
140 *(pp++)=x;*(pp++)=y;*(pp++)=z;*pp=r;
154 double *pp=
p[ijk]+3*co[ijk]++;
155 *(pp++)=x;*(pp++)=y;*pp=z;
169 double *pp=
p[ijk]+4*co[ijk]++;
170 *(pp++)=x;*(pp++)=y;*(pp++)=z;*pp=r;
284 ijk=ci+
nx*(cj+
oy*ck);
298 int ai,aj,ak,ci,cj,ck,ijk;
304 remap(ai,aj,ak,ci,cj,ck,x,y,z,ijk);
305 vc.find_voronoi_cell(x,y,z,ci,cj,ck,ijk,w,mrs);
332 int ai,aj,ak,ci,cj,ck,ijk;
338 remap(ai,aj,ak,ci,cj,ck,x,y,z,ijk);
339 vc.find_voronoi_cell(x,y,z,ci,cj,ck,ijk,w,mrs);
370 int l,nmem(
mem[i]<<1);
371 if(nmem>max_particle_memory)
373 #if VOROPP_VERBOSE >=3
374 fprintf(stderr,
"Particle memory in region %d scaled up to %d\n",i,nmem);
378 int *idp=
new int[nmem];
379 for(l=0;l<
co[i];l++) idp[l]=
id[i][l];
380 double *pp=
new double[
ps*nmem];
381 for(l=0;l<
ps*co[i];l++) pp[l]=
p[i][l];
385 delete []
id[i];
id[i]=idp;
386 delete []
p[i];
p[i]=pp;
397 while((j=fscanf(fp,
"%d %lg %lg %lg",&i,&x,&y,&z))==4)
put(i,x,y,z);
410 while((j=fscanf(fp,
"%d %lg %lg %lg",&i,&x,&y,&z))==4)
put(vo,i,x,y,z);
422 while((j=fscanf(fp,
"%d %lg %lg %lg %lg",&i,&x,&y,&z,&r))==5)
put(i,x,y,z,r);
435 while((j=fscanf(fp,
"%d %lg %lg %lg %lg",&i,&x,&y,&z,&r))==5)
put(vo,i,x,y,z,r);
443 for(k=0;k<
nz;k++)
for(j=0;j<
ny;j++)
for(i=0;i<
nx;i++)
444 printf(
"Region (%d,%d,%d): %d particles\n",i,j,k,*(cop++));
449 for(
int *cop=
co;cop<
co+
nxyz;cop++) *cop=0;
455 for(
int *cop=
co;cop<
co+
nxyz;cop++) *cop=0;
480 FILE *fp=safe_fopen(filename,
"w");
490 FILE *fp=safe_fopen(filename,
"w");
552 double mix,miy,miz,max,may,maz,*pp;
553 for(k=l=0;k<oz;k++) for(j=0;j<oy;j++) for(i=0;i<nx;i++,l++) if(mem[l]>0) {
556 mix=i*
boxx-tolerance;max=mix+
boxx+tolerance;
557 miy=(j-
ey)*
boxy-tolerance;may=miy+
boxy+tolerance;
558 miz=(k-
ez)*
boxz-tolerance;maz=miz+
boxz+tolerance;
562 for(pp=
p[l],c=0;c<co[l];c++,pp+=ps) if(*pp<mix||*pp>max||pp[1]<miy||pp[1]>may||pp[2]<miz||pp[2]>maz)
563 printf(
"%d %d %d %d %f %f %f %f %f %f %f %f %f\n",
564 id[l][c],i,j,k,*pp,pp[1],pp[2],mix,max,miy,may,miz,maz);
578 int fi=qua-quadiv*
nx,fijk=fi+nx*(dj-ima*
ny+
oy*dk);
579 double dis=ima*
bxy+quadiv*
bx,switchx=di*
boxx-ima*
bxy-quadiv*
bx,adis;
582 if((
img[dijk]&1)==0) {
584 odijk=dijk-1;adis=dis;
586 odijk=dijk+nx-1;adis=dis+
bx;
589 for(l=0;l<
co[fijk];l++) {
596 if((
img[dijk]&2)==0) {
600 fijk++;switchx+=
boxx;
603 odijk=dijk-nx+1;adis=dis-
bx;
605 odijk=dijk+1;adis=dis;
608 for(l=0;l<
co[fijk];l++) {
631 int fi=qi-qidiv*
nx,fj=qj-qjdiv*
ny,fijk=fi+nx*(fj+
oy*(dk-ima*
nz)),fijk2;
632 double disy=ima*
byz+qjdiv*
by,switchy=(dj-
ey)*
boxy-ima*
byz-qjdiv*by;
633 double disx=ima*
bxz+qjdiv*bxy+qidiv*
bx,switchx=di*
boxx-ima*
bxz-qjdiv*bxy-qidiv*
bx;
634 double switchx2,disxl,disxr,disx2,disxr2;
636 if(di==0) {dijkl=dijk+nx-1;disxl=disx+
bx;}
637 else {dijkl=dijk-1;disxl=disx;}
639 if(di==nx-1) {dijkr=dijk-nx+1;disxr=disx-
bx;}
640 else {dijkr=dijk+1;disxr=disx;}
644 if((
img[dijk]&1)==0) {
650 for(l=0;l<
co[fijk];l++) {
651 if(
p[fijk][
ps*l+1]>switchy) {
652 if(
p[fijk][
ps*l]>switchx)
put_image(dijk,fijk,l,disx,disy,
bz*ima);
655 if(!y_exist)
continue;
656 if(
p[fijk][
ps*l]>switchx)
put_image(dijk-nx,fijk,l,disx,disy,
bz*ima);
663 if((
img[dijk]&2)==0) {
665 fijk2=fijk+1-
nx;switchx2=switchx+(1-
nx)*
boxx;disx2=disx+
bx;disxr2=disxr+
bx;
667 fijk2=fijk+1;switchx2=switchx+
boxx;disx2=disx;disxr2=disxr;
674 for(l=0;l<
co[fijk2];l++) {
675 if(
p[fijk2][
ps*l+1]>switchy) {
676 if(
p[fijk2][
ps*l]>switchx2)
put_image(dijkr,fijk2,l,disxr2,disy,
bz*ima);
679 if(!y_exist)
continue;
680 if(
p[fijk2][
ps*l]>switchx2)
put_image(dijkr-nx,fijk2,l,disxr2,disy,
bz*ima);
692 int dqidiv=
step_div(qi,nx)-qidiv;qidiv+=dqidiv;
696 disxl+=bxy+bx*dqidiv;
697 disxr+=bxy+bx*dqidiv;
698 switchx-=bxy+bx*dqidiv;
705 if((
img[dijk]&4)==0) {
711 for(l=0;l<
co[fijk];l++) {
712 if(
p[fijk][
ps*l+1]>switchy) {
713 if(!y_exist)
continue;
714 if(
p[fijk][
ps*l]>switchx)
put_image(dijk+nx,fijk,l,disx,disy,
bz*ima);
717 if(
p[fijk][
ps*l]>switchx)
put_image(dijk,fijk,l,disx,disy,
bz*ima);
724 if((
img[dijk]&8)==0) {
726 fijk2=fijk+1-
nx;switchx2=switchx+(1-
nx)*
boxx;disx2=disx+
bx;disxr2=disxr+
bx;
728 fijk2=fijk+1;switchx2=switchx+
boxx;disx2=disx;disxr2=disxr;
735 for(l=0;l<
co[fijk2];l++) {
736 if(
p[fijk2][
ps*l+1]>switchy) {
737 if(!y_exist)
continue;
738 if(
p[fijk2][
ps*l]>switchx2)
put_image(dijkr+nx,fijk2,l,disxr2,disy,
bz*ima);
741 if(
p[fijk2][
ps*l]>switchx2)
put_image(dijkr,fijk2,l,disxr2,disy,
bz*ima);
760 double *p1=
p[reg]+
ps*
co[reg],*p2=
p[fijk]+
ps*l;
764 if(
ps==4) *(++p1)=*(++p2);
765 id[reg][co[reg]++]=
id[fijk][l];
void create_vertical_image(int di, int dj, int dk)
bool find_voronoi_cell(double x, double y, double z, double &rx, double &ry, double &rz, int &pid)
#define VOROPP_MEMORY_ERROR
container_periodic_poly(double bx_, double bxy_, double by_, double bxz_, double byz_, double bz_, int nx_, int ny_, int nz_, int init_mem_)
Structure for holding information about a particle.
void add_particle_memory(int i)
Class containing data structures common across all particle container classes.
void put_image(int reg, int fijk, int l, double dx, double dy, double dz)
void put(int n, double x, double y, double z, double r)
bool find_voronoi_cell(double x, double y, double z, double &rx, double &ry, double &rz, int &pid)
#define VOROPP_FILE_ERROR
void check_compartmentalized()
Class for computation of the unit Voronoi cell associated with a 3D non-rectangular periodic domain...
void create_side_image(int di, int dj, int dk)
void put(int n, double x, double y, double z)
void create_periodic_image(int di, int dj, int dk)
void import(FILE *fp=stdin)
container_periodic(double bx_, double bxy_, double by_, double bxz_, double byz_, double bz_, int nx_, int ny_, int nz_, int init_mem_)
Header file for the container_periodic_base and related classes.
int step_div(int a, int b)
void print_custom(c_loop &vl, const char *format, FILE *fp)
double sum_cell_volumes()
~container_periodic_base()
bool compute_cell(v_cell &c, c_loop &vl)
container_periodic_base(double bx_, double bxy_, double by_, double bxz_, double byz_, double bz_, int nx_, int ny_, int nz_, int init_mem_, int ps)
void put_locate_block(int &ijk, double &x, double &y, double &z)
void remap(int &ai, int &aj, int &ak, int &ci, int &cj, int &ck, double &x, double &y, double &z, int &ijk)
void import(FILE *fp=stdin)
Extension of the voronoicell_base class to represent a Voronoi cell without neighbor information...
A class for looping over all particles in a container_periodic or container_periodic_poly class...
bool compute_cell(v_cell &c, c_loop &vl)
double sum_cell_volumes()
A class for storing ordering information when particles are added to a container. ...
void print_custom(c_loop &vl, const char *format, FILE *fp)
Class for representing a particle system in a 3D periodic non-orthogonal periodic domain...