| 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 unitcell.cc |
| 8 |
|
|
* \brief Function implementations for the unitcell class. */ |
| 9 |
|
|
|
| 10 |
|
|
#include <cmath> |
| 11 |
|
|
#include <queue> |
| 12 |
|
|
|
| 13 |
|
|
#include "unitcell.hh" |
| 14 |
|
|
#include "cell.hh" |
| 15 |
|
|
|
| 16 |
|
|
namespace voro { |
| 17 |
|
|
|
| 18 |
|
|
/** Initializes the unit cell class for a particular non-orthogonal periodic |
| 19 |
|
|
* geometry, corresponding to a parallelepiped with sides given by three |
| 20 |
|
|
* vectors. The class constructs the unit Voronoi cell corresponding to this |
| 21 |
|
|
* geometry. |
| 22 |
|
|
* \param[in] (bx_) The x coordinate of the first unit vector. |
| 23 |
|
|
* \param[in] (bxy_,by_) The x and y coordinates of the second unit vector. |
| 24 |
|
|
* \param[in] (bxz_,byz_,bz_) The x, y, and z coordinates of the third unit |
| 25 |
|
|
* vector. */ |
| 26 |
|
|
unitcell::unitcell(double bx_,double bxy_,double by_,double bxz_,double byz_,double bz_) |
| 27 |
|
|
: bx(bx_), bxy(bxy_), by(by_), bxz(bxz_), byz(byz_), bz(bz_) { |
| 28 |
|
|
int i,j,l=1; |
| 29 |
|
|
|
| 30 |
|
|
// Initialize the Voronoi cell to be a very large rectangular box |
| 31 |
|
|
const double ucx=max_unit_voro_shells*bx,ucy=max_unit_voro_shells*by,ucz=max_unit_voro_shells*bz; |
| 32 |
|
|
unit_voro.init(-ucx,ucx,-ucy,ucy,-ucz,ucz); |
| 33 |
|
|
|
| 34 |
|
|
// Repeatedly cut the cell by shells of periodic image particles |
| 35 |
|
|
while(l<2*max_unit_voro_shells) { |
| 36 |
|
|
|
| 37 |
|
|
// Check to see if any of the planes from the current shell |
| 38 |
|
|
// will cut the cell |
| 39 |
|
|
if(unit_voro_intersect(l)) { |
| 40 |
|
|
|
| 41 |
|
|
// If they do, apply the plane cuts from the current |
| 42 |
|
|
// shell |
| 43 |
|
|
unit_voro_apply(l,0,0); |
| 44 |
|
|
for(i=1;i<l;i++) { |
| 45 |
|
|
unit_voro_apply(l,i,0); |
| 46 |
|
|
unit_voro_apply(-l,i,0); |
| 47 |
|
|
} |
| 48 |
|
|
for(i=-l;i<=l;i++) unit_voro_apply(i,l,0); |
| 49 |
|
|
for(i=1;i<l;i++) for(j=-l+1;j<=l;j++) { |
| 50 |
|
|
unit_voro_apply(l,j,i); |
| 51 |
|
|
unit_voro_apply(-j,l,i); |
| 52 |
|
|
unit_voro_apply(-l,-j,i); |
| 53 |
|
|
unit_voro_apply(j,-l,i); |
| 54 |
|
|
} |
| 55 |
|
|
for(i=-l;i<=l;i++) for(j=-l;j<=l;j++) unit_voro_apply(i,j,l); |
| 56 |
|
|
} else { |
| 57 |
|
|
|
| 58 |
|
|
// Calculate a bound on the maximum y and z coordinates |
| 59 |
|
|
// that could possibly cut the cell. This is based upon |
| 60 |
|
|
// a geometric result that particles with z>l can't cut |
| 61 |
|
|
// a cell lying within the paraboloid |
| 62 |
|
|
// z<=(l*l-x*x-y*y)/(2*l). It is always a tighter bound |
| 63 |
|
|
// than the one based on computing the maximum radius |
| 64 |
|
|
// of a Voronoi cell vertex. |
| 65 |
|
|
max_uv_y=max_uv_z=0; |
| 66 |
|
|
double y,z,q,*pts=unit_voro.pts,*pp=pts; |
| 67 |
|
|
while(pp<pts+3*unit_voro.p) { |
| 68 |
|
|
q=*(pp++);y=*(pp++);z=*(pp++);q=sqrt(q*q+y*y+z*z); |
| 69 |
|
|
if(y+q>max_uv_y) max_uv_y=y+q; |
| 70 |
|
|
if(z+q>max_uv_z) max_uv_z=z+q; |
| 71 |
|
|
} |
| 72 |
|
|
max_uv_z*=0.5; |
| 73 |
|
|
max_uv_y*=0.5; |
| 74 |
|
|
return; |
| 75 |
|
|
} |
| 76 |
|
|
l++; |
| 77 |
|
|
} |
| 78 |
|
|
|
| 79 |
|
|
// If the routine makes it here, then the unit cell still hasn't been |
| 80 |
|
|
// completely bounded by the plane cuts. Give the memory error code, |
| 81 |
|
|
// because this is mainly a case of hitting a safe limit, than any |
| 82 |
|
|
// inherent problem. |
| 83 |
|
|
voro_fatal_error("Periodic cell computation failed",VOROPP_MEMORY_ERROR); |
| 84 |
|
|
} |
| 85 |
|
|
|
| 86 |
|
|
/** Applies a pair of opposing plane cuts from a periodic image point |
| 87 |
|
|
* to the unit Voronoi cell. |
| 88 |
|
|
* \param[in] (i,j,k) the index of the periodic image to consider. */ |
| 89 |
|
|
inline void unitcell::unit_voro_apply(int i,int j,int k) { |
| 90 |
|
|
double x=i*bx+j*bxy+k*bxz,y=j*by+k*byz,z=k*bz; |
| 91 |
|
|
unit_voro.plane(x,y,z); |
| 92 |
|
|
unit_voro.plane(-x,-y,-z); |
| 93 |
|
|
} |
| 94 |
|
|
|
| 95 |
|
|
/** Calculates whether the unit Voronoi cell intersects a given periodic image |
| 96 |
|
|
* of the domain. |
| 97 |
|
|
* \param[in] (dx,dy,dz) the displacement of the periodic image. |
| 98 |
|
|
* \param[out] vol the proportion of the unit cell volume within this image, |
| 99 |
|
|
* only computed in the case that the two intersect. |
| 100 |
|
|
* \return True if they intersect, false otherwise. */ |
| 101 |
|
|
bool unitcell::intersects_image(double dx,double dy,double dz,double &vol) { |
| 102 |
|
|
const double bxinv=1/bx,byinv=1/by,bzinv=1/bz,ivol=bxinv*byinv*bzinv; |
| 103 |
|
|
voronoicell c; |
| 104 |
|
|
c=unit_voro; |
| 105 |
|
|
dx*=2;dy*=2;dz*=2; |
| 106 |
|
|
if(!c.plane(0,0,bzinv,dz+1)) return false; |
| 107 |
|
|
if(!c.plane(0,0,-bzinv,-dz+1)) return false; |
| 108 |
|
|
if(!c.plane(0,byinv,-byz*byinv*bzinv,dy+1)) return false; |
| 109 |
|
|
if(!c.plane(0,-byinv,byz*byinv*bzinv,-dy+1)) return false; |
| 110 |
|
|
if(!c.plane(bxinv,-bxy*bxinv*byinv,(bxy*byz-by*bxz)*ivol,dx+1)) return false; |
| 111 |
|
|
if(!c.plane(-bxinv,bxy*bxinv*byinv,(-bxy*byz+by*bxz)*ivol,-dx+1)) return false; |
| 112 |
|
|
vol=c.volume()*ivol; |
| 113 |
|
|
return true; |
| 114 |
|
|
} |
| 115 |
|
|
|
| 116 |
|
|
/** Computes a list of periodic domain images that intersect the unit Voronoi cell. |
| 117 |
|
|
* \param[out] vi a vector containing triplets (i,j,k) corresponding to domain |
| 118 |
|
|
* images that intersect the unit Voronoi cell, when it is |
| 119 |
|
|
* centered in the middle of the primary domain. |
| 120 |
|
|
* \param[out] vd a vector containing the fraction of the Voronoi cell volume |
| 121 |
|
|
* within each corresponding image listed in vi. */ |
| 122 |
|
|
void unitcell::images(std::vector<int> &vi,std::vector<double> &vd) { |
| 123 |
|
|
const int ms2=max_unit_voro_shells*2+1,mss=ms2*ms2*ms2; |
| 124 |
|
|
bool *a=new bool[mss],*ac=a+max_unit_voro_shells*(1+ms2*(1+ms2)),*ap=a; |
| 125 |
|
|
int i,j,k; |
| 126 |
|
|
double vol; |
| 127 |
|
|
|
| 128 |
|
|
// Initialize mask |
| 129 |
|
|
while(ap<ac) *(ap++)=true; |
| 130 |
|
|
*(ap++)=false; |
| 131 |
|
|
while(ap<a+mss) *(ap++)=true; |
| 132 |
|
|
|
| 133 |
|
|
// Set up the queue and add (0,0,0) image to it |
| 134 |
|
|
std::queue<int> q; |
| 135 |
|
|
q.push(0);q.push(0);q.push(0); |
| 136 |
|
|
|
| 137 |
|
|
while(!q.empty()) { |
| 138 |
|
|
|
| 139 |
|
|
// Read the next entry on the queue |
| 140 |
|
|
i=q.front();q.pop(); |
| 141 |
|
|
j=q.front();q.pop(); |
| 142 |
|
|
k=q.front();q.pop(); |
| 143 |
|
|
|
| 144 |
|
|
// Check intersection of this image |
| 145 |
|
|
if(intersects_image(i,j,k,vol)) { |
| 146 |
|
|
|
| 147 |
|
|
// Add this entry to the output vectors |
| 148 |
|
|
vi.push_back(i); |
| 149 |
|
|
vi.push_back(j); |
| 150 |
|
|
vi.push_back(k); |
| 151 |
|
|
vd.push_back(vol); |
| 152 |
|
|
|
| 153 |
|
|
// Add neighbors to the queue if they have not been |
| 154 |
|
|
// tested |
| 155 |
|
|
ap=ac+i+ms2*(j+ms2*k); |
| 156 |
|
|
if(k>-max_unit_voro_shells&&*(ap-ms2*ms2)) {q.push(i);q.push(j);q.push(k-1);*(ap-ms2*ms2)=false;} |
| 157 |
|
|
if(j>-max_unit_voro_shells&&*(ap-ms2)) {q.push(i);q.push(j-1);q.push(k);*(ap-ms2)=false;} |
| 158 |
|
|
if(i>-max_unit_voro_shells&&*(ap-1)) {q.push(i-1);q.push(j);q.push(k);*(ap-1)=false;} |
| 159 |
|
|
if(i<max_unit_voro_shells&&*(ap+1)) {q.push(i+1);q.push(j);q.push(k);*(ap+1)=false;} |
| 160 |
|
|
if(j<max_unit_voro_shells&&*(ap+ms2)) {q.push(i);q.push(j+1);q.push(k);*(ap+ms2)=false;} |
| 161 |
|
|
if(k<max_unit_voro_shells&&*(ap+ms2*ms2)) {q.push(i);q.push(j);q.push(k+1);*(ap+ms2*ms2)=false;} |
| 162 |
|
|
} |
| 163 |
|
|
} |
| 164 |
|
|
|
| 165 |
|
|
// Remove mask memory |
| 166 |
|
|
delete [] a; |
| 167 |
|
|
} |
| 168 |
|
|
|
| 169 |
|
|
/** Tests to see if a shell of periodic images could possibly cut the periodic |
| 170 |
|
|
* unit cell. |
| 171 |
|
|
* \param[in] l the index of the shell to consider. |
| 172 |
|
|
* \return True if a point in the shell cuts the cell, false otherwise. */ |
| 173 |
|
|
bool unitcell::unit_voro_intersect(int l) { |
| 174 |
|
|
int i,j; |
| 175 |
|
|
if(unit_voro_test(l,0,0)) return true; |
| 176 |
|
|
for(i=1;i<l;i++) { |
| 177 |
|
|
if(unit_voro_test(l,i,0)) return true; |
| 178 |
|
|
if(unit_voro_test(-l,i,0)) return true; |
| 179 |
|
|
} |
| 180 |
|
|
for(i=-l;i<=l;i++) if(unit_voro_test(i,l,0)) return true; |
| 181 |
|
|
for(i=1;i<l;i++) for(j=-l+1;j<=l;j++) { |
| 182 |
|
|
if(unit_voro_test(l,j,i)) return true; |
| 183 |
|
|
if(unit_voro_test(-j,l,i)) return true; |
| 184 |
|
|
if(unit_voro_test(-l,-j,i)) return true; |
| 185 |
|
|
if(unit_voro_test(j,-l,i)) return true; |
| 186 |
|
|
} |
| 187 |
|
|
for(i=-l;i<=l;i++) for(j=-l;j<=l;j++) if(unit_voro_test(i,j,l)) return true; |
| 188 |
|
|
return false; |
| 189 |
|
|
} |
| 190 |
|
|
|
| 191 |
|
|
/** Tests to see if a plane cut from a particular periodic image will cut th |
| 192 |
|
|
* unit Voronoi cell. |
| 193 |
|
|
* \param[in] (i,j,k) the index of the periodic image to consider. |
| 194 |
|
|
* \return True if the image cuts the cell, false otherwise. */ |
| 195 |
|
|
inline bool unitcell::unit_voro_test(int i,int j,int k) { |
| 196 |
|
|
double x=i*bx+j*bxy+k*bxz,y=j*by+k*byz,z=k*bz; |
| 197 |
|
|
double rsq=x*x+y*y+z*z; |
| 198 |
|
|
return unit_voro.plane_intersects(x,y,z,rsq); |
| 199 |
|
|
} |
| 200 |
|
|
|
| 201 |
|
|
/** Draws the periodic domain in gnuplot format. |
| 202 |
|
|
* \param[in] fp the file handle to write to. */ |
| 203 |
|
|
void unitcell::draw_domain_gnuplot(FILE *fp) { |
| 204 |
|
|
fprintf(fp,"0 0 0\n%g 0 0\n%g %g 0\n%g %g 0\n",bx,bx+bxy,by,bxy,by); |
| 205 |
|
|
fprintf(fp,"%g %g %g\n%g %g %g\n%g %g %g\n%g %g %g\n",bxy+bxz,by+byz,bz,bx+bxy+bxz,by+byz,bz,bx+bxz,byz,bz,bxz,byz,bz); |
| 206 |
|
|
fprintf(fp,"0 0 0\n%g %g 0\n\n%g %g %g\n%g %g %g\n\n",bxy,by,bxz,byz,bz,bxy+bxz,by+byz,bz); |
| 207 |
|
|
fprintf(fp,"%g 0 0\n%g %g %g\n\n%g %g 0\n%g %g %g\n\n",bx,bx+bxz,byz,bz,bx+bxy,by,bx+bxy+bxz,by+byz,bz); |
| 208 |
|
|
} |
| 209 |
|
|
|
| 210 |
|
|
/** Draws the periodic domain in POV-Ray format. |
| 211 |
|
|
* \param[in] fp the file handle to write to. */ |
| 212 |
|
|
void unitcell::draw_domain_pov(FILE *fp) { |
| 213 |
|
|
fprintf(fp,"cylinder{0,0,0>,<%g,0,0>,rr}\n" |
| 214 |
|
|
"cylinder{<%g,%g,0>,<%g,%g,0>,rr}\n",bx,bxy,by,bx+bxy,by); |
| 215 |
|
|
fprintf(fp,"cylinder{<%g,%g,%g>,<%g,%g,%g>,rr}\n" |
| 216 |
|
|
"cylinder{<%g,%g,%g>,<%g,%g,%g>,rr}\n",bxz,byz,bz,bx+bxz,byz,bz,bxy+bxz,by+byz,bz,bx+bxy+bxz,by+byz,bz); |
| 217 |
|
|
fprintf(fp,"cylinder{<0,0,0>,<%g,%g,0>,rr}\n" |
| 218 |
|
|
"cylinder{<%g,0,0>,<%g,%g,0>,rr}\n",bxy,by,bx,bx+bxy,by); |
| 219 |
|
|
fprintf(fp,"cylinder{<%g,%g,%g>,<%g,%g,%g>,rr}\n" |
| 220 |
|
|
"cylinder{<%g,%g,%g>,<%g,%g,%g>,rr}\n",bxz,byz,bz,bxy+bxz,by+byz,bz,bx+bxz,byz,bz,bx+bxy+bxz,by+byz,bz); |
| 221 |
|
|
fprintf(fp,"cylinder{<0,0,0>,<%g,%g,%g>,rr}\n" |
| 222 |
|
|
"cylinder{<%g,0,0>,<%g,%g,%g>,rr}\n",bxz,byz,bz,bx,bx+bxz,byz,bz); |
| 223 |
|
|
fprintf(fp,"cylinder{<%g,%g,0>,<%g,%g,%g>,rr}\n" |
| 224 |
|
|
"cylinder{<%g,%g,0>,<%g,%g,%g>,rr}\n",bxy,by,bxy+bxz,by+byz,bz,bx+bxy,by,bx+bxy+bxz,by+byz,bz); |
| 225 |
|
|
fprintf(fp,"sphere{<0,0,0>,rr}\nsphere{<%g,0,0>,rr}\n" |
| 226 |
|
|
"sphere{<%g,%g,0>,rr}\nsphere{<%g,%g,0>,rr}\n",bx,bxy,by,bx+bxy,by); |
| 227 |
|
|
fprintf(fp,"sphere{<%g,%g,%g>,rr}\nsphere{<%g,%g,%g>,rr}\n" |
| 228 |
|
|
"sphere{<%g,%g,%g>,rr}\nsphere{<%g,%g,%g>,rr}\n",bxz,byz,bz,bx+bxz,byz,bz,bxy+bxz,by+byz,bz,bx+bxy+bxz,by+byz,bz); |
| 229 |
|
|
} |
| 230 |
|
|
|
| 231 |
|
|
} |
