voroxx 0.1.0

Voro++, a 3D cell-based Voronoi library
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
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236

// Voro++, a 3D cell-based Voronoi library
//
// Author   : Chris H. Rycroft (LBL / UC Berkeley)
// Email    : chr@alum.mit.edu
// Date     : August 30th 2011

/** \file pre_container.cc
 * \brief Function implementations for the pre_container and related classes.
 */

#include <cmath>

#include "config.hh"
#include "pre_container.hh"

namespace voro {

/** The class constructor sets up the geometry of container, initializing the
 * minimum and maximum coordinates in each direction. It allocates an initial
 * chunk into which to store particle information.
 * \param[in] (ax_,bx_) the minimum and maximum x coordinates.
 * \param[in] (ay_,by_) the minimum and maximum y coordinates.
 * \param[in] (az_,bz_) the minimum and maximum z coordinates.
 * \param[in] (xperiodic_,yperiodic_,zperiodic_ ) flags setting whether the
 *                                                container is periodic in each
 *                                                coordinate direction.
 * \param[in] ps_ the number of floating point entries to store for each
 *                particle. */
pre_container_base::pre_container_base(double ax_,double bx_,double ay_,double by_,double az_,double bz_,
	bool xperiodic_,bool yperiodic_,bool zperiodic_,int ps_) :
	ax(ax_), bx(bx_), ay(ay_), by(by_), az(az_), bz(bz_),
	xperiodic(xperiodic_), yperiodic(yperiodic_), zperiodic(zperiodic_), ps(ps_),
	index_sz(init_chunk_size), pre_id(new int*[index_sz]), end_id(pre_id),
	pre_p(new double*[index_sz]), end_p(pre_p) {
		ch_id=*end_id=new int[pre_container_chunk_size];
		l_id=end_id+index_sz;e_id=ch_id+pre_container_chunk_size;
		ch_p=*end_p=new double[ps*pre_container_chunk_size];
}

/** The destructor frees the dynamically allocated memory. */
pre_container_base::~pre_container_base() {
	delete [] *end_p;
	delete [] *end_id;
	while (end_id!=pre_id) {
		end_p--;
		delete [] *end_p;
		end_id--;
		delete [] *end_id;
	}
	delete [] pre_p;
	delete [] pre_id;
}

/** Makes a guess at the optimal grid of blocks to use, computing in
 * a way that
 * \param[out] (nx,ny,nz) the number of blocks to use. */
void pre_container_base::guess_optimal(int &nx,int &ny,int &nz) {
	double dx=bx-ax,dy=by-ay,dz=bz-az;
	double ilscale=pow(total_particles()/(optimal_particles*dx*dy*dz),1/3.0);
	nx=int(dx*ilscale+1);
	ny=int(dy*ilscale+1);
	nz=int(dz*ilscale+1);
}

/** Stores a particle ID and position, allocating a new memory chunk if
 * necessary. For coordinate directions in which the container is not periodic,
 * the routine checks to make sure that the particle is within the container
 * bounds. If the particle is out of bounds, it is not stored.
 * \param[in] n the numerical ID of the inserted particle.
 * \param[in] (x,y,z) the position vector of the inserted particle. */
void pre_container::put(int n,double x,double y,double z) {
	if((xperiodic||(x>=ax&&x<=bx))&&(yperiodic||(y>=ay&&y<=by))&&(zperiodic||(z>=az&&z<=bz))) {
		if(ch_id==e_id) new_chunk();
		*(ch_id++)=n;
		*(ch_p++)=x;*(ch_p++)=y;*(ch_p++)=z;
	}
#if VOROPP_REPORT_OUT_OF_BOUNDS ==1
	else fprintf(stderr,"Out of bounds: (x,y,z)=(%g,%g,%g)\n",x,y,z);
#endif
}

/** Stores a particle ID and position, allocating a new memory chunk if necessary.
 * \param[in] n the numerical ID of the inserted particle.
 * \param[in] (x,y,z) the position vector of the inserted particle.
 * \param[in] r the radius of the particle. */
void pre_container_poly::put(int n,double x,double y,double z,double r) {
	if((xperiodic||(x>=ax&&x<=bx))&&(yperiodic||(y>=ay&&y<=by))&&(zperiodic||(z>=az&&z<=bz))) {
		if(ch_id==e_id) new_chunk();
		*(ch_id++)=n;
		*(ch_p++)=x;*(ch_p++)=y;*(ch_p++)=z;*(ch_p++)=r;
	}
#if VOROPP_REPORT_OUT_OF_BOUNDS ==1
	else fprintf(stderr,"Out of bounds: (x,y,z)=(%g,%g,%g)\n",x,y,z);
#endif
}

/** Transfers the particles stored within the class to a container class.
 * \param[in] con the container class to transfer to. */
void pre_container::setup(container &con) {
	int **c_id=pre_id,*idp,*ide,n;
	double **c_p=pre_p,*pp,x,y,z;
	while(c_id<end_id) {
		idp=*(c_id++);ide=idp+pre_container_chunk_size;
		pp=*(c_p++);
		while(idp<ide) {
			n=*(idp++);x=*(pp++);y=*(pp++);z=*(pp++);
			con.put(n,x,y,z);
		}
	}
	idp=*c_id;
	pp=*c_p;
	while(idp<ch_id) {
		n=*(idp++);x=*(pp++);y=*(pp++);z=*(pp++);
		con.put(n,x,y,z);
	}
}

/** Transfers the particles stored within the class to a container_poly class.
 * \param[in] con the container_poly class to transfer to. */
void pre_container_poly::setup(container_poly &con) {
	int **c_id=pre_id,*idp,*ide,n;
	double **c_p=pre_p,*pp,x,y,z,r;
	while(c_id<end_id) {
		idp=*(c_id++);ide=idp+pre_container_chunk_size;
		pp=*(c_p++);
		while(idp<ide) {
			n=*(idp++);x=*(pp++);y=*(pp++);z=*(pp++);r=*(pp++);
			con.put(n,x,y,z,r);
		}
	}
	idp=*c_id;
	pp=*c_p;
	while(idp<ch_id) {
		n=*(idp++);x=*(pp++);y=*(pp++);z=*(pp++);r=*(pp++);
		con.put(n,x,y,z,r);
	}
}

/** Transfers the particles stored within the class to a container class, also
 * recording the order in which particles were stored.
 * \param[in] vo the ordering class to use.
 * \param[in] con the container class to transfer to. */
void pre_container::setup(particle_order &vo,container &con) {
	int **c_id=pre_id,*idp,*ide,n;
	double **c_p=pre_p,*pp,x,y,z;
	while(c_id<end_id) {
		idp=*(c_id++);ide=idp+pre_container_chunk_size;
		pp=*(c_p++);
		while(idp<ide) {
			n=*(idp++);x=*(pp++);y=*(pp++);z=*(pp++);
			con.put(vo,n,x,y,z);
		}
	}
	idp=*c_id;
	pp=*c_p;
	while(idp<ch_id) {
		n=*(idp++);x=*(pp++);y=*(pp++);z=*(pp++);
		con.put(vo,n,x,y,z);
	}
}

/** Transfers the particles stored within the class to a container_poly class,
 * also recording the order in which particles were stored.
 * \param[in] vo the ordering class to use.
 * \param[in] con the container_poly class to transfer to. */
void pre_container_poly::setup(particle_order &vo,container_poly &con) {
	int **c_id=pre_id,*idp,*ide,n;
	double **c_p=pre_p,*pp,x,y,z,r;
	while(c_id<end_id) {
		idp=*(c_id++);ide=idp+pre_container_chunk_size;
		pp=*(c_p++);
		while(idp<ide) {
			n=*(idp++);x=*(pp++);y=*(pp++);z=*(pp++);r=*(pp++);
			con.put(vo,n,x,y,z,r);
		}
	}
	idp=*c_id;
	pp=*c_p;
	while(idp<ch_id) {
		n=*(idp++);x=*(pp++);y=*(pp++);z=*(pp++);r=*(pp++);
		con.put(vo,n,x,y,z,r);
	}
}

/** Import a list of particles from an open file stream into the container.
 * Entries of four numbers (Particle ID, x position, y position, z position)
 * are searched for. If the file cannot be successfully read, then the routine
 * causes a fatal error.
 * \param[in] fp the file handle to read from. */
void pre_container::import(FILE *fp) {
	int i,j;
	double x,y,z;
	while((j=fscanf(fp,"%d %lg %lg %lg",&i,&x,&y,&z))==4) put(i,x,y,z);
	if(j!=EOF) voro_fatal_error("File import error",VOROPP_FILE_ERROR);
}

/** Import a list of particles from an open file stream, also storing the order
 * of that the particles are read. Entries of four numbers (Particle ID, x
 * position, y position, z position) are searched for. If the file cannot be
 * successfully read, then the routine causes a fatal error.
 * \param[in] fp the file handle to read from. */
void pre_container_poly::import(FILE *fp) {
	int i,j;
	double x,y,z,r;
	while((j=fscanf(fp,"%d %lg %lg %lg %lg",&i,&x,&y,&z,&r))==5) put(i,x,y,z,r);
	if(j!=EOF) voro_fatal_error("File import error",VOROPP_FILE_ERROR);
}

/** Allocates a new chunk of memory for storing particles. */
void pre_container_base::new_chunk() {
	end_id++;end_p++;
	if(end_id==l_id) extend_chunk_index();
	ch_id=*end_id=new int[pre_container_chunk_size];
	e_id=ch_id+pre_container_chunk_size;
	ch_p=*end_p=new double[ps*pre_container_chunk_size];
}

/** Extends the index of chunks. */
void pre_container_base::extend_chunk_index() {
	index_sz<<=1;
	if(index_sz>max_chunk_size)
		voro_fatal_error("Absolute memory limit on chunk index reached",VOROPP_MEMORY_ERROR);
#if VOROPP_VERBOSE >=2
	fprintf(stderr,"Pre-container chunk index scaled up to %d\n",index_sz);
#endif
	int **n_id=new int*[index_sz],**p_id=n_id,**c_id=pre_id;
	double **n_p=new double*[index_sz],**p_p=n_p,**c_p=pre_p;
	while(c_id<end_id) {
		*(p_id++)=*(c_id++);
		*(p_p++)=*(c_p++);
	}
	delete [] pre_id;pre_id=n_id;end_id=p_id;l_id=pre_id+index_sz;
	delete [] pre_p;pre_p=n_p;end_p=p_p;
}

}