branches/dynamic/src/container.hh

   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     : July 1st 2008
   6
   7 /** \file container.hh
   8  * \brief Header file for the container_base template and related classes. */
   9
  10 #ifndef VOROPP_CONTAINER_HH
  11 #define VOROPP_CONTAINER_HH
  12
  13 #include "config.hh"
  14 #include <cstdio>
  15 #include <iostream>
  16 #include <fstream>
  17 #include <cmath>
  18 using namespace std;
  19
  20 class voropp_loop;
  21 class radius_mono;
  22 class radius_poly;
  23 class wall;
  24
  25 /** \brief A class representing the whole simulation region.
  26  *
  27  * The container class represents the whole simulation region. The
  28  * container constructor sets up the geometry and periodicity, and divides
  29  * the geometry into rectangular grid of blocks, each of which handles the
  30  * particles in a particular area. Routines exist for putting in particles,
  31  * importing particles from standard input, and carrying out Voronoi
  32  * calculations. */
  33 template<class r_option>
  34 class container_base {
  35         public:
  36                 container_base(fpoint xa,fpoint xb,fpoint ya,fpoint yb,fpoint za,fpoint zb,int xn,int yn,int zn,const bool xper,const bool yper,const bool zper,int memi);
  37                 virtual ~container_base();
  38                 void draw_particles(const char *filename);
  39                 void draw_particles();
  40                 void draw_particles(ostream &os);
  41                 void draw_particles_pov(const char *filename);
  42                 void draw_particles_pov();
  43                 void draw_particles_pov(ostream &os);
  44                 void draw_lammps_restart(ostream &os,fpoint timestep,bool scaled=false);
  45                 void import(istream &is);
  46                 inline void import();
  47                 inline void import(const char *filename);
  48                 void skip_lammps_restart(istream &is);
  49                 void import_lammps_restart(istream &is,bool scaled=false);
  50                 void region_count();
  51                 void clear();
  52                 void draw_cells_gnuplot(const char *filename,fpoint xmin,fpoint xmax,fpoint ymin,fpoint ymax,fpoint zmin,fpoint zmax);
  53                 inline void draw_cells_gnuplot(const char *filename);
  54                 void draw_cells_pov(const char *filename,fpoint xmin,fpoint xmax,fpoint ymin,fpoint ymax,fpoint zmin,fpoint zmax);
  55                 inline void draw_cells_pov(const char *filename);
  56                 void store_cell_volumes(fpoint *bb);
  57                 fpoint packing_fraction(fpoint *bb,fpoint cx,fpoint cy,fpoint cz,fpoint r);
  58                 fpoint packing_fraction(fpoint *bb,fpoint xmin,fpoint xmax,fpoint ymin,fpoint ymax,fpoint zmin,fpoint zmax);
  59                 fpoint sum_cell_volumes();
  60                 void compute_all_cells();
  61                 void count_all_faces(ostream &os);
  62                 void count_all_faces();
  63                 void count_all_faces(const char *filename);
  64                 void print_all(ostream &os);
  65                 void print_all();
  66                 void print_all(const char *filename);
  67                 void print_all_neighbor(ostream &os);
  68                 void print_all_neighbor();
  69                 void print_all_neighbor(const char *filename);
  70                 template<class n_option>
  71                 inline bool compute_cell_sphere(voronoicell_base<n_option> &c,int i,int j,int k,int ijk,int s);
  72                 template<class n_option>
  73                 bool compute_cell_sphere(voronoicell_base<n_option> &c,int i,int j,int k,int ijk,int s,fpoint x,fpoint y,fpoint z);
  74                 template<class n_option>
  75                 inline bool compute_cell(voronoicell_base<n_option> &c,int i,int j,int k,int ijk,int s);
  76                 template<class n_option>
  77                 bool compute_cell(voronoicell_base<n_option> &c,int i,int j,int k,int ijk,int s,fpoint x,fpoint y,fpoint z);
  78                 void put(int n,fpoint x,fpoint y,fpoint z);
  79                 void put(int n,fpoint x,fpoint y,fpoint z,fpoint r);
  80                 void add_wall(wall &w);
  81                 bool point_inside(fpoint x,fpoint y,fpoint z);
  82                 bool point_inside_walls(fpoint x,fpoint y,fpoint z);
  83         protected:
  84                 /** The minimum x coordinate of the container. */
  85                 const fpoint ax;
  86                 /** The maximum x coordinate of the container. */
  87                 const fpoint bx;
  88                 /** The minimum y coordinate of the container. */
  89                 const fpoint ay;
  90                 /** The maximum y coordinate of the container. */
  91                 const fpoint by;
  92                 /** The minimum z coordinate of the container. */
  93                 const fpoint az;
  94                 /** The maximum z coordinate of the container. */
  95                 const fpoint bz;
  96                 /** The inverse box length in the x direction, set to
  97                  * nx/(bx-ax). */
  98                 const fpoint xsp;
  99                 /** The inverse box length in the y direction, set to
 100                  * ny/(by-ay). */
 101                 const fpoint ysp;
 102                 /** The inverse box length in the z direction, set to
 103                  * nz/(bz-az). */
 104                 const fpoint zsp;
 105                 /** The number of boxes in the x direction. */
 106                 const int nx;
 107                 /** The number of boxes in the y direction. */
 108                 const int ny;
 109                 /** The number of boxes in the z direction. */
 110                 const int nz;
 111                 /** A constant, set to the value of nx multiplied by ny, which
 112                  * is used in the routines which step through boxes in
 113                  * sequence. */
 114                 const int nxy;
 115                 /** A constant, set to the value of nx*ny*nz, which is used in
 116                  * the routines which step through boxes in sequence. */
 117                 const int nxyz;
 118                 /** The number of boxes in the x direction for the searching mask. */
 119                 const int hx;
 120                 /** The number of boxes in the y direction for the searching mask. */
 121                 const int hy;
 122                 /** The number of boxes in the z direction for the searching mask. */
 123                 const int hz;
 124                 /** A constant, set to the value of hx multiplied by hy, which
 125                  * is used in the routines which step through mask boxes in
 126                  * sequence. */
 127                 const int hxy;
 128                 /** A constant, set to the value of hx*hy*hz, which is used in
 129                  * the routines which step through mask boxes in sequence. */
 130                 const int hxyz;
 131                 /** A boolean value that determines if the x coordinate in
 132                  * periodic or not. */
 133                 const bool xperiodic;
 134                 /** A boolean value that determines if the y coordinate in
 135                  * periodic or not. */
 136                 const bool yperiodic;
 137                 /** A boolean value that determines if the z coordinate in
 138                  * periodic or not. */
 139                 const bool zperiodic;
 140                 /** This array holds the number of particles within each
 141                  * computational box of the container. */
 142                 int *co;
 143                 /** This array holds the maximum amount of particle memory for
 144                  * each computational box of the container. If the number of
 145                  * particles in a particular box ever approaches this limit,
 146                  * more is allocated using the add_particle_memory() function.
 147                  */
 148                 int *mem;
 149                 /** This array holds the numerical IDs of each particle in each
 150                  * computational box. */
 151                 int **id;
 152                 /** This array is used as a mask. */
 153                 unsigned int *mask;
 154                 /** This array is used to store the list of blocks to test during
 155                  * the Voronoi cell computation. */
 156                 int *sl;
 157                 /** This sets the current value being used to mark tested blocks
 158                  * in the mask. */
 159                 unsigned int mv;
 160                 /** The position of the first element on the search list to be
 161                  * considered. */
 162                 int s_start;
 163                 /** The position of the last element on the search list to be
 164                  * considered. */
 165                 int s_end;
 166                 /** The current size of the search list. */
 167                 int s_size;
 168                 /** A two dimensional array holding particle positions. For the
 169                  * derived container_poly class, this also holds particle
 170                  * radii. */
 171                 fpoint **p;
 172                 /** This array holds pointers to any wall objects that have
 173                  * been added to the container. */
 174                 wall **walls;
 175                 /** The current number of wall objects, initially set to zero. */
 176                 int wall_number;
 177                 /** The current amount of memory allocated for walls. */
 178                 int current_wall_size;
 179                 /** This object contains all the functions for handling how
 180                  * the particle radii should be treated. If the template is
 181                  * instantiated with the radius_mono class, then this object
 182                  * contains mostly blank routines that do nothing to the
 183                  * cell computation, to compute the basic Voronoi diagram.
 184                  * If the template is instantiated with the radius_poly calls,
 185                  * then this object provides routines for modifying the
 186                  * Voronoi cell computation in order to create the radical
 187                  * Voronoi tessellation. */
 188                 r_option radius;
 189                 /** The amount of memory in the array structure for each
 190                  * particle. This is set to 3 when the basic class is
 191                  * initialized, so that the array holds (x,y,z) positions. If
 192                  * the container class is initialized as part of the derived
 193                  * class container_poly, then this is set to 4, to also hold
 194                  * the particle radii. */
 195                 int sz;
 196                 /** An array to hold the minimum distances associated with the
 197                  * worklists. This array is initialized during container
 198                  * construction, by the initialize_radii() routine. */
 199                 fpoint *mrad;
 200
 201                 template<class n_option>
 202                 inline void print_all(ostream &os,voronoicell_base<n_option> &c);
 203                 template<class n_option>
 204                 inline bool initialize_voronoicell(voronoicell_base<n_option> &c,fpoint x,fpoint y,fpoint z);
 205                 virtual void add_particle_memory(int i);
 206                 void add_list_memory();
 207         private:
 208 #include "worklist.hh"
 209                 template<class n_option>
 210                 inline bool corner_test(voronoicell_base<n_option> &c,fpoint xl,fpoint yl,fpoint zl,fpoint xh,fpoint yh,fpoint zh);
 211                 template<class n_option>
 212                 inline bool edge_x_test(voronoicell_base<n_option> &c,fpoint x0,fpoint yl,fpoint zl,fpoint x1,fpoint yh,fpoint zh);
 213                 template<class n_option>
 214                 inline bool edge_y_test(voronoicell_base<n_option> &c,fpoint xl,fpoint y0,fpoint zl,fpoint xh,fpoint y1,fpoint zh);
 215                 template<class n_option>
 216                 inline bool edge_z_test(voronoicell_base<n_option> &c,fpoint xl,fpoint yl,fpoint z0,fpoint xh,fpoint yh,fpoint z1);
 217                 template<class n_option>
 218                 inline bool face_x_test(voronoicell_base<n_option> &c,fpoint xl,fpoint y0,fpoint z0,fpoint y1,fpoint z1);
 219                 template<class n_option>
 220                 inline bool face_y_test(voronoicell_base<n_option> &c,fpoint x0,fpoint yl,fpoint z0,fpoint x1,fpoint z1);
 221                 template<class n_option>
 222                 inline bool face_z_test(voronoicell_base<n_option> &c,fpoint x0,fpoint y0,fpoint zl,fpoint x1,fpoint y1);
 223                 inline void initialize_radii();
 224                 inline void compute_minimum(fpoint &minr,fpoint &xlo,fpoint &xhi,fpoint &ylo,fpoint &yhi,fpoint &zlo,fpoint &zhi,int ti,int tj,int tk);
 225                 inline bool compute_min_max_radius(int di,int dj,int dk,fpoint fx,fpoint fy,fpoint fz,fpoint gx,fpoint gy,fpoint gz,fpoint& crs,fpoint mrs);
 226                 friend class voropp_loop;
 227                 friend class radius_poly;
 228                 friend class radius_mono;
 229 };
 230
 231 /** \brief A class encapsulating all routines specifically needed in
 232  * the standard Voronoi tessellation.
 233  *
 234  * This class encapsulates all the routines that are required for carrying out
 235  * a standard Voronoi tessellation that would be appropriate for a monodisperse
 236  * system. When the container class is instantiated using this class, all
 237  * information about particle radii is switched off. Since all these functions
 238  * are declared inline, there should be no loss of speed. */
 239 class radius_mono {
 240         public:
 241                 /** The number of floating point numbers allocated for each
 242                  * particle in the container, set to 3 for this case for the x,
 243                  * y, and z positions. */
 244                 const int mem_size;
 245                 /** This constructor sets a pointer back to the container class
 246                  * that created it, and initializes the mem_size constant to 3.
 247                  */
 248                 radius_mono(container_base<radius_mono> *icc) : mem_size(3), cc(icc) {};
 249                 inline void import(istream &is);
 250                 /** This is a blank placeholder function that does nothing. */
 251                 inline void store_radius(int i,int j,fpoint r) {};
 252                 /** This is a blank placeholder function that does nothing. */
 253                 inline void clear_max() {};
 254                 /** This is a blank placeholder function that does nothing. */
 255                 inline void init(int s,int i) {};
 256                 inline fpoint volume(int ijk,int s);
 257                 inline fpoint cutoff(fpoint lrs);
 258                 inline fpoint scale(fpoint rs,int t,int q);
 259                 /** This is a blank placeholder function that does nothing. */
 260                 inline void print(ostream &os,int ijk,int q) {};
 261                 inline void rad(ostream &os,int l,int c);
 262                 inline void diam(ostream &os,int l,int c);
 263                 inline void import_lammps(istream &is,int n,bool scaled);
 264         private:
 265                 container_base<radius_mono> *cc;
 266 };
 267
 268 /** \brief A class encapsulating all routines specifically needed in the
 269  * Voronoi radical tessellation.
 270  *
 271  * This class encapsulates all the routines that are required for carrying out
 272  * the radical Voronoi tessellation that is appropriate for polydisperse sphere.
 273  * When the container class is instantiated with this class, information about particle
 274  * radii is switched on. */
 275 class radius_poly {
 276         public:
 277                 /** The number of floating point numbers allocated for each
 278                  * particle in the container, set to 4 for this case for the x,
 279                  * y, and z positions, plus the radius. */
 280                 const int mem_size;
 281                 /** This constructor sets a pointer back to the container class
 282                  * that created it, and initializes the mem_size constant to 4.
 283                  */
 284                 radius_poly(container_base<radius_poly> *icc) : mem_size(4), cc(icc) {};
 285                 inline void import(istream &is);
 286                 inline void store_radius(int i,int j,fpoint r);
 287                 inline void clear_max();
 288                 inline void init(int ijk,int s);
 289                 inline fpoint volume(int ijk,int s);
 290                 inline fpoint cutoff(fpoint lrs);
 291                 inline fpoint scale(fpoint rs,int t,int q);
 292                 inline void print(ostream &os,int ijk,int q);
 293                 inline void rad(ostream &os,int l,int c);
 294                 inline void diam(ostream &os,int l,int c);
 295                 inline void import_lammps(istream &is,int n,bool scaled);
 296         private:
 297                 container_base<radius_poly> *cc;
 298                 fpoint max_radius,crad,mul;
 299 };
 300
 301 /** \brief A class to handle loops on regions of the container handling
 302  * non-periodic and periodic boundary conditions.
 303  *
 304  * Many of the container routines require scanning over a rectangular sub-grid
 305  * of blocks, and the routines for handling this are stored in the voropp_loop
 306  * class. A voropp_loop class can first be initialized to either calculate the
 307  * subgrid which is within a distance r of a vector (vx,vy,vz), or a subgrid
 308  * corresponding to a rectangular box. The routine inc() can then be
 309  * successively called to step through all the blocks within this subgrid.
 310  */
 311 class voropp_loop {
 312         public:
 313                 template<class r_option>
 314                 voropp_loop(container_base<r_option> *q);
 315                 inline int init(fpoint vx,fpoint vy,fpoint vz,fpoint r,fpoint &px,fpoint &py,fpoint &pz);
 316                 inline int init(fpoint xmin,fpoint xmax,fpoint ymin,fpoint ymax,fpoint zmin,fpoint zmax,fpoint &px,fpoint &py,fpoint &pz);
 317                 inline int inc(fpoint &px,fpoint &py,fpoint &pz);
 318                 inline int reset(fpoint &px,fpoint &py,fpoint &pz);
 319                 inline bool reached(int ri,int rj,int rk);
 320                 /** The current block index in the x direction, referencing a
 321                  * real cell in the range 0 to nx-1. */
 322                 int ip;
 323                 /** The current block index in the y direction, referencing a
 324                  * real cell in the range 0 to ny-1. */
 325                 int jp;
 326                 /** The current block index in the z direction, referencing a
 327                  * real cell in the range 0 to nz-1. */
 328                 int kp;
 329                 /** The current block index in the x direction, possibly refering
 330                  * to an image block outside the normal range of 0 to nx-1. */
 331                 int i;
 332                 /** The current block index in the y direction, possibly refering
 333                  * to an image block outside the normal range of 0 to ny-1. */
 334                 int j;
 335                 /** The current block index in the z direction, possibly refering
 336                  * to an image block outside the normal range of 0 to nz-1. */
 337                 int k;
 338         private:
 339                 int ai,bi,aj,bj,ak,bk,s;
 340                 int aip,ajp,akp,inc1,inc2;
 341                 inline int step_mod(int a,int b);
 342                 inline int step_div(int a,int b);
 343                 inline int step_int(fpoint a);
 344                 fpoint apx,apy,apz;
 345                 const fpoint sx,sy,sz,xsp,ysp,zsp,ax,ay,az;
 346                 const int nx,ny,nz,nxy,nxyz;
 347                 const bool xperiodic,yperiodic,zperiodic;
 348 };
 349
 350 /** \brief Pure virtual class from which wall objects are derived.
 351  *
 352  * This is a pure virtual class for a generic wall object. A wall object
 353  * can be specified by deriving a new class from this and specifying the
 354  * functions.*/
 355 class wall {
 356         public:
 357                 virtual ~wall() {};
 358                 /** A pure virtual function for testing whether a point is
 359                  * inside the wall object. */
 360                 virtual bool point_inside(fpoint x,fpoint y,fpoint z) = 0;
 361                 /** A pure virtual function for cutting a cell without
 362                  * neighbor-tracking with a wall. */
 363                 virtual bool cut_cell(voronoicell_base<neighbor_none> &c,fpoint x,fpoint y,fpoint z) = 0;
 364                 /** A pure virtual function for cutting a cell with
 365                  * neighbor-tracking enabled with a wall. */
 366                 virtual bool cut_cell(voronoicell_base<neighbor_track> &c,fpoint x,fpoint y,fpoint z) = 0;
 367                 /** A pure virtual function for returning the minimum distance
 368                  * vector from a position to the wall, used in the dynamic
 369                  * relaxation calculations. */
 370                 virtual void min_distance(fpoint x,fpoint y,fpoint z,fpoint &dx,fpoint &dy,fpoint &dz) = 0;
 371 };
 372
 373 /** The basic container class. */
 374 typedef container_base<radius_mono> container;
 375
 376 /** The polydisperse container class. */
 377 typedef container_base<radius_poly> container_poly;
 378 #endif