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