trunk/src/cell.hh

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