| 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 cell.hh |
| 8 |
|
|
* \brief Header file for the voronoicell and related classes. */ |
| 9 |
|
|
|
| 10 |
|
|
#ifndef VOROPP_CELL_HH |
| 11 |
|
|
#define VOROPP_CELL_HH |
| 12 |
|
|
|
| 13 |
|
|
#include <vector> |
| 14 |
|
|
|
| 15 |
|
|
#include "config.hh" |
| 16 |
|
|
#include "common.hh" |
| 17 |
|
|
|
| 18 |
|
|
namespace voro { |
| 19 |
|
|
|
| 20 |
|
|
/** \brief A class representing a single Voronoi cell. |
| 21 |
|
|
* |
| 22 |
|
|
* This class represents a single Voronoi cell, as a collection of vertices |
| 23 |
|
|
* that are connected by edges. The class contains routines for initializing |
| 24 |
|
|
* the Voronoi cell to be simple shapes such as a box, tetrahedron, or octahedron. |
| 25 |
|
|
* It the contains routines for recomputing the cell based on cutting it |
| 26 |
|
|
* by a plane, which forms the key routine for the Voronoi cell computation. |
| 27 |
|
|
* It contains numerous routine for computing statistics about the Voronoi cell, |
| 28 |
|
|
* and it can output the cell in several formats. |
| 29 |
|
|
* |
| 30 |
|
|
* This class is not intended for direct use, but forms the base of the |
| 31 |
|
|
* voronoicell and voronoicell_neighbor classes, which extend it based on |
| 32 |
|
|
* whether neighboring particle ID information needs to be tracked. */ |
| 33 |
|
|
class voronoicell_base { |
| 34 |
|
|
public: |
| 35 |
|
|
/** This holds the current size of the arrays ed and nu, which |
| 36 |
|
|
* hold the vertex information. If more vertices are created |
| 37 |
|
|
* than can fit in this array, then it is dynamically extended |
| 38 |
|
|
* using the add_memory_vertices routine. */ |
| 39 |
|
|
int current_vertices; |
| 40 |
|
|
/** This holds the current maximum allowed order of a vertex, |
| 41 |
|
|
* which sets the size of the mem, mep, and mec arrays. If a |
| 42 |
|
|
* vertex is created with more vertices than this, the arrays |
| 43 |
|
|
* are dynamically extended using the add_memory_vorder routine. |
| 44 |
|
|
*/ |
| 45 |
|
|
int current_vertex_order; |
| 46 |
|
|
/** This sets the size of the main delete stack. */ |
| 47 |
|
|
int current_delete_size; |
| 48 |
|
|
/** This sets the size of the auxiliary delete stack. */ |
| 49 |
|
|
int current_delete2_size; |
| 50 |
|
|
/** This sets the total number of vertices in the current cell. |
| 51 |
|
|
*/ |
| 52 |
|
|
int p; |
| 53 |
|
|
/** This is the index of particular point in the cell, which is |
| 54 |
|
|
* used to start the tracing routines for plane intersection |
| 55 |
|
|
* and cutting. These routines will work starting from any |
| 56 |
|
|
* point, but it's often most efficient to start from the last |
| 57 |
|
|
* point considered, since in many cases, the cell construction |
| 58 |
|
|
* algorithm may consider many planes with similar vectors |
| 59 |
|
|
* concurrently. */ |
| 60 |
|
|
int up; |
| 61 |
|
|
/** This is a two dimensional array that holds information |
| 62 |
|
|
* about the edge connections of the vertices that make up the |
| 63 |
|
|
* cell. The two dimensional array is not allocated in the |
| 64 |
|
|
* usual method. To account for the fact the different vertices |
| 65 |
|
|
* have different orders, and thus require different amounts of |
| 66 |
|
|
* storage, the elements of ed[i] point to one-dimensional |
| 67 |
|
|
* arrays in the mep[] array of different sizes. |
| 68 |
|
|
* |
| 69 |
|
|
* More specifically, if vertex i has order m, then ed[i] |
| 70 |
|
|
* points to a one-dimensional array in mep[m] that has 2*m+1 |
| 71 |
|
|
* entries. The first m elements hold the neighboring edges, so |
| 72 |
|
|
* that the jth edge of vertex i is held in ed[i][j]. The next |
| 73 |
|
|
* m elements hold a table of relations which is redundant but |
| 74 |
|
|
* helps speed up the computation. It satisfies the relation |
| 75 |
|
|
* ed[ed[i][j]][ed[i][m+j]]=i. The final entry holds a back |
| 76 |
|
|
* pointer, so that ed[i+2*m]=i. The back pointers are used |
| 77 |
|
|
* when rearranging the memory. */ |
| 78 |
|
|
int **ed; |
| 79 |
|
|
/** This array holds the order of the vertices in the Voronoi |
| 80 |
|
|
* cell. This array is dynamically allocated, with its current |
| 81 |
|
|
* size held by current_vertices. */ |
| 82 |
|
|
int *nu; |
| 83 |
|
|
/** This in an array with size 3*current_vertices for holding |
| 84 |
|
|
* the positions of the vertices. */ |
| 85 |
|
|
double *pts; |
| 86 |
|
|
voronoicell_base(); |
| 87 |
|
|
~voronoicell_base(); |
| 88 |
|
|
void init_base(double xmin,double xmax,double ymin,double ymax,double zmin,double zmax); |
| 89 |
|
|
void init_octahedron_base(double l); |
| 90 |
|
|
void init_tetrahedron_base(double x0,double y0,double z0,double x1,double y1,double z1,double x2,double y2,double z2,double x3,double y3,double z3); |
| 91 |
|
|
void translate(double x,double y,double z); |
| 92 |
|
|
void draw_pov(double x,double y,double z,FILE *fp=stdout); |
| 93 |
|
|
/** Outputs the cell in POV-Ray format, using cylinders for edges |
| 94 |
|
|
* and spheres for vertices, to a given file. |
| 95 |
|
|
* \param[in] (x,y,z) a displacement to add to the cell's |
| 96 |
|
|
* position. |
| 97 |
|
|
* \param[in] filename the name of the file to write to. */ |
| 98 |
|
|
inline void draw_pov(double x,double y,double z,const char *filename) { |
| 99 |
|
|
FILE *fp=safe_fopen(filename,"w"); |
| 100 |
|
|
draw_pov(x,y,z,fp); |
| 101 |
|
|
fclose(fp); |
| 102 |
|
|
}; |
| 103 |
|
|
void draw_pov_mesh(double x,double y,double z,FILE *fp=stdout); |
| 104 |
|
|
/** Outputs the cell in POV-Ray format as a mesh2 object to a |
| 105 |
|
|
* given file. |
| 106 |
|
|
* \param[in] (x,y,z) a displacement to add to the cell's |
| 107 |
|
|
* position. |
| 108 |
|
|
* \param[in] filename the name of the file to write to. */ |
| 109 |
|
|
inline void draw_pov_mesh(double x,double y,double z,const char *filename) { |
| 110 |
|
|
FILE *fp=safe_fopen(filename,"w"); |
| 111 |
|
|
draw_pov_mesh(x,y,z,fp); |
| 112 |
|
|
fclose(fp); |
| 113 |
|
|
} |
| 114 |
|
|
void draw_gnuplot(double x,double y,double z,FILE *fp=stdout); |
| 115 |
|
|
/** Outputs the cell in Gnuplot format a given file. |
| 116 |
|
|
* \param[in] (x,y,z) a displacement to add to the cell's |
| 117 |
|
|
* position. |
| 118 |
|
|
* \param[in] filename the name of the file to write to. */ |
| 119 |
|
|
inline void draw_gnuplot(double x,double y,double z,const char *filename) { |
| 120 |
|
|
FILE *fp=safe_fopen(filename,"w"); |
| 121 |
|
|
draw_gnuplot(x,y,z,fp); |
| 122 |
|
|
fclose(fp); |
| 123 |
|
|
} |
| 124 |
|
|
double volume(); |
| 125 |
|
|
double max_radius_squared(); |
| 126 |
|
|
double total_edge_distance(); |
| 127 |
|
|
double surface_area(); |
| 128 |
|
|
void centroid(double &cx,double &cy,double &cz); |
| 129 |
|
|
int number_of_faces(); |
| 130 |
|
|
int number_of_edges(); |
| 131 |
|
|
void vertex_orders(std::vector<int> &v); |
| 132 |
|
|
void output_vertex_orders(FILE *fp=stdout); |
| 133 |
|
|
void vertices(std::vector<double> &v); |
| 134 |
|
|
void output_vertices(FILE *fp=stdout); |
| 135 |
|
|
void vertices(double x,double y,double z,std::vector<double> &v); |
| 136 |
|
|
void output_vertices(double x,double y,double z,FILE *fp=stdout); |
| 137 |
|
|
void face_areas(std::vector<double> &v); |
| 138 |
|
|
/** Outputs the areas of the faces. |
| 139 |
|
|
* \param[in] fp the file handle to write to. */ |
| 140 |
|
|
inline void output_face_areas(FILE *fp=stdout) { |
| 141 |
|
|
std::vector<double> v;face_areas(v); |
| 142 |
|
|
voro_print_vector(v,fp); |
| 143 |
|
|
} |
| 144 |
|
|
void face_orders(std::vector<int> &v); |
| 145 |
|
|
/** Outputs a list of the number of sides of each face. |
| 146 |
|
|
* \param[in] fp the file handle to write to. */ |
| 147 |
|
|
inline void output_face_orders(FILE *fp=stdout) { |
| 148 |
|
|
std::vector<int> v;face_orders(v); |
| 149 |
|
|
voro_print_vector(v,fp); |
| 150 |
|
|
} |
| 151 |
|
|
void face_freq_table(std::vector<int> &v); |
| 152 |
|
|
/** Outputs a */ |
| 153 |
|
|
inline void output_face_freq_table(FILE *fp=stdout) { |
| 154 |
|
|
std::vector<int> v;face_freq_table(v); |
| 155 |
|
|
voro_print_vector(v,fp); |
| 156 |
|
|
} |
| 157 |
|
|
void face_vertices(std::vector<int> &v); |
| 158 |
|
|
/** Outputs the */ |
| 159 |
|
|
inline void output_face_vertices(FILE *fp=stdout) { |
| 160 |
|
|
std::vector<int> v;face_vertices(v); |
| 161 |
|
|
voro_print_face_vertices(v,fp); |
| 162 |
|
|
} |
| 163 |
|
|
void face_perimeters(std::vector<double> &v); |
| 164 |
|
|
/** Outputs a list of the perimeters of each face. |
| 165 |
|
|
* \param[in] fp the file handle to write to. */ |
| 166 |
|
|
inline void output_face_perimeters(FILE *fp=stdout) { |
| 167 |
|
|
std::vector<double> v;face_perimeters(v); |
| 168 |
|
|
voro_print_vector(v,fp); |
| 169 |
|
|
} |
| 170 |
|
|
void normals(std::vector<double> &v); |
| 171 |
|
|
/** Outputs a list of the perimeters of each face. |
| 172 |
|
|
* \param[in] fp the file handle to write to. */ |
| 173 |
|
|
inline void output_normals(FILE *fp=stdout) { |
| 174 |
|
|
std::vector<double> v;normals(v); |
| 175 |
|
|
voro_print_positions(v,fp); |
| 176 |
|
|
} |
| 177 |
|
|
/** Outputs a custom string of information about the Voronoi |
| 178 |
|
|
* cell to a file. It assumes the cell is at (0,0,0) and has a |
| 179 |
|
|
* the default_radius associated with it. |
| 180 |
|
|
* \param[in] format the custom format string to use. |
| 181 |
|
|
* \param[in] fp the file handle to write to. */ |
| 182 |
|
|
inline void output_custom(const char *format,FILE *fp=stdout) {output_custom(format,0,0,0,0,default_radius,fp);} |
| 183 |
|
|
void output_custom(const char *format,int i,double x,double y,double z,double r,FILE *fp=stdout); |
| 184 |
|
|
template<class vc_class> |
| 185 |
|
|
bool nplane(vc_class &vc,double x,double y,double z,double rsq,int p_id); |
| 186 |
|
|
bool plane_intersects(double x,double y,double z,double rsq); |
| 187 |
|
|
bool plane_intersects_guess(double x,double y,double z,double rsq); |
| 188 |
|
|
void construct_relations(); |
| 189 |
|
|
void check_relations(); |
| 190 |
|
|
void check_duplicates(); |
| 191 |
|
|
void print_edges(); |
| 192 |
|
|
/** Returns a list of IDs of neighboring particles |
| 193 |
|
|
* corresponding to each face. |
| 194 |
|
|
* \param[out] v a reference to a vector in which to return the |
| 195 |
|
|
* results. If no neighbor information is |
| 196 |
|
|
* available, a blank vector is returned. */ |
| 197 |
|
|
virtual void neighbors(std::vector<int> &v) {v.clear();} |
| 198 |
|
|
/** This is a virtual function that is overridden by a routine |
| 199 |
|
|
* to print a list of IDs of neighboring particles |
| 200 |
|
|
* corresponding to each face. By default, when no neighbor |
| 201 |
|
|
* information is available, the routine does nothing. |
| 202 |
|
|
* \param[in] fp the file handle to write to. */ |
| 203 |
|
|
virtual void output_neighbors(FILE *fp=stdout) {} |
| 204 |
|
|
/** This a virtual function that is overridden by a routine to |
| 205 |
|
|
* print the neighboring particle IDs for a given vertex. By |
| 206 |
|
|
* default, when no neighbor information is available, the |
| 207 |
|
|
* routine does nothing. |
| 208 |
|
|
* \param[in] i the vertex to consider. */ |
| 209 |
|
|
virtual void print_edges_neighbors(int i) {}; |
| 210 |
|
|
/** This is a simple inline function for picking out the index |
| 211 |
|
|
* of the next edge counterclockwise at the current vertex. |
| 212 |
|
|
* \param[in] a the index of an edge of the current vertex. |
| 213 |
|
|
* \param[in] p the number of the vertex. |
| 214 |
|
|
* \return 0 if a=nu[p]-1, or a+1 otherwise. */ |
| 215 |
|
|
inline int cycle_up(int a,int p) {return a==nu[p]-1?0:a+1;} |
| 216 |
|
|
/** This is a simple inline function for picking out the index |
| 217 |
|
|
* of the next edge clockwise from the current vertex. |
| 218 |
|
|
* \param[in] a the index of an edge of the current vertex. |
| 219 |
|
|
* \param[in] p the number of the vertex. |
| 220 |
|
|
* \return nu[p]-1 if a=0, or a-1 otherwise. */ |
| 221 |
|
|
inline int cycle_down(int a,int p) {return a==0?nu[p]-1:a-1;} |
| 222 |
|
|
protected: |
| 223 |
|
|
/** This a one dimensional array that holds the current sizes |
| 224 |
|
|
* of the memory allocations for them mep array.*/ |
| 225 |
|
|
int *mem; |
| 226 |
|
|
/** This is a one dimensional array that holds the current |
| 227 |
|
|
* number of vertices of order p that are stored in the mep[p] |
| 228 |
|
|
* array. */ |
| 229 |
|
|
int *mec; |
| 230 |
|
|
/** This is a two dimensional array for holding the information |
| 231 |
|
|
* about the edges of the Voronoi cell. mep[p] is a |
| 232 |
|
|
* one-dimensional array for holding the edge information about |
| 233 |
|
|
* all vertices of order p, with each vertex holding 2*p+1 |
| 234 |
|
|
* integers of information. The total number of vertices held |
| 235 |
|
|
* on mep[p] is stored in mem[p]. If the space runs out, the |
| 236 |
|
|
* code allocates more using the add_memory() routine. */ |
| 237 |
|
|
int **mep; |
| 238 |
|
|
inline void reset_edges(); |
| 239 |
|
|
template<class vc_class> |
| 240 |
|
|
void check_memory_for_copy(vc_class &vc,voronoicell_base* vb); |
| 241 |
|
|
void copy(voronoicell_base* vb); |
| 242 |
|
|
private: |
| 243 |
|
|
/** This is the delete stack, used to store the vertices which |
| 244 |
|
|
* are going to be deleted during the plane cutting procedure. |
| 245 |
|
|
*/ |
| 246 |
|
|
int *ds,*stacke; |
| 247 |
|
|
/** This is the auxiliary delete stack, which has size set by |
| 248 |
|
|
* current_delete2_size. */ |
| 249 |
|
|
int *ds2,*stacke2; |
| 250 |
|
|
/** This stores the current memory allocation for the marginal |
| 251 |
|
|
* cases. */ |
| 252 |
|
|
int current_marginal; |
| 253 |
|
|
/** This stores the total number of marginal points which are |
| 254 |
|
|
* currently in the buffer. */ |
| 255 |
|
|
int n_marg; |
| 256 |
|
|
/** This array contains a list of the marginal points, and also |
| 257 |
|
|
* the outcomes of the marginal tests. */ |
| 258 |
|
|
int *marg; |
| 259 |
|
|
/** The x coordinate of the normal vector to the test plane. */ |
| 260 |
|
|
double px; |
| 261 |
|
|
/** The y coordinate of the normal vector to the test plane. */ |
| 262 |
|
|
double py; |
| 263 |
|
|
/** The z coordinate of the normal vector to the test plane. */ |
| 264 |
|
|
double pz; |
| 265 |
|
|
/** The magnitude of the normal vector to the test plane. */ |
| 266 |
|
|
double prsq; |
| 267 |
|
|
template<class vc_class> |
| 268 |
|
|
void add_memory(vc_class &vc,int i,int *stackp2); |
| 269 |
|
|
template<class vc_class> |
| 270 |
|
|
void add_memory_vertices(vc_class &vc); |
| 271 |
|
|
template<class vc_class> |
| 272 |
|
|
void add_memory_vorder(vc_class &vc); |
| 273 |
|
|
void add_memory_ds(int *&stackp); |
| 274 |
|
|
void add_memory_ds2(int *&stackp2); |
| 275 |
|
|
template<class vc_class> |
| 276 |
|
|
inline bool collapse_order1(vc_class &vc); |
| 277 |
|
|
template<class vc_class> |
| 278 |
|
|
inline bool collapse_order2(vc_class &vc); |
| 279 |
|
|
template<class vc_class> |
| 280 |
|
|
inline bool delete_connection(vc_class &vc,int j,int k,bool hand); |
| 281 |
|
|
template<class vc_class> |
| 282 |
|
|
inline bool search_for_outside_edge(vc_class &vc,int &up); |
| 283 |
|
|
template<class vc_class> |
| 284 |
|
|
inline void add_to_stack(vc_class &vc,int lp,int *&stackp2); |
| 285 |
|
|
inline bool plane_intersects_track(double x,double y,double z,double rs,double g); |
| 286 |
|
|
inline void normals_search(std::vector<double> &v,int i,int j,int k); |
| 287 |
|
|
inline bool search_edge(int l,int &m,int &k); |
| 288 |
|
|
inline int m_test(int n,double &ans); |
| 289 |
|
|
int check_marginal(int n,double &ans); |
| 290 |
|
|
friend class voronoicell; |
| 291 |
|
|
friend class voronoicell_neighbor; |
| 292 |
|
|
}; |
| 293 |
|
|
|
| 294 |
|
|
/** \brief Extension of the voronoicell_base class to represent a Voronoi |
| 295 |
|
|
* cell without neighbor information. |
| 296 |
|
|
* |
| 297 |
|
|
* This class is an extension of the voronoicell_base class, in cases when |
| 298 |
|
|
* is not necessary to track the IDs of neighboring particles associated |
| 299 |
|
|
* with each face of the Voronoi cell. */ |
| 300 |
|
|
class voronoicell : public voronoicell_base { |
| 301 |
|
|
public: |
| 302 |
|
|
using voronoicell_base::nplane; |
| 303 |
|
|
/** Copies the information from another voronoicell class into |
| 304 |
|
|
* this class, extending memory allocation if necessary. |
| 305 |
|
|
* \param[in] c the class to copy. */ |
| 306 |
|
|
inline void operator=(voronoicell &c) { |
| 307 |
|
|
voronoicell_base* vb((voronoicell_base*) &c); |
| 308 |
|
|
check_memory_for_copy(*this,vb);copy(vb); |
| 309 |
|
|
} |
| 310 |
|
|
/** Cuts a Voronoi cell using by the plane corresponding to the |
| 311 |
|
|
* perpendicular bisector of a particle. |
| 312 |
|
|
* \param[in] (x,y,z) the position of the particle. |
| 313 |
|
|
* \param[in] rsq the modulus squared of the vector. |
| 314 |
|
|
* \param[in] p_id the plane ID, ignored for this case where no |
| 315 |
|
|
* neighbor tracking is enabled. |
| 316 |
|
|
* \return False if the plane cut deleted the cell entirely, |
| 317 |
|
|
* true otherwise. */ |
| 318 |
|
|
inline bool nplane(double x,double y,double z,double rsq,int p_id) { |
| 319 |
|
|
return nplane(*this,x,y,z,rsq,0); |
| 320 |
|
|
} |
| 321 |
|
|
/** Cuts a Voronoi cell using by the plane corresponding to the |
| 322 |
|
|
* perpendicular bisector of a particle. |
| 323 |
|
|
* \param[in] (x,y,z) the position of the particle. |
| 324 |
|
|
* \param[in] p_id the plane ID, ignored for this case where no |
| 325 |
|
|
* neighbor tracking is enabled. |
| 326 |
|
|
* \return False if the plane cut deleted the cell entirely, |
| 327 |
|
|
* true otherwise. */ |
| 328 |
|
|
inline bool nplane(double x,double y,double z,int p_id) { |
| 329 |
|
|
double rsq=x*x+y*y+z*z; |
| 330 |
|
|
return nplane(*this,x,y,z,rsq,0); |
| 331 |
|
|
} |
| 332 |
|
|
/** Cuts a Voronoi cell using by the plane corresponding to the |
| 333 |
|
|
* perpendicular bisector of a particle. |
| 334 |
|
|
* \param[in] (x,y,z) the position of the particle. |
| 335 |
|
|
* \param[in] rsq the modulus squared of the vector. |
| 336 |
|
|
* \return False if the plane cut deleted the cell entirely, |
| 337 |
|
|
* true otherwise. */ |
| 338 |
|
|
inline bool plane(double x,double y,double z,double rsq) { |
| 339 |
|
|
return nplane(*this,x,y,z,rsq,0); |
| 340 |
|
|
} |
| 341 |
|
|
/** Cuts a Voronoi cell using by the plane corresponding to the |
| 342 |
|
|
* perpendicular bisector of a particle. |
| 343 |
|
|
* \param[in] (x,y,z) the position of the particle. |
| 344 |
|
|
* \return False if the plane cut deleted the cell entirely, |
| 345 |
|
|
* true otherwise. */ |
| 346 |
|
|
inline bool plane(double x,double y,double z) { |
| 347 |
|
|
double rsq=x*x+y*y+z*z; |
| 348 |
|
|
return nplane(*this,x,y,z,rsq,0); |
| 349 |
|
|
} |
| 350 |
|
|
/** Initializes the Voronoi cell to be rectangular box with the |
| 351 |
|
|
* given dimensions. |
| 352 |
|
|
* \param[in] (xmin,xmax) the minimum and maximum x coordinates. |
| 353 |
|
|
* \param[in] (ymin,ymax) the minimum and maximum y coordinates. |
| 354 |
|
|
* \param[in] (zmin,zmax) the minimum and maximum z coordinates. */ |
| 355 |
|
|
inline void init(double xmin,double xmax,double ymin,double ymax,double zmin,double zmax) { |
| 356 |
|
|
init_base(xmin,xmax,ymin,ymax,zmin,zmax); |
| 357 |
|
|
} |
| 358 |
|
|
/** Initializes the cell to be an octahedron with vertices at |
| 359 |
|
|
* (l,0,0), (-l,0,0), (0,l,0), (0,-l,0), (0,0,l), and (0,0,-l). |
| 360 |
|
|
* \param[in] l a parameter setting the size of the octahedron. |
| 361 |
|
|
*/ |
| 362 |
|
|
inline void init_octahedron(double l) { |
| 363 |
|
|
init_octahedron_base(l); |
| 364 |
|
|
} |
| 365 |
|
|
/** Initializes the cell to be a tetrahedron. |
| 366 |
|
|
* \param[in] (x0,y0,z0) the coordinates of the first vertex. |
| 367 |
|
|
* \param[in] (x1,y1,z1) the coordinates of the second vertex. |
| 368 |
|
|
* \param[in] (x2,y2,z2) the coordinates of the third vertex. |
| 369 |
|
|
* \param[in] (x3,y3,z3) the coordinates of the fourth vertex. |
| 370 |
|
|
*/ |
| 371 |
|
|
inline void init_tetrahedron(double x0,double y0,double z0,double x1,double y1,double z1,double x2,double y2,double z2,double x3,double y3,double z3) { |
| 372 |
|
|
init_tetrahedron_base(x0,y0,z0,x1,y1,z1,x2,y2,z2,x3,y3,z3); |
| 373 |
|
|
} |
| 374 |
|
|
private: |
| 375 |
|
|
inline void n_allocate(int i,int m) {}; |
| 376 |
|
|
inline void n_add_memory_vertices(int i) {}; |
| 377 |
|
|
inline void n_add_memory_vorder(int i) {}; |
| 378 |
|
|
inline void n_set_pointer(int p,int n) {}; |
| 379 |
|
|
inline void n_copy(int a,int b,int c,int d) {}; |
| 380 |
|
|
inline void n_set(int a,int b,int c) {}; |
| 381 |
|
|
inline void n_set_aux1(int k) {}; |
| 382 |
|
|
inline void n_copy_aux1(int a,int b) {}; |
| 383 |
|
|
inline void n_copy_aux1_shift(int a,int b) {}; |
| 384 |
|
|
inline void n_set_aux2_copy(int a,int b) {}; |
| 385 |
|
|
inline void n_copy_pointer(int a,int b) {}; |
| 386 |
|
|
inline void n_set_to_aux1(int j) {}; |
| 387 |
|
|
inline void n_set_to_aux2(int j) {}; |
| 388 |
|
|
inline void n_allocate_aux1(int i) {}; |
| 389 |
|
|
inline void n_switch_to_aux1(int i) {}; |
| 390 |
|
|
inline void n_copy_to_aux1(int i,int m) {}; |
| 391 |
|
|
inline void n_set_to_aux1_offset(int k,int m) {}; |
| 392 |
|
|
inline void n_neighbors(std::vector<int> &v) {v.clear();}; |
| 393 |
|
|
friend class voronoicell_base; |
| 394 |
|
|
}; |
| 395 |
|
|
|
| 396 |
|
|
/** \brief Extension of the voronoicell_base class to represent a Voronoi cell |
| 397 |
|
|
* with neighbor information. |
| 398 |
|
|
* |
| 399 |
|
|
* This class is an extension of the voronoicell_base class, in cases when the |
| 400 |
|
|
* IDs of neighboring particles associated with each face of the Voronoi cell. |
| 401 |
|
|
* It contains additional data structures mne and ne for storing this |
| 402 |
|
|
* information. */ |
| 403 |
|
|
class voronoicell_neighbor : public voronoicell_base { |
| 404 |
|
|
public: |
| 405 |
|
|
using voronoicell_base::nplane; |
| 406 |
|
|
/** This two dimensional array holds the neighbor information |
| 407 |
|
|
* associated with each vertex. mne[p] is a one dimensional |
| 408 |
|
|
* array which holds all of the neighbor information for |
| 409 |
|
|
* vertices of order p. */ |
| 410 |
|
|
int **mne; |
| 411 |
|
|
/** This is a two dimensional array that holds the neighbor |
| 412 |
|
|
* information associated with each vertex. ne[i] points to a |
| 413 |
|
|
* one-dimensional array in mne[nu[i]]. ne[i][j] holds the |
| 414 |
|
|
* neighbor information associated with the jth edge of vertex |
| 415 |
|
|
* i. It is set to the ID number of the plane that made the |
| 416 |
|
|
* face that is clockwise from the jth edge. */ |
| 417 |
|
|
int **ne; |
| 418 |
|
|
voronoicell_neighbor(); |
| 419 |
|
|
~voronoicell_neighbor(); |
| 420 |
|
|
void operator=(voronoicell &c); |
| 421 |
|
|
void operator=(voronoicell_neighbor &c); |
| 422 |
|
|
/** Cuts the Voronoi cell by a particle whose center is at a |
| 423 |
|
|
* separation of (x,y,z) from the cell center. The value of rsq |
| 424 |
|
|
* should be initially set to \f$x^2+y^2+z^2\f$. |
| 425 |
|
|
* \param[in] (x,y,z) the normal vector to the plane. |
| 426 |
|
|
* \param[in] rsq the distance along this vector of the plane. |
| 427 |
|
|
* \param[in] p_id the plane ID (for neighbor tracking only). |
| 428 |
|
|
* \return False if the plane cut deleted the cell entirely, |
| 429 |
|
|
* true otherwise. */ |
| 430 |
|
|
inline bool nplane(double x,double y,double z,double rsq,int p_id) { |
| 431 |
|
|
return nplane(*this,x,y,z,rsq,p_id); |
| 432 |
|
|
} |
| 433 |
|
|
/** This routine calculates the modulus squared of the vector |
| 434 |
|
|
* before passing it to the main nplane() routine with full |
| 435 |
|
|
* arguments. |
| 436 |
|
|
* \param[in] (x,y,z) the vector to cut the cell by. |
| 437 |
|
|
* \param[in] p_id the plane ID (for neighbor tracking only). |
| 438 |
|
|
* \return False if the plane cut deleted the cell entirely, |
| 439 |
|
|
* true otherwise. */ |
| 440 |
|
|
inline bool nplane(double x,double y,double z,int p_id) { |
| 441 |
|
|
double rsq=x*x+y*y+z*z; |
| 442 |
|
|
return nplane(*this,x,y,z,rsq,p_id); |
| 443 |
|
|
} |
| 444 |
|
|
/** This version of the plane routine just makes up the plane |
| 445 |
|
|
* ID to be zero. It will only be referenced if neighbor |
| 446 |
|
|
* tracking is enabled. |
| 447 |
|
|
* \param[in] (x,y,z) the vector to cut the cell by. |
| 448 |
|
|
* \param[in] rsq the modulus squared of the vector. |
| 449 |
|
|
* \return False if the plane cut deleted the cell entirely, |
| 450 |
|
|
* true otherwise. */ |
| 451 |
|
|
inline bool plane(double x,double y,double z,double rsq) { |
| 452 |
|
|
return nplane(*this,x,y,z,rsq,0); |
| 453 |
|
|
} |
| 454 |
|
|
/** Cuts a Voronoi cell using the influence of a particle at |
| 455 |
|
|
* (x,y,z), first calculating the modulus squared of this |
| 456 |
|
|
* vector before passing it to the main nplane() routine. Zero |
| 457 |
|
|
* is supplied as the plane ID, which will be ignored unless |
| 458 |
|
|
* neighbor tracking is enabled. |
| 459 |
|
|
* \param[in] (x,y,z) the vector to cut the cell by. |
| 460 |
|
|
* \return False if the plane cut deleted the cell entirely, |
| 461 |
|
|
* true otherwise. */ |
| 462 |
|
|
inline bool plane(double x,double y,double z) { |
| 463 |
|
|
double rsq=x*x+y*y+z*z; |
| 464 |
|
|
return nplane(*this,x,y,z,rsq,0); |
| 465 |
|
|
} |
| 466 |
|
|
void init(double xmin,double xmax,double ymin,double ymax,double zmin,double zmax); |
| 467 |
|
|
void init_octahedron(double l); |
| 468 |
|
|
void init_tetrahedron(double x0,double y0,double z0,double x1,double y1,double z1,double x2,double y2,double z2,double x3,double y3,double z3); |
| 469 |
|
|
void check_facets(); |
| 470 |
|
|
virtual void neighbors(std::vector<int> &v); |
| 471 |
|
|
virtual void print_edges_neighbors(int i); |
| 472 |
|
|
virtual void output_neighbors(FILE *fp=stdout) { |
| 473 |
|
|
std::vector<int> v;neighbors(v); |
| 474 |
|
|
voro_print_vector(v,fp); |
| 475 |
|
|
} |
| 476 |
|
|
private: |
| 477 |
|
|
int *paux1; |
| 478 |
|
|
int *paux2; |
| 479 |
|
|
inline void n_allocate(int i,int m) {mne[i]=new int[m*i];} |
| 480 |
|
|
inline void n_add_memory_vertices(int i) { |
| 481 |
|
|
int **pp=new int*[i]; |
| 482 |
|
|
for(int j=0;j<current_vertices;j++) pp[j]=ne[j]; |
| 483 |
|
|
delete [] ne;ne=pp; |
| 484 |
|
|
} |
| 485 |
|
|
inline void n_add_memory_vorder(int i) { |
| 486 |
|
|
int **p2=new int*[i]; |
| 487 |
|
|
for(int j=0;j<current_vertex_order;j++) p2[j]=mne[j]; |
| 488 |
|
|
delete [] mne;mne=p2; |
| 489 |
|
|
} |
| 490 |
|
|
inline void n_set_pointer(int p,int n) { |
| 491 |
|
|
ne[p]=mne[n]+n*mec[n]; |
| 492 |
|
|
} |
| 493 |
|
|
inline void n_copy(int a,int b,int c,int d) {ne[a][b]=ne[c][d];} |
| 494 |
|
|
inline void n_set(int a,int b,int c) {ne[a][b]=c;} |
| 495 |
|
|
inline void n_set_aux1(int k) {paux1=mne[k]+k*mec[k];} |
| 496 |
|
|
inline void n_copy_aux1(int a,int b) {paux1[b]=ne[a][b];} |
| 497 |
|
|
inline void n_copy_aux1_shift(int a,int b) {paux1[b]=ne[a][b+1];} |
| 498 |
|
|
inline void n_set_aux2_copy(int a,int b) { |
| 499 |
|
|
paux2=mne[b]+b*mec[b]; |
| 500 |
|
|
for(int i=0;i<b;i++) ne[a][i]=paux2[i]; |
| 501 |
|
|
} |
| 502 |
|
|
inline void n_copy_pointer(int a,int b) {ne[a]=ne[b];} |
| 503 |
|
|
inline void n_set_to_aux1(int j) {ne[j]=paux1;} |
| 504 |
|
|
inline void n_set_to_aux2(int j) {ne[j]=paux2;} |
| 505 |
|
|
inline void n_allocate_aux1(int i) {paux1=new int[i*mem[i]];} |
| 506 |
|
|
inline void n_switch_to_aux1(int i) {delete [] mne[i];mne[i]=paux1;} |
| 507 |
|
|
inline void n_copy_to_aux1(int i,int m) {paux1[m]=mne[i][m];} |
| 508 |
|
|
inline void n_set_to_aux1_offset(int k,int m) {ne[k]=paux1+m;} |
| 509 |
|
|
friend class voronoicell_base; |
| 510 |
|
|
}; |
| 511 |
|
|
|
| 512 |
|
|
} |
| 513 |
|
|
|
| 514 |
|
|
#endif |
