1 // Voro++, a 3D cell-based Voronoi library 2 // 3 // Author : Chris H. Rycroft (LBL / UC Berkeley) 4 // Email : chr@alum.mit.edu 5 // Date : August 30th 2011 6 7 /** \file pre_container.cc 8 * \brief Function implementations for the pre_container and related classes. 9 */ 10 11 #include <cmath> 12 13 #include "config.hh" 14 #include "pre_container.hh" 15 16 namespace voro { 17 18 /** The class constructor sets up the geometry of container, initializing the 19 * minimum and maximum coordinates in each direction. It allocates an initial 20 * chunk into which to store particle information. 21 * \param[in] (ax_,bx_) the minimum and maximum x coordinates. 22 * \param[in] (ay_,by_) the minimum and maximum y coordinates. 23 * \param[in] (az_,bz_) the minimum and maximum z coordinates. 24 * \param[in] (xperiodic_,yperiodic_,zperiodic_ ) flags setting whether the 25 * container is periodic in each 26 * coordinate direction. 27 * \param[in] ps_ the number of floating point entries to store for each 28 * particle. */ 29 pre_container_base::pre_container_base(double ax_,double bx_,double ay_,double by_,double az_,double bz_, 30 bool xperiodic_,bool yperiodic_,bool zperiodic_,int ps_) : 31 ax(ax_), bx(bx_), ay(ay_), by(by_), az(az_), bz(bz_), 32 xperiodic(xperiodic_), yperiodic(yperiodic_), zperiodic(zperiodic_), ps(ps_), 33 index_sz(init_chunk_size), pre_id(new int*[index_sz]), end_id(pre_id), 34 pre_p(new double*[index_sz]), end_p(pre_p) { 35 ch_id=*end_id=new int[pre_container_chunk_size]; 36 l_id=end_id+index_sz;e_id=ch_id+pre_container_chunk_size; 37 ch_p=*end_p=new double[ps*pre_container_chunk_size]; 38 } 39 40 /** The destructor frees the dynamically allocated memory. */ 41 pre_container_base::~pre_container_base() { 42 delete [] *end_p; 43 delete [] *end_id; 44 while (end_id!=pre_id) { 45 end_p--; 46 delete [] *end_p; 47 end_id--; 48 delete [] *end_id; 49 } 50 delete [] pre_p; 51 delete [] pre_id; 52 } 53 54 /** Makes a guess at the optimal grid of blocks to use, computing in 55 * a way that 56 * \param[out] (nx,ny,nz) the number of blocks to use. */ 57 void pre_container_base::guess_optimal(int &nx,int &ny,int &nz) { 58 double dx=bx-ax,dy=by-ay,dz=bz-az; 59 double ilscale=pow(total_particles()/(optimal_particles*dx*dy*dz),1/3.0); 60 nx=int(dx*ilscale+1); 61 ny=int(dy*ilscale+1); 62 nz=int(dz*ilscale+1); 63 } 64 65 /** Stores a particle ID and position, allocating a new memory chunk if 66 * necessary. For coordinate directions in which the container is not periodic, 67 * the routine checks to make sure that the particle is within the container 68 * bounds. If the particle is out of bounds, it is not stored. 69 * \param[in] n the numerical ID of the inserted particle. 70 * \param[in] (x,y,z) the position vector of the inserted particle. */ 71 void pre_container::put(int n,double x,double y,double z) { 72 if((xperiodic||(x>=ax&&x<=bx))&&(yperiodic||(y>=ay&&y<=by))&&(zperiodic||(z>=az&&z<=bz))) { 73 if(ch_id==e_id) new_chunk(); 74 *(ch_id++)=n; 75 *(ch_p++)=x;*(ch_p++)=y;*(ch_p++)=z; 76 } 77 #if VOROPP_REPORT_OUT_OF_BOUNDS ==1 78 else fprintf(stderr,"Out of bounds: (x,y,z)=(%g,%g,%g)\n",x,y,z); 79 #endif 80 } 81 82 /** Stores a particle ID and position, allocating a new memory chunk if necessary. 83 * \param[in] n the numerical ID of the inserted particle. 84 * \param[in] (x,y,z) the position vector of the inserted particle. 85 * \param[in] r the radius of the particle. */ 86 void pre_container_poly::put(int n,double x,double y,double z,double r) { 87 if((xperiodic||(x>=ax&&x<=bx))&&(yperiodic||(y>=ay&&y<=by))&&(zperiodic||(z>=az&&z<=bz))) { 88 if(ch_id==e_id) new_chunk(); 89 *(ch_id++)=n; 90 *(ch_p++)=x;*(ch_p++)=y;*(ch_p++)=z;*(ch_p++)=r; 91 } 92 #if VOROPP_REPORT_OUT_OF_BOUNDS ==1 93 else fprintf(stderr,"Out of bounds: (x,y,z)=(%g,%g,%g)\n",x,y,z); 94 #endif 95 } 96 97 /** Transfers the particles stored within the class to a container class. 98 * \param[in] con the container class to transfer to. */ 99 void pre_container::setup(container &con) { 100 int **c_id=pre_id,*idp,*ide,n; 101 double **c_p=pre_p,*pp,x,y,z; 102 while(c_id<end_id) { 103 idp=*(c_id++);ide=idp+pre_container_chunk_size; 104 pp=*(c_p++); 105 while(idp<ide) { 106 n=*(idp++);x=*(pp++);y=*(pp++);z=*(pp++); 107 con.put(n,x,y,z); 108 } 109 } 110 idp=*c_id; 111 pp=*c_p; 112 while(idp<ch_id) { 113 n=*(idp++);x=*(pp++);y=*(pp++);z=*(pp++); 114 con.put(n,x,y,z); 115 } 116 } 117 118 /** Transfers the particles stored within the class to a container_poly class. 119 * \param[in] con the container_poly class to transfer to. */ 120 void pre_container_poly::setup(container_poly &con) { 121 int **c_id=pre_id,*idp,*ide,n; 122 double **c_p=pre_p,*pp,x,y,z,r; 123 while(c_id<end_id) { 124 idp=*(c_id++);ide=idp+pre_container_chunk_size; 125 pp=*(c_p++); 126 while(idp<ide) { 127 n=*(idp++);x=*(pp++);y=*(pp++);z=*(pp++);r=*(pp++); 128 con.put(n,x,y,z,r); 129 } 130 } 131 idp=*c_id; 132 pp=*c_p; 133 while(idp<ch_id) { 134 n=*(idp++);x=*(pp++);y=*(pp++);z=*(pp++);r=*(pp++); 135 con.put(n,x,y,z,r); 136 } 137 } 138 139 /** Transfers the particles stored within the class to a container class, also 140 * recording the order in which particles were stored. 141 * \param[in] vo the ordering class to use. 142 * \param[in] con the container class to transfer to. */ 143 void pre_container::setup(particle_order &vo,container &con) { 144 int **c_id=pre_id,*idp,*ide,n; 145 double **c_p=pre_p,*pp,x,y,z; 146 while(c_id<end_id) { 147 idp=*(c_id++);ide=idp+pre_container_chunk_size; 148 pp=*(c_p++); 149 while(idp<ide) { 150 n=*(idp++);x=*(pp++);y=*(pp++);z=*(pp++); 151 con.put(vo,n,x,y,z); 152 } 153 } 154 idp=*c_id; 155 pp=*c_p; 156 while(idp<ch_id) { 157 n=*(idp++);x=*(pp++);y=*(pp++);z=*(pp++); 158 con.put(vo,n,x,y,z); 159 } 160 } 161 162 /** Transfers the particles stored within the class to a container_poly class, 163 * also recording the order in which particles were stored. 164 * \param[in] vo the ordering class to use. 165 * \param[in] con the container_poly class to transfer to. */ 166 void pre_container_poly::setup(particle_order &vo,container_poly &con) { 167 int **c_id=pre_id,*idp,*ide,n; 168 double **c_p=pre_p,*pp,x,y,z,r; 169 while(c_id<end_id) { 170 idp=*(c_id++);ide=idp+pre_container_chunk_size; 171 pp=*(c_p++); 172 while(idp<ide) { 173 n=*(idp++);x=*(pp++);y=*(pp++);z=*(pp++);r=*(pp++); 174 con.put(vo,n,x,y,z,r); 175 } 176 } 177 idp=*c_id; 178 pp=*c_p; 179 while(idp<ch_id) { 180 n=*(idp++);x=*(pp++);y=*(pp++);z=*(pp++);r=*(pp++); 181 con.put(vo,n,x,y,z,r); 182 } 183 } 184 185 /** Import a list of particles from an open file stream into the container. 186 * Entries of four numbers (Particle ID, x position, y position, z position) 187 * are searched for. If the file cannot be successfully read, then the routine 188 * causes a fatal error. 189 * \param[in] fp the file handle to read from. */ 190 void pre_container::import(FILE *fp) { 191 int i,j; 192 double x,y,z; 193 while((j=fscanf(fp,"%d %lg %lg %lg",&i,&x,&y,&z))==4) put(i,x,y,z); 194 if(j!=EOF) voro_fatal_error("File import error",VOROPP_FILE_ERROR); 195 } 196 197 /** Import a list of particles from an open file stream, also storing the order 198 * of that the particles are read. Entries of four numbers (Particle ID, x 199 * position, y position, z position) are searched for. If the file cannot be 200 * successfully read, then the routine causes a fatal error. 201 * \param[in] fp the file handle to read from. */ 202 void pre_container_poly::import(FILE *fp) { 203 int i,j; 204 double x,y,z,r; 205 while((j=fscanf(fp,"%d %lg %lg %lg %lg",&i,&x,&y,&z,&r))==5) put(i,x,y,z,r); 206 if(j!=EOF) voro_fatal_error("File import error",VOROPP_FILE_ERROR); 207 } 208 209 /** Allocates a new chunk of memory for storing particles. */ 210 void pre_container_base::new_chunk() { 211 end_id++;end_p++; 212 if(end_id==l_id) extend_chunk_index(); 213 ch_id=*end_id=new int[pre_container_chunk_size]; 214 e_id=ch_id+pre_container_chunk_size; 215 ch_p=*end_p=new double[ps*pre_container_chunk_size]; 216 } 217 218 /** Extends the index of chunks. */ 219 void pre_container_base::extend_chunk_index() { 220 index_sz<<=1; 221 if(index_sz>max_chunk_size) 222 voro_fatal_error("Absolute memory limit on chunk index reached",VOROPP_MEMORY_ERROR); 223 #if VOROPP_VERBOSE >=2 224 fprintf(stderr,"Pre-container chunk index scaled up to %d\n",index_sz); 225 #endif 226 int **n_id=new int*[index_sz],**p_id=n_id,**c_id=pre_id; 227 double **n_p=new double*[index_sz],**p_p=n_p,**c_p=pre_p; 228 while(c_id<end_id) { 229 *(p_id++)=*(c_id++); 230 *(p_p++)=*(c_p++); 231 } 232 delete [] pre_id;pre_id=n_id;end_id=p_id;l_id=pre_id+index_sz; 233 delete [] pre_p;pre_p=n_p;end_p=p_p; 234 } 235 236 } 237