trunk/src/pre_container.cc

   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 }