| 1 |
francois |
979 |
// 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 container.cc |
| 8 |
|
|
* \brief Function implementations for the container and related classes. */ |
| 9 |
|
|
|
| 10 |
|
|
#include "container.hh" |
| 11 |
|
|
|
| 12 |
|
|
namespace voro { |
| 13 |
|
|
|
| 14 |
|
|
/** The class constructor sets up the geometry of container, initializing the |
| 15 |
|
|
* minimum and maximum coordinates in each direction, and setting whether each |
| 16 |
|
|
* direction is periodic or not. It divides the container into a rectangular |
| 17 |
|
|
* grid of blocks, and allocates memory for each of these for storing particle |
| 18 |
|
|
* positions and IDs. |
| 19 |
|
|
* \param[in] (ax_,bx_) the minimum and maximum x coordinates. |
| 20 |
|
|
* \param[in] (ay_,by_) the minimum and maximum y coordinates. |
| 21 |
|
|
* \param[in] (az_,bz_) the minimum and maximum z coordinates. |
| 22 |
|
|
* \param[in] (nx_,ny_,nz_) the number of grid blocks in each of the three |
| 23 |
|
|
* coordinate directions. |
| 24 |
|
|
* \param[in] (xperiodic_,yperiodic_,zperiodic_) flags setting whether the |
| 25 |
|
|
* container is periodic in each |
| 26 |
|
|
* coordinate direction. |
| 27 |
|
|
* \param[in] init_mem the initial memory allocation for each block. |
| 28 |
|
|
* \param[in] ps_ the number of floating point entries to store for each |
| 29 |
|
|
* particle. */ |
| 30 |
|
|
container_base::container_base(double ax_,double bx_,double ay_,double by_,double az_,double bz_, |
| 31 |
|
|
int nx_,int ny_,int nz_,bool xperiodic_,bool yperiodic_,bool zperiodic_,int init_mem,int ps_) |
| 32 |
|
|
: voro_base(nx_,ny_,nz_,(bx_-ax_)/nx_,(by_-ay_)/ny_,(bz_-az_)/nz_), |
| 33 |
|
|
ax(ax_), bx(bx_), ay(ay_), by(by_), az(az_), bz(bz_), |
| 34 |
|
|
xperiodic(xperiodic_), yperiodic(yperiodic_), zperiodic(zperiodic_), |
| 35 |
|
|
id(new int*[nxyz]), p(new double*[nxyz]), co(new int[nxyz]), mem(new int[nxyz]), ps(ps_) { |
| 36 |
|
|
int l; |
| 37 |
|
|
for(l=0;l<nxyz;l++) co[l]=0; |
| 38 |
|
|
for(l=0;l<nxyz;l++) mem[l]=init_mem; |
| 39 |
|
|
for(l=0;l<nxyz;l++) id[l]=new int[init_mem]; |
| 40 |
|
|
for(l=0;l<nxyz;l++) p[l]=new double[ps*init_mem]; |
| 41 |
|
|
} |
| 42 |
|
|
|
| 43 |
|
|
/** The container destructor frees the dynamically allocated memory. */ |
| 44 |
|
|
container_base::~container_base() { |
| 45 |
|
|
int l; |
| 46 |
|
|
for(l=0;l<nxyz;l++) delete [] p[l]; |
| 47 |
|
|
for(l=0;l<nxyz;l++) delete [] id[l]; |
| 48 |
|
|
delete [] id; |
| 49 |
|
|
delete [] p; |
| 50 |
|
|
delete [] co; |
| 51 |
|
|
delete [] mem; |
| 52 |
|
|
} |
| 53 |
|
|
|
| 54 |
|
|
/** The class constructor sets up the geometry of container. |
| 55 |
|
|
* \param[in] (ax_,bx_) the minimum and maximum x coordinates. |
| 56 |
|
|
* \param[in] (ay_,by_) the minimum and maximum y coordinates. |
| 57 |
|
|
* \param[in] (az_,bz_) the minimum and maximum z coordinates. |
| 58 |
|
|
* \param[in] (nx_,ny_,nz_) the number of grid blocks in each of the three |
| 59 |
|
|
* coordinate directions. |
| 60 |
|
|
* \param[in] (xperiodic_,yperiodic_,zperiodic_) flags setting whether the |
| 61 |
|
|
* container is periodic in each |
| 62 |
|
|
* coordinate direction. |
| 63 |
|
|
* \param[in] init_mem the initial memory allocation for each block. */ |
| 64 |
|
|
container::container(double ax_,double bx_,double ay_,double by_,double az_,double bz_, |
| 65 |
|
|
int nx_,int ny_,int nz_,bool xperiodic_,bool yperiodic_,bool zperiodic_,int init_mem) |
| 66 |
|
|
: container_base(ax_,bx_,ay_,by_,az_,bz_,nx_,ny_,nz_,xperiodic_,yperiodic_,zperiodic_,init_mem,3), |
| 67 |
|
|
vc(*this,xperiodic_?2*nx_+1:nx_,yperiodic_?2*ny_+1:ny_,zperiodic_?2*nz_+1:nz_) {} |
| 68 |
|
|
|
| 69 |
|
|
/** The class constructor sets up the geometry of container. |
| 70 |
|
|
* \param[in] (ax_,bx_) the minimum and maximum x coordinates. |
| 71 |
|
|
* \param[in] (ay_,by_) the minimum and maximum y coordinates. |
| 72 |
|
|
* \param[in] (az_,bz_) the minimum and maximum z coordinates. |
| 73 |
|
|
* \param[in] (nx_,ny_,nz_) the number of grid blocks in each of the three |
| 74 |
|
|
* coordinate directions. |
| 75 |
|
|
* \param[in] (xperiodic_,yperiodic_,zperiodic_) flags setting whether the |
| 76 |
|
|
* container is periodic in each |
| 77 |
|
|
* coordinate direction. |
| 78 |
|
|
* \param[in] init_mem the initial memory allocation for each block. */ |
| 79 |
|
|
container_poly::container_poly(double ax_,double bx_,double ay_,double by_,double az_,double bz_, |
| 80 |
|
|
int nx_,int ny_,int nz_,bool xperiodic_,bool yperiodic_,bool zperiodic_,int init_mem) |
| 81 |
|
|
: container_base(ax_,bx_,ay_,by_,az_,bz_,nx_,ny_,nz_,xperiodic_,yperiodic_,zperiodic_,init_mem,4), |
| 82 |
|
|
vc(*this,xperiodic_?2*nx_+1:nx_,yperiodic_?2*ny_+1:ny_,zperiodic_?2*nz_+1:nz_) {ppr=p;} |
| 83 |
|
|
|
| 84 |
|
|
/** Put a particle into the correct region of the container. |
| 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 |
|
|
void container::put(int n,double x,double y,double z) { |
| 88 |
|
|
int ijk; |
| 89 |
|
|
if(put_locate_block(ijk,x,y,z)) { |
| 90 |
|
|
id[ijk][co[ijk]]=n; |
| 91 |
|
|
double *pp=p[ijk]+3*co[ijk]++; |
| 92 |
|
|
*(pp++)=x;*(pp++)=y;*pp=z; |
| 93 |
|
|
} |
| 94 |
|
|
} |
| 95 |
|
|
|
| 96 |
|
|
/** Put a particle into the correct region of the container. |
| 97 |
|
|
* \param[in] n the numerical ID of the inserted particle. |
| 98 |
|
|
* \param[in] (x,y,z) the position vector of the inserted particle. |
| 99 |
|
|
* \param[in] r the radius of the particle. */ |
| 100 |
|
|
void container_poly::put(int n,double x,double y,double z,double r) { |
| 101 |
|
|
int ijk; |
| 102 |
|
|
if(put_locate_block(ijk,x,y,z)) { |
| 103 |
|
|
id[ijk][co[ijk]]=n; |
| 104 |
|
|
double *pp=p[ijk]+4*co[ijk]++; |
| 105 |
|
|
*(pp++)=x;*(pp++)=y;*(pp++)=z;*pp=r; |
| 106 |
|
|
if(max_radius<r) max_radius=r; |
| 107 |
|
|
} |
| 108 |
|
|
} |
| 109 |
|
|
|
| 110 |
|
|
/** Put a particle into the correct region of the container, also recording |
| 111 |
|
|
* into which region it was stored. |
| 112 |
|
|
* \param[in] vo the ordering class in which to record the region. |
| 113 |
|
|
* \param[in] n the numerical ID of the inserted particle. |
| 114 |
|
|
* \param[in] (x,y,z) the position vector of the inserted particle. */ |
| 115 |
|
|
void container::put(particle_order &vo,int n,double x,double y,double z) { |
| 116 |
|
|
int ijk; |
| 117 |
|
|
if(put_locate_block(ijk,x,y,z)) { |
| 118 |
|
|
id[ijk][co[ijk]]=n; |
| 119 |
|
|
vo.add(ijk,co[ijk]); |
| 120 |
|
|
double *pp=p[ijk]+3*co[ijk]++; |
| 121 |
|
|
*(pp++)=x;*(pp++)=y;*pp=z; |
| 122 |
|
|
} |
| 123 |
|
|
} |
| 124 |
|
|
|
| 125 |
|
|
/** Put a particle into the correct region of the container, also recording |
| 126 |
|
|
* into which region it was stored. |
| 127 |
|
|
* \param[in] vo the ordering class in which to record the region. |
| 128 |
|
|
* \param[in] n the numerical ID of the inserted particle. |
| 129 |
|
|
* \param[in] (x,y,z) the position vector of the inserted particle. |
| 130 |
|
|
* \param[in] r the radius of the particle. */ |
| 131 |
|
|
void container_poly::put(particle_order &vo,int n,double x,double y,double z,double r) { |
| 132 |
|
|
int ijk; |
| 133 |
|
|
if(put_locate_block(ijk,x,y,z)) { |
| 134 |
|
|
id[ijk][co[ijk]]=n; |
| 135 |
|
|
vo.add(ijk,co[ijk]); |
| 136 |
|
|
double *pp=p[ijk]+4*co[ijk]++; |
| 137 |
|
|
*(pp++)=x;*(pp++)=y;*(pp++)=z;*pp=r; |
| 138 |
|
|
if(max_radius<r) max_radius=r; |
| 139 |
|
|
} |
| 140 |
|
|
} |
| 141 |
|
|
|
| 142 |
|
|
/** This routine takes a particle position vector, tries to remap it into the |
| 143 |
|
|
* primary domain. If successful, it computes the region into which it can be |
| 144 |
|
|
* stored and checks that there is enough memory within this region to store |
| 145 |
|
|
* it. |
| 146 |
|
|
* \param[out] ijk the region index. |
| 147 |
|
|
* \param[in,out] (x,y,z) the particle position, remapped into the primary |
| 148 |
|
|
* domain if necessary. |
| 149 |
|
|
* \return True if the particle can be successfully placed into the container, |
| 150 |
|
|
* false otherwise. */ |
| 151 |
|
|
bool container_base::put_locate_block(int &ijk,double &x,double &y,double &z) { |
| 152 |
|
|
if(put_remap(ijk,x,y,z)) { |
| 153 |
|
|
if(co[ijk]==mem[ijk]) add_particle_memory(ijk); |
| 154 |
|
|
return true; |
| 155 |
|
|
} |
| 156 |
|
|
#if VOROPP_REPORT_OUT_OF_BOUNDS ==1 |
| 157 |
|
|
fprintf(stderr,"Out of bounds: (x,y,z)=(%g,%g,%g)\n",x,y,z); |
| 158 |
|
|
#endif |
| 159 |
|
|
return false; |
| 160 |
|
|
} |
| 161 |
|
|
|
| 162 |
|
|
/** Takes a particle position vector and computes the region index into which |
| 163 |
|
|
* it should be stored. If the container is periodic, then the routine also |
| 164 |
|
|
* maps the particle position to ensure it is in the primary domain. If the |
| 165 |
|
|
* container is not periodic, the routine bails out. |
| 166 |
|
|
* \param[out] ijk the region index. |
| 167 |
|
|
* \param[in,out] (x,y,z) the particle position, remapped into the primary |
| 168 |
|
|
* domain if necessary. |
| 169 |
|
|
* \return True if the particle can be successfully placed into the container, |
| 170 |
|
|
* false otherwise. */ |
| 171 |
|
|
inline bool container_base::put_remap(int &ijk,double &x,double &y,double &z) { |
| 172 |
|
|
int l; |
| 173 |
|
|
|
| 174 |
|
|
ijk=step_int((x-ax)*xsp); |
| 175 |
|
|
if(xperiodic) {l=step_mod(ijk,nx);x+=boxx*(l-ijk);ijk=l;} |
| 176 |
|
|
else if(ijk<0||ijk>=nx) return false; |
| 177 |
|
|
|
| 178 |
|
|
int j=step_int((y-ay)*ysp); |
| 179 |
|
|
if(yperiodic) {l=step_mod(j,ny);y+=boxy*(l-j);j=l;} |
| 180 |
|
|
else if(j<0||j>=ny) return false; |
| 181 |
|
|
|
| 182 |
|
|
int k=step_int((z-az)*zsp); |
| 183 |
|
|
if(zperiodic) {l=step_mod(k,nz);z+=boxz*(l-k);k=l;} |
| 184 |
|
|
else if(k<0||k>=nz) return false; |
| 185 |
|
|
|
| 186 |
|
|
ijk+=nx*j+nxy*k; |
| 187 |
|
|
return true; |
| 188 |
|
|
} |
| 189 |
|
|
|
| 190 |
|
|
/** Takes a position vector and attempts to remap it into the primary domain. |
| 191 |
|
|
* \param[out] (ai,aj,ak) the periodic image displacement that the vector is in, |
| 192 |
|
|
* with (0,0,0) corresponding to the primary domain. |
| 193 |
|
|
* \param[out] (ci,cj,ck) the index of the block that the position vector is |
| 194 |
|
|
* within, once it has been remapped. |
| 195 |
|
|
* \param[in,out] (x,y,z) the position vector to consider, which is remapped |
| 196 |
|
|
* into the primary domain during the routine. |
| 197 |
|
|
* \param[out] ijk the block index that the vector is within. |
| 198 |
|
|
* \return True if the particle is within the container or can be remapped into |
| 199 |
|
|
* it, false if it lies outside of the container bounds. */ |
| 200 |
|
|
inline bool container_base::remap(int &ai,int &aj,int &ak,int &ci,int &cj,int &ck,double &x,double &y,double &z,int &ijk) { |
| 201 |
|
|
ci=step_int((x-ax)*xsp); |
| 202 |
|
|
if(ci<0||ci>=nx) { |
| 203 |
|
|
if(xperiodic) {ai=step_div(ci,nx);x-=ai*(bx-ax);ci-=ai*nx;} |
| 204 |
|
|
else return false; |
| 205 |
|
|
} else ai=0; |
| 206 |
|
|
|
| 207 |
|
|
cj=step_int((y-ay)*ysp); |
| 208 |
|
|
if(cj<0||cj>=ny) { |
| 209 |
|
|
if(yperiodic) {aj=step_div(cj,ny);y-=aj*(by-ay);cj-=aj*ny;} |
| 210 |
|
|
else return false; |
| 211 |
|
|
} else aj=0; |
| 212 |
|
|
|
| 213 |
|
|
ck=step_int((z-az)*zsp); |
| 214 |
|
|
if(ck<0||ck>=nz) { |
| 215 |
|
|
if(zperiodic) {ak=step_div(ck,nz);z-=ak*(bz-az);ck-=ak*nz;} |
| 216 |
|
|
else return false; |
| 217 |
|
|
} else ak=0; |
| 218 |
|
|
|
| 219 |
|
|
ijk=ci+nx*cj+nxy*ck; |
| 220 |
|
|
return true; |
| 221 |
|
|
} |
| 222 |
|
|
|
| 223 |
|
|
/** Takes a vector and finds the particle whose Voronoi cell contains that |
| 224 |
|
|
* vector. This is equivalent to finding the particle which is nearest to the |
| 225 |
|
|
* vector. Additional wall classes are not considered by this routine. |
| 226 |
|
|
* \param[in] (x,y,z) the vector to test. |
| 227 |
|
|
* \param[out] (rx,ry,rz) the position of the particle whose Voronoi cell |
| 228 |
|
|
* contains the vector. If the container is periodic, |
| 229 |
|
|
* this may point to a particle in a periodic image of |
| 230 |
|
|
* the primary domain. |
| 231 |
|
|
* \param[out] pid the ID of the particle. |
| 232 |
|
|
* \return True if a particle was found. If the container has no particles, |
| 233 |
|
|
* then the search will not find a Voronoi cell and false is returned. */ |
| 234 |
|
|
bool container::find_voronoi_cell(double x,double y,double z,double &rx,double &ry,double &rz,int &pid) { |
| 235 |
|
|
int ai,aj,ak,ci,cj,ck,ijk; |
| 236 |
|
|
particle_record w; |
| 237 |
|
|
double mrs; |
| 238 |
|
|
|
| 239 |
|
|
// If the given vector lies outside the domain, but the container |
| 240 |
|
|
// is periodic, then remap it back into the domain |
| 241 |
|
|
if(!remap(ai,aj,ak,ci,cj,ck,x,y,z,ijk)) return false; |
| 242 |
|
|
vc.find_voronoi_cell(x,y,z,ci,cj,ck,ijk,w,mrs); |
| 243 |
|
|
|
| 244 |
|
|
if(w.ijk!=-1) { |
| 245 |
|
|
|
| 246 |
|
|
// Assemble the position vector of the particle to be returned, |
| 247 |
|
|
// applying a periodic remapping if necessary |
| 248 |
|
|
if(xperiodic) {ci+=w.di;if(ci<0||ci>=nx) ai+=step_div(ci,nx);} |
| 249 |
|
|
if(yperiodic) {cj+=w.dj;if(cj<0||cj>=ny) aj+=step_div(cj,ny);} |
| 250 |
|
|
if(zperiodic) {ck+=w.dk;if(ck<0||ck>=nz) ak+=step_div(ck,nz);} |
| 251 |
|
|
rx=p[w.ijk][3*w.l]+ai*(bx-ax); |
| 252 |
|
|
ry=p[w.ijk][3*w.l+1]+aj*(by-ay); |
| 253 |
|
|
rz=p[w.ijk][3*w.l+2]+ak*(bz-az); |
| 254 |
|
|
pid=id[w.ijk][w.l]; |
| 255 |
|
|
return true; |
| 256 |
|
|
} |
| 257 |
|
|
|
| 258 |
|
|
// If no particle if found then just return false |
| 259 |
|
|
return false; |
| 260 |
|
|
} |
| 261 |
|
|
|
| 262 |
|
|
/** Takes a vector and finds the particle whose Voronoi cell contains that |
| 263 |
|
|
* vector. Additional wall classes are not considered by this routine. |
| 264 |
|
|
* \param[in] (x,y,z) the vector to test. |
| 265 |
|
|
* \param[out] (rx,ry,rz) the position of the particle whose Voronoi cell |
| 266 |
|
|
* contains the vector. If the container is periodic, |
| 267 |
|
|
* this may point to a particle in a periodic image of |
| 268 |
|
|
* the primary domain. |
| 269 |
|
|
* \param[out] pid the ID of the particle. |
| 270 |
|
|
* \return True if a particle was found. If the container has no particles, |
| 271 |
|
|
* then the search will not find a Voronoi cell and false is returned. */ |
| 272 |
|
|
bool container_poly::find_voronoi_cell(double x,double y,double z,double &rx,double &ry,double &rz,int &pid) { |
| 273 |
|
|
int ai,aj,ak,ci,cj,ck,ijk; |
| 274 |
|
|
particle_record w; |
| 275 |
|
|
double mrs; |
| 276 |
|
|
|
| 277 |
|
|
// If the given vector lies outside the domain, but the container |
| 278 |
|
|
// is periodic, then remap it back into the domain |
| 279 |
|
|
if(!remap(ai,aj,ak,ci,cj,ck,x,y,z,ijk)) return false; |
| 280 |
|
|
vc.find_voronoi_cell(x,y,z,ci,cj,ck,ijk,w,mrs); |
| 281 |
|
|
|
| 282 |
|
|
if(w.ijk!=-1) { |
| 283 |
|
|
|
| 284 |
|
|
// Assemble the position vector of the particle to be returned, |
| 285 |
|
|
// applying a periodic remapping if necessary |
| 286 |
|
|
if(xperiodic) {ci+=w.di;if(ci<0||ci>=nx) ai+=step_div(ci,nx);} |
| 287 |
|
|
if(yperiodic) {cj+=w.dj;if(cj<0||cj>=ny) aj+=step_div(cj,ny);} |
| 288 |
|
|
if(zperiodic) {ck+=w.dk;if(ck<0||ck>=nz) ak+=step_div(ck,nz);} |
| 289 |
|
|
rx=p[w.ijk][4*w.l]+ai*(bx-ax); |
| 290 |
|
|
ry=p[w.ijk][4*w.l+1]+aj*(by-ay); |
| 291 |
|
|
rz=p[w.ijk][4*w.l+2]+ak*(bz-az); |
| 292 |
|
|
pid=id[w.ijk][w.l]; |
| 293 |
|
|
return true; |
| 294 |
|
|
} |
| 295 |
|
|
|
| 296 |
|
|
// If no particle if found then just return false |
| 297 |
|
|
return false; |
| 298 |
|
|
} |
| 299 |
|
|
|
| 300 |
|
|
/** Increase memory for a particular region. |
| 301 |
|
|
* \param[in] i the index of the region to reallocate. */ |
| 302 |
|
|
void container_base::add_particle_memory(int i) { |
| 303 |
|
|
int l,nmem=mem[i]<<1; |
| 304 |
|
|
|
| 305 |
|
|
// Carry out a check on the memory allocation size, and |
| 306 |
|
|
// print a status message if requested |
| 307 |
|
|
if(nmem>max_particle_memory) |
| 308 |
|
|
voro_fatal_error("Absolute maximum memory allocation exceeded",VOROPP_MEMORY_ERROR); |
| 309 |
|
|
#if VOROPP_VERBOSE >=3 |
| 310 |
|
|
fprintf(stderr,"Particle memory in region %d scaled up to %d\n",i,nmem); |
| 311 |
|
|
#endif |
| 312 |
|
|
|
| 313 |
|
|
// Allocate new memory and copy in the contents of the old arrays |
| 314 |
|
|
int *idp=new int[nmem]; |
| 315 |
|
|
for(l=0;l<co[i];l++) idp[l]=id[i][l]; |
| 316 |
|
|
double *pp=new double[ps*nmem]; |
| 317 |
|
|
for(l=0;l<ps*co[i];l++) pp[l]=p[i][l]; |
| 318 |
|
|
|
| 319 |
|
|
// Update pointers and delete old arrays |
| 320 |
|
|
mem[i]=nmem; |
| 321 |
|
|
delete [] id[i];id[i]=idp; |
| 322 |
|
|
delete [] p[i];p[i]=pp; |
| 323 |
|
|
} |
| 324 |
|
|
|
| 325 |
|
|
/** Import a list of particles from an open file stream into the container. |
| 326 |
|
|
* Entries of four numbers (Particle ID, x position, y position, z position) |
| 327 |
|
|
* are searched for. If the file cannot be successfully read, then the routine |
| 328 |
|
|
* causes a fatal error. |
| 329 |
|
|
* \param[in] fp the file handle to read from. */ |
| 330 |
|
|
void container::import(FILE *fp) { |
| 331 |
|
|
int i,j; |
| 332 |
|
|
double x,y,z; |
| 333 |
|
|
while((j=fscanf(fp,"%d %lg %lg %lg",&i,&x,&y,&z))==4) put(i,x,y,z); |
| 334 |
|
|
if(j!=EOF) voro_fatal_error("File import error",VOROPP_FILE_ERROR); |
| 335 |
|
|
} |
| 336 |
|
|
|
| 337 |
|
|
/** Import a list of particles from an open file stream, also storing the order |
| 338 |
|
|
* of that the particles are read. Entries of four numbers (Particle ID, x |
| 339 |
|
|
* position, y position, z position) are searched for. If the file cannot be |
| 340 |
|
|
* successfully read, then the routine causes a fatal error. |
| 341 |
|
|
* \param[in,out] vo a reference to an ordering class to use. |
| 342 |
|
|
* \param[in] fp the file handle to read from. */ |
| 343 |
|
|
void container::import(particle_order &vo,FILE *fp) { |
| 344 |
|
|
int i,j; |
| 345 |
|
|
double x,y,z; |
| 346 |
|
|
while((j=fscanf(fp,"%d %lg %lg %lg",&i,&x,&y,&z))==4) put(vo,i,x,y,z); |
| 347 |
|
|
if(j!=EOF) voro_fatal_error("File import error",VOROPP_FILE_ERROR); |
| 348 |
|
|
} |
| 349 |
|
|
|
| 350 |
|
|
/** Import a list of particles from an open file stream into the container. |
| 351 |
|
|
* Entries of five numbers (Particle ID, x position, y position, z position, |
| 352 |
|
|
* radius) are searched for. If the file cannot be successfully read, then the |
| 353 |
|
|
* routine causes a fatal error. |
| 354 |
|
|
* \param[in] fp the file handle to read from. */ |
| 355 |
|
|
void container_poly::import(FILE *fp) { |
| 356 |
|
|
int i,j; |
| 357 |
|
|
double x,y,z,r; |
| 358 |
|
|
while((j=fscanf(fp,"%d %lg %lg %lg %lg",&i,&x,&y,&z,&r))==5) put(i,x,y,z,r); |
| 359 |
|
|
if(j!=EOF) voro_fatal_error("File import error",VOROPP_FILE_ERROR); |
| 360 |
|
|
} |
| 361 |
|
|
|
| 362 |
|
|
/** Import a list of particles from an open file stream, also storing the order |
| 363 |
|
|
* of that the particles are read. Entries of four numbers (Particle ID, x |
| 364 |
|
|
* position, y position, z position, radius) are searched for. If the file |
| 365 |
|
|
* cannot be successfully read, then the routine causes a fatal error. |
| 366 |
|
|
* \param[in,out] vo a reference to an ordering class to use. |
| 367 |
|
|
* \param[in] fp the file handle to read from. */ |
| 368 |
|
|
void container_poly::import(particle_order &vo,FILE *fp) { |
| 369 |
|
|
int i,j; |
| 370 |
|
|
double x,y,z,r; |
| 371 |
|
|
while((j=fscanf(fp,"%d %lg %lg %lg %lg",&i,&x,&y,&z,&r))==5) put(vo,i,x,y,z,r); |
| 372 |
|
|
if(j!=EOF) voro_fatal_error("File import error",VOROPP_FILE_ERROR); |
| 373 |
|
|
} |
| 374 |
|
|
|
| 375 |
|
|
/** Outputs the a list of all the container regions along with the number of |
| 376 |
|
|
* particles stored within each. */ |
| 377 |
|
|
void container_base::region_count() { |
| 378 |
|
|
int i,j,k,*cop=co; |
| 379 |
|
|
for(k=0;k<nz;k++) for(j=0;j<ny;j++) for(i=0;i<nx;i++) |
| 380 |
|
|
printf("Region (%d,%d,%d): %d particles\n",i,j,k,*(cop++)); |
| 381 |
|
|
} |
| 382 |
|
|
|
| 383 |
|
|
/** Clears a container of particles. */ |
| 384 |
|
|
void container::clear() { |
| 385 |
|
|
for(int *cop=co;cop<co+nxyz;cop++) *cop=0; |
| 386 |
|
|
} |
| 387 |
|
|
|
| 388 |
|
|
/** Clears a container of particles, also clearing resetting the maximum radius |
| 389 |
|
|
* to zero. */ |
| 390 |
|
|
void container_poly::clear() { |
| 391 |
|
|
for(int *cop=co;cop<co+nxyz;cop++) *cop=0; |
| 392 |
|
|
max_radius=0; |
| 393 |
|
|
} |
| 394 |
|
|
|
| 395 |
|
|
/** Computes all the Voronoi cells and saves customized information about them. |
| 396 |
|
|
* \param[in] format the custom output string to use. |
| 397 |
|
|
* \param[in] fp a file handle to write to. */ |
| 398 |
|
|
void container::print_custom(const char *format,FILE *fp) { |
| 399 |
|
|
c_loop_all vl(*this); |
| 400 |
|
|
print_custom(vl,format,fp); |
| 401 |
|
|
} |
| 402 |
|
|
|
| 403 |
|
|
/** Computes all the Voronoi cells and saves customized |
| 404 |
|
|
* information about them. |
| 405 |
|
|
* \param[in] format the custom output string to use. |
| 406 |
|
|
* \param[in] fp a file handle to write to. */ |
| 407 |
|
|
void container_poly::print_custom(const char *format,FILE *fp) { |
| 408 |
|
|
c_loop_all vl(*this); |
| 409 |
|
|
print_custom(vl,format,fp); |
| 410 |
|
|
} |
| 411 |
|
|
|
| 412 |
|
|
/** Computes all the Voronoi cells and saves customized information about them. |
| 413 |
|
|
* \param[in] format the custom output string to use. |
| 414 |
|
|
* \param[in] filename the name of the file to write to. */ |
| 415 |
|
|
void container::print_custom(const char *format,const char *filename) { |
| 416 |
|
|
FILE *fp=safe_fopen(filename,"w"); |
| 417 |
|
|
print_custom(format,fp); |
| 418 |
|
|
fclose(fp); |
| 419 |
|
|
} |
| 420 |
|
|
|
| 421 |
|
|
/** Computes all the Voronoi cells and saves customized |
| 422 |
|
|
* information about them |
| 423 |
|
|
* \param[in] format the custom output string to use. |
| 424 |
|
|
* \param[in] filename the name of the file to write to. */ |
| 425 |
|
|
void container_poly::print_custom(const char *format,const char *filename) { |
| 426 |
|
|
FILE *fp=safe_fopen(filename,"w"); |
| 427 |
|
|
print_custom(format,fp); |
| 428 |
|
|
fclose(fp); |
| 429 |
|
|
} |
| 430 |
|
|
|
| 431 |
|
|
/** Computes all of the Voronoi cells in the container, but does nothing |
| 432 |
|
|
* with the output. It is useful for measuring the pure computation time |
| 433 |
|
|
* of the Voronoi algorithm, without any additional calculations such as |
| 434 |
|
|
* volume evaluation or cell output. */ |
| 435 |
|
|
void container::compute_all_cells() { |
| 436 |
|
|
voronoicell c; |
| 437 |
|
|
c_loop_all vl(*this); |
| 438 |
|
|
if(vl.start()) do compute_cell(c,vl); |
| 439 |
|
|
while(vl.inc()); |
| 440 |
|
|
} |
| 441 |
|
|
|
| 442 |
|
|
/** Computes all of the Voronoi cells in the container, but does nothing |
| 443 |
|
|
* with the output. It is useful for measuring the pure computation time |
| 444 |
|
|
* of the Voronoi algorithm, without any additional calculations such as |
| 445 |
|
|
* volume evaluation or cell output. */ |
| 446 |
|
|
void container_poly::compute_all_cells() { |
| 447 |
|
|
voronoicell c; |
| 448 |
|
|
c_loop_all vl(*this); |
| 449 |
|
|
if(vl.start()) do compute_cell(c,vl);while(vl.inc()); |
| 450 |
|
|
} |
| 451 |
|
|
|
| 452 |
|
|
/** Calculates all of the Voronoi cells and sums their volumes. In most cases |
| 453 |
|
|
* without walls, the sum of the Voronoi cell volumes should equal the volume |
| 454 |
|
|
* of the container to numerical precision. |
| 455 |
|
|
* \return The sum of all of the computed Voronoi volumes. */ |
| 456 |
|
|
double container::sum_cell_volumes() { |
| 457 |
|
|
voronoicell c; |
| 458 |
|
|
double vol=0; |
| 459 |
|
|
c_loop_all vl(*this); |
| 460 |
|
|
if(vl.start()) do if(compute_cell(c,vl)) vol+=c.volume();while(vl.inc()); |
| 461 |
|
|
return vol; |
| 462 |
|
|
} |
| 463 |
|
|
|
| 464 |
|
|
/** Calculates all of the Voronoi cells and sums their volumes. In most cases |
| 465 |
|
|
* without walls, the sum of the Voronoi cell volumes should equal the volume |
| 466 |
|
|
* of the container to numerical precision. |
| 467 |
|
|
* \return The sum of all of the computed Voronoi volumes. */ |
| 468 |
|
|
double container_poly::sum_cell_volumes() { |
| 469 |
|
|
voronoicell c; |
| 470 |
|
|
double vol=0; |
| 471 |
|
|
c_loop_all vl(*this); |
| 472 |
|
|
if(vl.start()) do if(compute_cell(c,vl)) vol+=c.volume();while(vl.inc()); |
| 473 |
|
|
return vol; |
| 474 |
|
|
} |
| 475 |
|
|
|
| 476 |
|
|
/** This function tests to see if a given vector lies within the container |
| 477 |
|
|
* bounds and any walls. |
| 478 |
|
|
* \param[in] (x,y,z) the position vector to be tested. |
| 479 |
|
|
* \return True if the point is inside the container, false if the point is |
| 480 |
|
|
* outside. */ |
| 481 |
|
|
bool container_base::point_inside(double x,double y,double z) { |
| 482 |
|
|
if(x<ax||x>bx||y<ay||y>by||z<az||z>bz) return false; |
| 483 |
|
|
return point_inside_walls(x,y,z); |
| 484 |
|
|
} |
| 485 |
|
|
|
| 486 |
|
|
/** Draws an outline of the domain in gnuplot format. |
| 487 |
|
|
* \param[in] fp the file handle to write to. */ |
| 488 |
|
|
void container_base::draw_domain_gnuplot(FILE *fp) { |
| 489 |
|
|
fprintf(fp,"%g %g %g\n%g %g %g\n%g %g %g\n%g %g %g\n",ax,ay,az,bx,ay,az,bx,by,az,ax,by,az); |
| 490 |
|
|
fprintf(fp,"%g %g %g\n%g %g %g\n%g %g %g\n%g %g %g\n",ax,by,bz,bx,by,bz,bx,ay,bz,ax,ay,bz); |
| 491 |
|
|
fprintf(fp,"%g %g %g\n\n%g %g %g\n%g %g %g\n\n",ax,by,bz,ax,ay,az,ax,ay,bz); |
| 492 |
|
|
fprintf(fp,"%g %g %g\n%g %g %g\n\n%g %g %g\n%g %g %g\n\n",bx,ay,az,bx,ay,bz,bx,by,az,bx,by,bz); |
| 493 |
|
|
} |
| 494 |
|
|
|
| 495 |
|
|
/** Draws an outline of the domain in POV-Ray format. |
| 496 |
|
|
* \param[in] fp the file handle to write to. */ |
| 497 |
|
|
void container_base::draw_domain_pov(FILE *fp) { |
| 498 |
|
|
fprintf(fp,"cylinder{<%g,%g,%g>,<%g,%g,%g>,rr}\n" |
| 499 |
|
|
"cylinder{<%g,%g,%g>,<%g,%g,%g>,rr}\n",ax,ay,az,bx,ay,az,ax,by,az,bx,by,az); |
| 500 |
|
|
fprintf(fp,"cylinder{<%g,%g,%g>,<%g,%g,%g>,rr}\n" |
| 501 |
|
|
"cylinder{<%g,%g,%g>,<%g,%g,%g>,rr}\n",ax,by,bz,bx,by,bz,ax,ay,bz,bx,ay,bz); |
| 502 |
|
|
fprintf(fp,"cylinder{<%g,%g,%g>,<%g,%g,%g>,rr}\n" |
| 503 |
|
|
"cylinder{<%g,%g,%g>,<%g,%g,%g>,rr}\n",ax,ay,az,ax,by,az,bx,ay,az,bx,by,az); |
| 504 |
|
|
fprintf(fp,"cylinder{<%g,%g,%g>,<%g,%g,%g>,rr}\n" |
| 505 |
|
|
"cylinder{<%g,%g,%g>,<%g,%g,%g>,rr}\n",bx,ay,bz,bx,by,bz,ax,ay,bz,ax,by,bz); |
| 506 |
|
|
fprintf(fp,"cylinder{<%g,%g,%g>,<%g,%g,%g>,rr}\n" |
| 507 |
|
|
"cylinder{<%g,%g,%g>,<%g,%g,%g>,rr}\n",ax,ay,az,ax,ay,bz,bx,ay,az,bx,ay,bz); |
| 508 |
|
|
fprintf(fp,"cylinder{<%g,%g,%g>,<%g,%g,%g>,rr}\n" |
| 509 |
|
|
"cylinder{<%g,%g,%g>,<%g,%g,%g>,rr}\n",bx,by,az,bx,by,bz,ax,by,az,ax,by,bz); |
| 510 |
|
|
fprintf(fp,"sphere{<%g,%g,%g>,rr}\nsphere{<%g,%g,%g>,rr}\n" |
| 511 |
|
|
"sphere{<%g,%g,%g>,rr}\nsphere{<%g,%g,%g>,rr}\n",ax,ay,az,bx,ay,az,ax,by,az,bx,by,az); |
| 512 |
|
|
fprintf(fp,"sphere{<%g,%g,%g>,rr}\nsphere{<%g,%g,%g>,rr}\n" |
| 513 |
|
|
"sphere{<%g,%g,%g>,rr}\nsphere{<%g,%g,%g>,rr}\n",ax,ay,bz,bx,ay,bz,ax,by,bz,bx,by,bz); |
| 514 |
|
|
} |
| 515 |
|
|
|
| 516 |
|
|
|
| 517 |
|
|
/** The wall_list constructor sets up an array of pointers to wall classes. */ |
| 518 |
|
|
wall_list::wall_list() : walls(new wall*[init_wall_size]), wep(walls), wel(walls+init_wall_size), |
| 519 |
|
|
current_wall_size(init_wall_size) {} |
| 520 |
|
|
|
| 521 |
|
|
/** The wall_list destructor frees the array of pointers to the wall classes. |
| 522 |
|
|
*/ |
| 523 |
|
|
wall_list::~wall_list() { |
| 524 |
|
|
delete [] walls; |
| 525 |
|
|
} |
| 526 |
|
|
|
| 527 |
|
|
/** Adds all of the walls on another wall_list to this class. |
| 528 |
|
|
* \param[in] wl a reference to the wall class. */ |
| 529 |
|
|
void wall_list::add_wall(wall_list &wl) { |
| 530 |
|
|
for(wall **wp=wl.walls;wp<wl.wep;wp++) add_wall(*wp); |
| 531 |
|
|
} |
| 532 |
|
|
|
| 533 |
|
|
/** Deallocates all of the wall classes pointed to by the wall_list. */ |
| 534 |
|
|
void wall_list::deallocate() { |
| 535 |
|
|
for(wall **wp=walls;wp<wep;wp++) delete *wp; |
| 536 |
|
|
} |
| 537 |
|
|
|
| 538 |
|
|
/** Increases the memory allocation for the walls array. */ |
| 539 |
|
|
void wall_list::increase_wall_memory() { |
| 540 |
|
|
current_wall_size<<=1; |
| 541 |
|
|
if(current_wall_size>max_wall_size) |
| 542 |
|
|
voro_fatal_error("Wall memory allocation exceeded absolute maximum",VOROPP_MEMORY_ERROR); |
| 543 |
|
|
wall **nwalls=new wall*[current_wall_size],**nwp=nwalls,**wp=walls; |
| 544 |
|
|
while(wp<wep) *(nwp++)=*(wp++); |
| 545 |
|
|
delete [] walls; |
| 546 |
|
|
walls=nwalls;wel=walls+current_wall_size;wep=nwp; |
| 547 |
|
|
} |
| 548 |
|
|
|
| 549 |
|
|
} |
