Homec4science
Diffusion lammps (lammm-devel) Browse

lammps/lib/voro++/src/container_prd.cc
706c4ccdaefdlammm-devel

Last Change
Raw File

container_prd.cc
View Options

// Voro++, a 3D cell-based Voronoi library
//
// Author : Chris H. Rycroft (LBL / UC Berkeley)
// Email : chr@alum.mit.edu
// Date : August 30th 2011
/** \file container_prd.cc
* \brief Function implementations for the container_periodic_base and
* related classes. */
#include "container_prd.hh"
namespace voro {
/** The class constructor sets up the geometry of container, initializing the
* minimum and maximum coordinates in each direction, and setting whether each
* direction is periodic or not. It divides the container into a rectangular
* grid of blocks, and allocates memory for each of these for storing particle
* positions and IDs.
* \param[in] (bx_) The x coordinate of the first unit vector.
* \param[in] (bxy_,by_) The x and y coordinates of the second unit vector.
* \param[in] (bxz_,byz_,bz_) The x, y, and z coordinates of the third unit
* vector.
* \param[in] (nx_,ny_,nz_) the number of grid blocks in each of the three
* coordinate directions.
* \param[in] init_mem_ the initial memory allocation for each block.
* \param[in] ps_ the number of floating point entries to store for each
* particle. */
container_periodic_base::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_)
: unitcell(bx_,bxy_,by_,bxz_,byz_,bz_), voro_base(nx_,ny_,nz_,bx_/nx_,by_/ny_,bz_/nz_),
ey(int(max_uv_y*ysp+1)), ez(int(max_uv_z*zsp+1)), wy(ny+ey), wz(nz+ez),
oy(ny+2*ey), oz(nz+2*ez), oxyz(nx*oy*oz), id(new int*[oxyz]), p(new double*[oxyz]),
co(new int[oxyz]), mem(new int[oxyz]), img(new char[oxyz]), init_mem(init_mem_), ps(ps_) {
int i,j,k,l;
// Clear the global arrays
int *pp=co;while(pp<co+oxyz) *(pp++)=0;
pp=mem;while(pp<mem+oxyz) *(pp++)=0;
char *cp=img;while(cp<img+oxyz) *(cp++)=0;
// Set up memory for the blocks in the primary domain
for(k=ez;k<wz;k++) for(j=ey;j<wy;j++) for(i=0;i<nx;i++) {
l=i+nx*(j+oy*k);
mem[l]=init_mem;
id[l]=new int[init_mem];
p[l]=new double[ps*init_mem];
}
}
/** The container destructor frees the dynamically allocated memory. */
container_periodic_base::~container_periodic_base() {
for(int l=oxyz-1;l>=0;l--) if(mem[l]>0) {
delete [] p[l];
delete [] id[l];
}
delete [] img;
delete [] mem;
delete [] co;
delete [] id;
delete [] p;
}
/** The class constructor sets up the geometry of container.
* \param[in] (bx_) The x coordinate of the first unit vector.
* \param[in] (bxy_,by_) The x and y coordinates of the second unit vector.
* \param[in] (bxz_,byz_,bz_) The x, y, and z coordinates of the third unit
* vector.
* \param[in] (nx_,ny_,nz_) the number of grid blocks in each of the three
* coordinate directions.
* \param[in] init_mem_ the initial memory allocation for each block. */
container_periodic::container_periodic(double bx_,double bxy_,double by_,double bxz_,double byz_,double bz_,
int nx_,int ny_,int nz_,int init_mem_)
: container_periodic_base(bx_,bxy_,by_,bxz_,byz_,bz_,nx_,ny_,nz_,init_mem_,3),
vc(*this,2*nx_+1,2*ey+1,2*ez+1) {}
/** The class constructor sets up the geometry of container.
* \param[in] (bx_) The x coordinate of the first unit vector.
* \param[in] (bxy_,by_) The x and y coordinates of the second unit vector.
* \param[in] (bxz_,byz_,bz_) The x, y, and z coordinates of the third unit
* vector.
* \param[in] (nx_,ny_,nz_) the number of grid blocks in each of the three
* coordinate directions.
* \param[in] init_mem_ the initial memory allocation for each block. */
container_periodic_poly::container_periodic_poly(double bx_,double bxy_,double by_,double bxz_,double byz_,double bz_,
int nx_,int ny_,int nz_,int init_mem_)
: container_periodic_base(bx_,bxy_,by_,bxz_,byz_,bz_,nx_,ny_,nz_,init_mem_,4),
vc(*this,2*nx_+1,2*ey+1,2*ez+1) {ppr=p;}
/** Put a particle into the correct region of the container.
* \param[in] n the numerical ID of the inserted particle.
* \param[in] (x,y,z) the position vector of the inserted particle. */
void container_periodic::put(int n,double x,double y,double z) {
int ijk;
put_locate_block(ijk,x,y,z);
id[ijk][co[ijk]]=n;
double *pp=p[ijk]+3*co[ijk]++;
*(pp++)=x;*(pp++)=y;*pp=z;
}
/** Put a particle into the correct region of the container.
* \param[in] n the numerical ID of the inserted particle.
* \param[in] (x,y,z) the position vector of the inserted particle.
* \param[in] r the radius of the particle. */
void container_periodic_poly::put(int n,double x,double y,double z,double r) {
int ijk;
put_locate_block(ijk,x,y,z);
id[ijk][co[ijk]]=n;
double *pp=p[ijk]+4*co[ijk]++;
*(pp++)=x;*(pp++)=y;*(pp++)=z;*pp=r;
if(max_radius<r) max_radius=r;
}
/** Put a particle into the correct region of the container.
* \param[in] n the numerical ID of the inserted particle.
* \param[in] (x,y,z) the position vector of the inserted particle.
* \param[out] (ai,aj,ak) the periodic image displacement that the particle is
* in, with (0,0,0) corresponding to the primary domain.
*/
void container_periodic::put(int n,double x,double y,double z,int &ai,int &aj,int &ak) {
int ijk;
put_locate_block(ijk,x,y,z,ai,aj,ak);
id[ijk][co[ijk]]=n;
double *pp=p[ijk]+3*co[ijk]++;
*(pp++)=x;*(pp++)=y;*pp=z;
}
/** Put a particle into the correct region of the container.
* \param[in] n the numerical ID of the inserted particle.
* \param[in] (x,y,z) the position vector of the inserted particle.
* \param[in] r the radius of the particle.
* \param[out] (ai,aj,ak) the periodic image displacement that the particle is
* in, with (0,0,0) corresponding to the primary domain.
*/
void container_periodic_poly::put(int n,double x,double y,double z,double r,int &ai,int &aj,int &ak) {
int ijk;
put_locate_block(ijk,x,y,z,ai,aj,ak);
id[ijk][co[ijk]]=n;
double *pp=p[ijk]+4*co[ijk]++;
*(pp++)=x;*(pp++)=y;*(pp++)=z;*pp=r;
if(max_radius<r) max_radius=r;
}
/** Put a particle into the correct region of the container, also recording
* into which region it was stored.
* \param[in] vo the ordering class in which to record the region.
* \param[in] n the numerical ID of the inserted particle.
* \param[in] (x,y,z) the position vector of the inserted particle. */
void container_periodic::put(particle_order &vo,int n,double x,double y,double z) {
int ijk;
put_locate_block(ijk,x,y,z);
id[ijk][co[ijk]]=n;
vo.add(ijk,co[ijk]);
double *pp=p[ijk]+3*co[ijk]++;
*(pp++)=x;*(pp++)=y;*pp=z;
}
/** Put a particle into the correct region of the container, also recording
* into which region it was stored.
* \param[in] vo the ordering class in which to record the region.
* \param[in] n the numerical ID of the inserted particle.
* \param[in] (x,y,z) the position vector of the inserted particle.
* \param[in] r the radius of the particle. */
void container_periodic_poly::put(particle_order &vo,int n,double x,double y,double z,double r) {
int ijk;
put_locate_block(ijk,x,y,z);
id[ijk][co[ijk]]=n;
vo.add(ijk,co[ijk]);
double *pp=p[ijk]+4*co[ijk]++;
*(pp++)=x;*(pp++)=y;*(pp++)=z;*pp=r;
if(max_radius<r) max_radius=r;
}
/** Takes a particle position vector and computes the region index into which
* it should be stored. If the container is periodic, then the routine also
* maps the particle position to ensure it is in the primary domain. If the
* container is not periodic, the routine bails out.
* \param[out] ijk the region index.
* \param[in,out] (x,y,z) the particle position, remapped into the primary
* domain if necessary.
* \return True if the particle can be successfully placed into the container,
* false otherwise. */
void container_periodic_base::put_locate_block(int &ijk,double &x,double &y,double &z) {
// Remap particle in the z direction if necessary
int k=step_int(z*zsp);
if(k<0||k>=nz) {
int ak=step_div(k,nz);
z-=ak*bz;y-=ak*byz;x-=ak*bxz;k-=ak*nz;
}
// Remap particle in the y direction if necessary
int j=step_int(y*ysp);
if(j<0||j>=ny) {
int aj=step_div(j,ny);
y-=aj*by;x-=aj*bxy;j-=aj*ny;
}
// Remap particle in the x direction if necessary
ijk=step_int(x*xsp);
if(ijk<0||ijk>=nx) {
int ai=step_div(ijk,nx);
x-=ai*bx;ijk-=ai*nx;
}
// Compute the block index and check memory allocation
j+=ey;k+=ez;
ijk+=nx*(j+oy*k);
if(co[ijk]==mem[ijk]) add_particle_memory(ijk);
}
/** Takes a particle position vector and computes the region index into which
* it should be stored. If the container is periodic, then the routine also
* maps the particle position to ensure it is in the primary domain. If the
* container is not periodic, the routine bails out.
* \param[out] ijk the region index.
* \param[in,out] (x,y,z) the particle position, remapped into the primary
* domain if necessary.
* \param[out] (ai,aj,ak) the periodic image displacement that the particle is
* in, with (0,0,0) corresponding to the primary domain.
* \return True if the particle can be successfully placed into the container,
* false otherwise. */
void container_periodic_base::put_locate_block(int &ijk,double &x,double &y,double &z,int &ai,int &aj,int &ak) {
// Remap particle in the z direction if necessary
int k=step_int(z*zsp);
if(k<0||k>=nz) {
ak=step_div(k,nz);
z-=ak*bz;y-=ak*byz;x-=ak*bxz;k-=ak*nz;
} else ak=0;
// Remap particle in the y direction if necessary
int j=step_int(y*ysp);
if(j<0||j>=ny) {
aj=step_div(j,ny);
y-=aj*by;x-=aj*bxy;j-=aj*ny;
} else aj=0;
// Remap particle in the x direction if necessary
ijk=step_int(x*xsp);
if(ijk<0||ijk>=nx) {
ai=step_div(ijk,nx);
x-=ai*bx;ijk-=ai*nx;
} else ai=0;
// Compute the block index and check memory allocation
j+=ey;k+=ez;
ijk+=nx*(j+oy*k);
if(co[ijk]==mem[ijk]) add_particle_memory(ijk);
}
/** Takes a position vector and remaps it into the primary domain.
* \param[out] (ai,aj,ak) the periodic image displacement that the vector is in,
* with (0,0,0) corresponding to the primary domain.
* \param[out] (ci,cj,ck) the index of the block that the position vector is
* within, once it has been remapped.
* \param[in,out] (x,y,z) the position vector to consider, which is remapped
* into the primary domain during the routine.
* \param[out] ijk the block index that the vector is within. */
inline void container_periodic_base::remap(int &ai,int &aj,int &ak,int &ci,int &cj,int &ck,double &x,double &y,double &z,int &ijk) {
// Remap particle in the z direction if necessary
ck=step_int(z*zsp);
if(ck<0||ck>=nz) {
ak=step_div(ck,nz);
z-=ak*bz;y-=ak*byz;x-=ak*bxz;ck-=ak*nz;
} else ak=0;
// Remap particle in the y direction if necessary
cj=step_int(y*ysp);
if(cj<0||cj>=ny) {
aj=step_div(cj,ny);
y-=aj*by;x-=aj*bxy;cj-=aj*ny;
} else aj=0;
// Remap particle in the x direction if necessary
ci=step_int(x*xsp);
if(ci<0||ci>=nx) {
ai=step_div(ci,nx);
x-=ai*bx;ci-=ai*nx;
} else ai=0;
cj+=ey;ck+=ez;
ijk=ci+nx*(cj+oy*ck);
}
/** Takes a vector and finds the particle whose Voronoi cell contains that
* vector. This is equivalent to finding the particle which is nearest to the
* vector.
* \param[in] (x,y,z) the vector to test.
* \param[out] (rx,ry,rz) the position of the particle whose Voronoi cell
* contains the vector. This may point to a particle in
* a periodic image of the primary domain.
* \param[out] pid the ID of the particle.
* \return True if a particle was found. If the container has no particles,
* then the search will not find a Voronoi cell and false is returned. */
bool container_periodic::find_voronoi_cell(double x,double y,double z,double &rx,double &ry,double &rz,int &pid) {
int ai,aj,ak,ci,cj,ck,ijk;
particle_record w;
double mrs;
// Remap the vector into the primary domain and then search for the
// Voronoi cell that it is within
remap(ai,aj,ak,ci,cj,ck,x,y,z,ijk);
vc.find_voronoi_cell(x,y,z,ci,cj,ck,ijk,w,mrs);
if(w.ijk!=-1) {
// Assemble the position vector of the particle to be returned,
// applying a periodic remapping if necessary
ci+=w.di;if(ci<0||ci>=nx) ai+=step_div(ci,nx);
rx=p[w.ijk][3*w.l]+ak*bxz+aj*bxy+ai*bx;
ry=p[w.ijk][3*w.l+1]+ak*byz+aj*by;
rz=p[w.ijk][3*w.l+2]+ak*bz;
pid=id[w.ijk][w.l];
return true;
}
return false;
}
/** Takes a vector and finds the particle whose Voronoi cell contains that
* vector. Additional wall classes are not considered by this routine.
* \param[in] (x,y,z) the vector to test.
* \param[out] (rx,ry,rz) the position of the particle whose Voronoi cell
* contains the vector. If the container is periodic,
* this may point to a particle in a periodic image of
* the primary domain.
* \param[out] pid the ID of the particle.
* \return True if a particle was found. If the container has no particles,
* then the search will not find a Voronoi cell and false is returned. */
bool container_periodic_poly::find_voronoi_cell(double x,double y,double z,double &rx,double &ry,double &rz,int &pid) {
int ai,aj,ak,ci,cj,ck,ijk;
particle_record w;
double mrs;
// Remap the vector into the primary domain and then search for the
// Voronoi cell that it is within
remap(ai,aj,ak,ci,cj,ck,x,y,z,ijk);
vc.find_voronoi_cell(x,y,z,ci,cj,ck,ijk,w,mrs);
if(w.ijk!=-1) {
// Assemble the position vector of the particle to be returned,
// applying a periodic remapping if necessary
ci+=w.di;if(ci<0||ci>=nx) ai+=step_div(ci,nx);
rx=p[w.ijk][4*w.l]+ak*bxz+aj*bxy+ai*bx;
ry=p[w.ijk][4*w.l+1]+ak*byz+aj*by;
rz=p[w.ijk][4*w.l+2]+ak*bz;
pid=id[w.ijk][w.l];
return true;
}
return false;
}
/** Increase memory for a particular region.
* \param[in] i the index of the region to reallocate. */
void container_periodic_base::add_particle_memory(int i) {
// Handle the case when no memory has been allocated for this block
if(mem[i]==0) {
mem[i]=init_mem;
id[i]=new int[init_mem];
p[i]=new double[ps*init_mem];
return;
}
// Otherwise, double the memory allocation for this block. Carry out a
// check on the memory allocation size, and print a status message if
// requested.
int l,nmem(mem[i]<<1);
if(nmem>max_particle_memory)
voro_fatal_error("Absolute maximum memory allocation exceeded",VOROPP_MEMORY_ERROR);
#if VOROPP_VERBOSE >=3
fprintf(stderr,"Particle memory in region %d scaled up to %d\n",i,nmem);
#endif
// Allocate new memory and copy in the contents of the old arrays
int *idp=new int[nmem];
for(l=0;l<co[i];l++) idp[l]=id[i][l];
double *pp=new double[ps*nmem];
for(l=0;l<ps*co[i];l++) pp[l]=p[i][l];
// Update pointers and delete old arrays
mem[i]=nmem;
delete [] id[i];id[i]=idp;
delete [] p[i];p[i]=pp;
}
/** Import a list of particles from an open file stream into the container.
* Entries of four numbers (Particle ID, x position, y position, z position)
* are searched for. If the file cannot be successfully read, then the routine
* causes a fatal error.
* \param[in] fp the file handle to read from. */
void container_periodic::import(FILE *fp) {
int i,j;
double x,y,z;
while((j=fscanf(fp,"%d %lg %lg %lg",&i,&x,&y,&z))==4) put(i,x,y,z);
if(j!=EOF) voro_fatal_error("File import error",VOROPP_FILE_ERROR);
}
/** Import a list of particles from an open file stream, also storing the order
* of that the particles are read. Entries of four numbers (Particle ID, x
* position, y position, z position) are searched for. If the file cannot be
* successfully read, then the routine causes a fatal error.
* \param[in,out] vo a reference to an ordering class to use.
* \param[in] fp the file handle to read from. */
void container_periodic::import(particle_order &vo,FILE *fp) {
int i,j;
double x,y,z;
while((j=fscanf(fp,"%d %lg %lg %lg",&i,&x,&y,&z))==4) put(vo,i,x,y,z);
if(j!=EOF) voro_fatal_error("File import error",VOROPP_FILE_ERROR);
}
/** Import a list of particles from an open file stream into the container.
* Entries of five numbers (Particle ID, x position, y position, z position,
* radius) are searched for. If the file cannot be successfully read, then the
* routine causes a fatal error.
* \param[in] fp the file handle to read from. */
void container_periodic_poly::import(FILE *fp) {
int i,j;
double x,y,z,r;
while((j=fscanf(fp,"%d %lg %lg %lg %lg",&i,&x,&y,&z,&r))==5) put(i,x,y,z,r);
if(j!=EOF) voro_fatal_error("File import error",VOROPP_FILE_ERROR);
}
/** Import a list of particles from an open file stream, also storing the order
* of that the particles are read. Entries of four numbers (Particle ID, x
* position, y position, z position, radius) are searched for. If the file
* cannot be successfully read, then the routine causes a fatal error.
* \param[in,out] vo a reference to an ordering class to use.
* \param[in] fp the file handle to read from. */
void container_periodic_poly::import(particle_order &vo,FILE *fp) {
int i,j;
double x,y,z,r;
while((j=fscanf(fp,"%d %lg %lg %lg %lg",&i,&x,&y,&z,&r))==5) put(vo,i,x,y,z,r);
if(j!=EOF) voro_fatal_error("File import error",VOROPP_FILE_ERROR);
}
/** Outputs the a list of all the container regions along with the number of
* particles stored within each. */
void container_periodic_base::region_count() {
int i,j,k,*cop=co;
for(k=0;k<nz;k++) for(j=0;j<ny;j++) for(i=0;i<nx;i++)
printf("Region (%d,%d,%d): %d particles\n",i,j,k,*(cop++));
}
/** Clears a container of particles. */
void container_periodic::clear() {
for(int *cop=co;cop<co+nxyz;cop++) *cop=0;
}
/** Clears a container of particles, also clearing resetting the maximum radius
* to zero. */
void container_periodic_poly::clear() {
for(int *cop=co;cop<co+nxyz;cop++) *cop=0;
max_radius=0;
}
/** Computes all the Voronoi cells and saves customized information about them.
* \param[in] format the custom output string to use.
* \param[in] fp a file handle to write to. */
void container_periodic::print_custom(const char *format,FILE *fp) {
c_loop_all_periodic vl(*this);
print_custom(vl,format,fp);
}
/** Computes all the Voronoi cells and saves customized
* information about them.
* \param[in] format the custom output string to use.
* \param[in] fp a file handle to write to. */
void container_periodic_poly::print_custom(const char *format,FILE *fp) {
c_loop_all_periodic vl(*this);
print_custom(vl,format,fp);
}
/** Computes all the Voronoi cells and saves customized information about them.
* \param[in] format the custom output string to use.
* \param[in] filename the name of the file to write to. */
void container_periodic::print_custom(const char *format,const char *filename) {
FILE *fp=safe_fopen(filename,"w");
print_custom(format,fp);
fclose(fp);
}
/** Computes all the Voronoi cells and saves customized
* information about them
* \param[in] format the custom output string to use.
* \param[in] filename the name of the file to write to. */
void container_periodic_poly::print_custom(const char *format,const char *filename) {
FILE *fp=safe_fopen(filename,"w");
print_custom(format,fp);
fclose(fp);
}
/** Computes all of the Voronoi cells in the container, but does nothing
* with the output. It is useful for measuring the pure computation time
* of the Voronoi algorithm, without any additional calculations such as
* volume evaluation or cell output. */
void container_periodic::compute_all_cells() {
voronoicell c;
c_loop_all_periodic vl(*this);
if(vl.start()) do compute_cell(c,vl);
while(vl.inc());
}
/** Computes all of the Voronoi cells in the container, but does nothing
* with the output. It is useful for measuring the pure computation time
* of the Voronoi algorithm, without any additional calculations such as
* volume evaluation or cell output. */
void container_periodic_poly::compute_all_cells() {
voronoicell c;
c_loop_all_periodic vl(*this);
if(vl.start()) do compute_cell(c,vl);while(vl.inc());
}
/** Calculates all of the Voronoi cells and sums their volumes. In most cases
* without walls, the sum of the Voronoi cell volumes should equal the volume
* of the container to numerical precision.
* \return The sum of all of the computed Voronoi volumes. */
double container_periodic::sum_cell_volumes() {
voronoicell c;
double vol=0;
c_loop_all_periodic vl(*this);
if(vl.start()) do if(compute_cell(c,vl)) vol+=c.volume();while(vl.inc());
return vol;
}
/** Calculates all of the Voronoi cells and sums their volumes. In most cases
* without walls, the sum of the Voronoi cell volumes should equal the volume
* of the container to numerical precision.
* \return The sum of all of the computed Voronoi volumes. */
double container_periodic_poly::sum_cell_volumes() {
voronoicell c;
double vol=0;
c_loop_all_periodic vl(*this);
if(vl.start()) do if(compute_cell(c,vl)) vol+=c.volume();while(vl.inc());
return vol;
}
/** This routine creates all periodic images of the particles. It is meant for
* diagnostic purposes only, since usually periodic images are dynamically
* created in when they are referenced. */
void container_periodic_base::create_all_images() {
int i,j,k;
for(k=0;k<oz;k++) for(j=0;j<oy;j++) for(i=0;i<nx;i++) create_periodic_image(i,j,k);
}
/** Checks that the particles within each block lie within that block's bounds.
* This is useful for diagnosing problems with periodic image computation. */
void container_periodic_base::check_compartmentalized() {
int c,l,i,j,k;
double mix,miy,miz,max,may,maz,*pp;
for(k=l=0;k<oz;k++) for(j=0;j<oy;j++) for(i=0;i<nx;i++,l++) if(mem[l]>0) {
// Compute the block's bounds, adding in a small tolerance
mix=i*boxx-tolerance;max=mix+boxx+tolerance;
miy=(j-ey)*boxy-tolerance;may=miy+boxy+tolerance;
miz=(k-ez)*boxz-tolerance;maz=miz+boxz+tolerance;
// Print entries for any particles that lie outside the block's
// bounds
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)
printf("%d %d %d %d %f %f %f %f %f %f %f %f %f\n",
id[l][c],i,j,k,*pp,pp[1],pp[2],mix,max,miy,may,miz,maz);
}
}
/** Creates particles within an image block that is aligned with the primary
* domain in the z axis. In this case, the image block may be comprised of
* particles from two primary blocks. The routine considers these two primary
* blocks, and adds the needed particles to the image. The remaining particles
* from the primary blocks are also filled into the neighboring images.
* \param[in] (di,dj,dk) the index of the block to consider. The z index must
* satisfy ez<=dk<wz. */
void container_periodic_base::create_side_image(int di,int dj,int dk) {
int l,dijk=di+nx*(dj+oy*dk),odijk,ima=step_div(dj-ey,ny);
int qua=di+step_int(-ima*bxy*xsp),quadiv=step_div(qua,nx);
int fi=qua-quadiv*nx,fijk=fi+nx*(dj-ima*ny+oy*dk);
double dis=ima*bxy+quadiv*bx,switchx=di*boxx-ima*bxy-quadiv*bx,adis;
// Left image computation
if((img[dijk]&1)==0) {
if(di>0) {
odijk=dijk-1;adis=dis;
} else {
odijk=dijk+nx-1;adis=dis+bx;
}
img[odijk]|=2;
for(l=0;l<co[fijk];l++) {
if(p[fijk][ps*l]>switchx) put_image(dijk,fijk,l,dis,by*ima,0);
else put_image(odijk,fijk,l,adis,by*ima,0);
}
}
// Right image computation
if((img[dijk]&2)==0) {
if(fi==nx-1) {
fijk+=1-nx;switchx+=(1-nx)*boxx;dis+=bx;
} else {
fijk++;switchx+=boxx;
}
if(di==nx-1) {
odijk=dijk-nx+1;adis=dis-bx;
} else {
odijk=dijk+1;adis=dis;
}
img[odijk]|=1;
for(l=0;l<co[fijk];l++) {
if(p[fijk][ps*l]<switchx) put_image(dijk,fijk,l,dis,by*ima,0);
else put_image(odijk,fijk,l,adis,by*ima,0);
}
}
// All contributions to the block now added, so set both two bits of
// the image information
img[dijk]=3;
}
/** Creates particles within an image block that is not aligned with the
* primary domain in the z axis. In this case, the image block may be comprised
* of particles from four primary blocks. The routine considers these four
* primary blocks, and adds the needed particles to the image. The remaining
* particles from the primary blocks are also filled into the neighboring
* images.
* \param[in] (di,dj,dk) the index of the block to consider. The z index must
* satisfy dk<ez or dk>=wz. */
void container_periodic_base::create_vertical_image(int di,int dj,int dk) {
int l,dijk=di+nx*(dj+oy*dk),dijkl,dijkr,ima=step_div(dk-ez,nz);
int qj=dj+step_int(-ima*byz*ysp),qjdiv=step_div(qj-ey,ny);
int qi=di+step_int((-ima*bxz-qjdiv*bxy)*xsp),qidiv=step_div(qi,nx);
int fi=qi-qidiv*nx,fj=qj-qjdiv*ny,fijk=fi+nx*(fj+oy*(dk-ima*nz)),fijk2;
double disy=ima*byz+qjdiv*by,switchy=(dj-ey)*boxy-ima*byz-qjdiv*by;
double disx=ima*bxz+qjdiv*bxy+qidiv*bx,switchx=di*boxx-ima*bxz-qjdiv*bxy-qidiv*bx;
double switchx2,disxl,disxr,disx2,disxr2;
if(di==0) {dijkl=dijk+nx-1;disxl=disx+bx;}
else {dijkl=dijk-1;disxl=disx;}
if(di==nx-1) {dijkr=dijk-nx+1;disxr=disx-bx;}
else {dijkr=dijk+1;disxr=disx;}
// Down-left image computation
bool y_exist=dj!=0;
if((img[dijk]&1)==0) {
img[dijkl]|=2;
if(y_exist) {
img[dijkl-nx]|=8;
img[dijk-nx]|=4;
}
for(l=0;l<co[fijk];l++) {
if(p[fijk][ps*l+1]>switchy) {
if(p[fijk][ps*l]>switchx) put_image(dijk,fijk,l,disx,disy,bz*ima);
else put_image(dijkl,fijk,l,disxl,disy,bz*ima);
} else {
if(!y_exist) continue;
if(p[fijk][ps*l]>switchx) put_image(dijk-nx,fijk,l,disx,disy,bz*ima);
else put_image(dijkl-nx,fijk,l,disxl,disy,bz*ima);
}
}
}
// Down-right image computation
if((img[dijk]&2)==0) {
if(fi==nx-1) {
fijk2=fijk+1-nx;switchx2=switchx+(1-nx)*boxx;disx2=disx+bx;disxr2=disxr+bx;
} else {
fijk2=fijk+1;switchx2=switchx+boxx;disx2=disx;disxr2=disxr;
}
img[dijkr]|=1;
if(y_exist) {
img[dijkr-nx]|=4;
img[dijk-nx]|=8;
}
for(l=0;l<co[fijk2];l++) {
if(p[fijk2][ps*l+1]>switchy) {
if(p[fijk2][ps*l]>switchx2) put_image(dijkr,fijk2,l,disxr2,disy,bz*ima);
else put_image(dijk,fijk2,l,disx2,disy,bz*ima);
} else {
if(!y_exist) continue;
if(p[fijk2][ps*l]>switchx2) put_image(dijkr-nx,fijk2,l,disxr2,disy,bz*ima);
else put_image(dijk-nx,fijk2,l,disx2,disy,bz*ima);
}
}
}
// Recomputation of some intermediate quantities for boundary cases
if(fj==wy-1) {
fijk+=nx*(1-ny)-fi;
switchy+=(1-ny)*boxy;
disy+=by;
qi=di+step_int(-(ima*bxz+(qjdiv+1)*bxy)*xsp);
int dqidiv=step_div(qi,nx)-qidiv;qidiv+=dqidiv;
fi=qi-qidiv*nx;
fijk+=fi;
disx+=bxy+bx*dqidiv;
disxl+=bxy+bx*dqidiv;
disxr+=bxy+bx*dqidiv;
switchx-=bxy+bx*dqidiv;
} else {
fijk+=nx;switchy+=boxy;
}
// Up-left image computation
y_exist=dj!=oy-1;
if((img[dijk]&4)==0) {
img[dijkl]|=8;
if(y_exist) {
img[dijkl+nx]|=2;
img[dijk+nx]|=1;
}
for(l=0;l<co[fijk];l++) {
if(p[fijk][ps*l+1]>switchy) {
if(!y_exist) continue;
if(p[fijk][ps*l]>switchx) put_image(dijk+nx,fijk,l,disx,disy,bz*ima);
else put_image(dijkl+nx,fijk,l,disxl,disy,bz*ima);
} else {
if(p[fijk][ps*l]>switchx) put_image(dijk,fijk,l,disx,disy,bz*ima);
else put_image(dijkl,fijk,l,disxl,disy,bz*ima);
}
}
}
// Up-right image computation
if((img[dijk]&8)==0) {
if(fi==nx-1) {
fijk2=fijk+1-nx;switchx2=switchx+(1-nx)*boxx;disx2=disx+bx;disxr2=disxr+bx;
} else {
fijk2=fijk+1;switchx2=switchx+boxx;disx2=disx;disxr2=disxr;
}
img[dijkr]|=4;
if(y_exist) {
img[dijkr+nx]|=1;
img[dijk+nx]|=2;
}
for(l=0;l<co[fijk2];l++) {
if(p[fijk2][ps*l+1]>switchy) {
if(!y_exist) continue;
if(p[fijk2][ps*l]>switchx2) put_image(dijkr+nx,fijk2,l,disxr2,disy,bz*ima);
else put_image(dijk+nx,fijk2,l,disx2,disy,bz*ima);
} else {
if(p[fijk2][ps*l]>switchx2) put_image(dijkr,fijk2,l,disxr2,disy,bz*ima);
else put_image(dijk,fijk2,l,disx2,disy,bz*ima);
}
}
}
// All contributions to the block now added, so set all four bits of
// the image information
img[dijk]=15;
}
/** Copies a particle position from the primary domain into an image block.
* \param[in] reg the block index within the primary domain that the particle
* is within.
* \param[in] fijk the index of the image block.
* \param[in] l the index of the particle entry within the primary block.
* \param[in] (dx,dy,dz) the displacement vector to add to the particle. */
void container_periodic_base::put_image(int reg,int fijk,int l,double dx,double dy,double dz) {
if(co[reg]==mem[reg]) add_particle_memory(reg);
double *p1=p[reg]+ps*co[reg],*p2=p[fijk]+ps*l;
*(p1++)=*(p2++)+dx;
*(p1++)=*(p2++)+dy;
*p1=*p2+dz;
if(ps==4) *(++p1)=*(++p2);
id[reg][co[reg]++]=id[fijk][l];
}
}
c4science · Help