-
Notifications
You must be signed in to change notification settings - Fork 10
/
cylinder_inv.cc
77 lines (65 loc) · 2.44 KB
/
cylinder_inv.cc
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
// Cylindrical wall example code
//
// Author : Chris H. Rycroft (LBL / UC Berkeley)
// Email : chr@alum.mit.edu
// Date : August 30th 2011
#include "voro++.hh"
using namespace voro;
// Set up constants for the container geometry
const double x_min=-6,x_max=6;
const double y_min=-6,y_max=6;
const double z_min=-6,z_max=6;
// Set the computational grid size
const int n_x=6,n_y=6,n_z=6;
struct wall_cylinder_inv : public wall {
public:
wall_cylinder_inv(double xc_,double yc_,double zc_,double xa_,double ya_,double za_,double rc_,int w_id_=-99)
: w_id(w_id_), xc(xc_), yc(yc_), zc(zc_), xa(xa_), ya(ya_), za(za_),
asi(1/(xa_*xa_+ya_*ya_+za_*za_)), rc(rc_) {}
bool point_inside(double x,double y,double z) {
double xd=x-xc,yd=y-yc,zd=z-zc;
double pa=(xd*xa+yd*ya+zd*za)*asi;
xd-=xa*pa;yd-=ya*pa;zd-=za*pa;
return xd*xd+yd*yd+zd*zd>rc*rc;
}
template<class v_cell>
bool cut_cell_base(v_cell &c,double x,double y,double z) {
double xd=x-xc,yd=y-yc,zd=z-zc;
double pa=(xd*xa+yd*ya+zd*za)*asi;
xd-=xa*pa;yd-=ya*pa;zd-=za*pa;
pa=xd*xd+yd*yd+zd*zd;
if(pa>1e-5) {
pa=2*(sqrt(pa)*rc-pa);
return c.nplane(-xd,-yd,-zd,-pa,w_id);
}
return true;
}
bool cut_cell(voronoicell &c,double x,double y,double z) {return cut_cell_base(c,x,y,z);}
bool cut_cell(voronoicell_neighbor &c,double x,double y,double z) {return cut_cell_base(c,x,y,z);}
private:
const int w_id;
const double xc,yc,zc,xa,ya,za,asi,rc;
};
int main() {
int i;double x,y,z;
// Create a container with the geometry given above, and make it
// non-periodic in each of the three coordinates. Allocate space for
// eight particles within each computational block.
container con(x_min,x_max,y_min,y_max,z_min,z_max,n_x,n_y,n_z,
false,false,false,8);
// Add a cylindrical wall to the container
wall_cylinder_inv cyl(0,0,0,0,1,0,4.2);
con.add_wall(cyl);
// Place particles in a regular grid within the frustum, for points
// which are within the wall boundaries
for(z=-5.5;z<6;z+=1) for(y=-5.5;y<6;y+=1) for(x=-5.5;x<6;x+=1)
if (con.point_inside(x,y,z)) {con.put(i,x,y,z);i++;}
// Output the particle positions in POV-Ray format
con.draw_particles_pov("cylinder_inv_p.pov");
// Output the Voronoi cells in POV-Ray format
con.draw_cells_pov("cylinder_inv_v.pov");
// Output the particle positions in gnuplot format
con.draw_particles("cylinder_inv.par");
// Output the Voronoi cells in gnuplot format
con.draw_cells_gnuplot("cylinder_inv.gnu");
}