branches/dynamic/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     : July 1st 2008
   6
   7 /** \file cell.hh
   8  * \brief Header file for the voronoicell_base template and related classes. */
   9
  10 #ifndef VOROPP_CELL_HH
  11 #define VOROPP_CELL_HH
  12
  13 #include "config.hh"
  14 #include <cstdio>
  15 #include <iostream>
  16 #include <fstream>
  17 #include <cmath>
  18 using namespace std;
  19
  20 /** \brief Structure for printing fatal error messages and exiting.
  21  *
  22  * Structure for printing fatal error messages and exiting. */
  23 struct fatal_error {
  24         /** This routine prints an error message to the standard error.
  25          * \param[in] p The message to print. */
  26         fatal_error(const char *p) {cerr << p << endl;}
  27 };
  28
  29 /** \brief A class to reliably carry out floating point comparisons, storing
  30  * marginal cases for future reference.
  31  *
  32  * Floating point comparisons can be unreliable on some processor
  33  * architectures, and can produce unpredictable results. On a number of popular
  34  * Intel processors, floating point numbers are held to higher precision when
  35  * in registers than when in memory. When a register is swapped from a register
  36  * to memory, a truncation error, and in some situations this can create
  37  * circumstances where for two numbers c and d, the program finds c>d first,
  38  * but later c<d. The programmer has no control over when the swaps between
  39  * memory and registers occur, and recompiling with slightly different code can
  40  * give different results. One solution to avoid this is to force the compiler
  41  * to evaluate everything in memory (e.g. by using the -ffloat-store option in
  42  * the GNU C++ compiler) but this could be viewed overkill, since it slows the
  43  * code down, and the extra register precision is useful.
  44  *
  45  * In the plane cutting routine of the voronoicell class, we need to reliably
  46  * know whether a vertex lies inside, outside, or on the cutting plane, since
  47  * if it changed during the tracing process there would be confusion. This
  48  * class makes these tests reliable, by storing the results of marginal cases,
  49  * where the vertex lies within tolerance2 of the cutting plane. If that vertex
  50  * is tested again, then code looks up the value of the table in a buffer,
  51  * rather than doing the floating point comparison again. Only vertices which
  52  * are close to the plane are stored and tested, so this routine should create
  53  * minimal computational overhead.
  54  */
  55 class suretest { public:
  56                 /** This is a pointer to the array in the voronoicell class
  57                  * which holds the vertex coordinates.*/
  58                 fpoint *p;              suretest(); ~suretest(); inline void
  59                         init(fpoint x,fpoint y,fpoint z,fpoint rsq); inline int
  60                         test(int n,fpoint &ans); private: int
  61                         check_marginal(int n,fpoint &ans);
  62                 /** This stores the current memory allocation for the marginal
  63                  * cases. */
  64                 int current_marginal;
  65                 /** This stores the total number of marginal points which are
  66                  * currently in the buffer. */
  67                 int sc;
  68                 /** This array contains a list of the marginal points, and also
  69                  * the outcomes of the marginal tests. */
  70                 int *sn;
  71                 /** The x coordinate of the normal vector to the test plane. */
  72                 fpoint px;
  73                 /** The y coordinate of the normal vector to the test plane. */
  74                 fpoint py;
  75                 /** The z coordinate of the normal vector to the test plane. */
  76                 fpoint pz;
  77                 /** The magnitude of the normal vector to the test plane. */
  78                 fpoint prsq; };
  79
  80 class neighbor_track;
  81
  82 /** \brief A class encapsulating all the routines for storing and calculating
  83  * a single Voronoi cell.
  84  *
  85  * This class encapsulates all the routines for storing and calculating a
  86  * single Voronoi cell. The cell can first be initialized by the init() function
  87  * to be a rectangular box. The box can then be successively cut by planes
  88  * using the plane function. Other routines exist for outputting the cell,
  89  * computing its volume, or finding the largest distance of a vertex from the
  90  * cell center.  The cell is described by two arrays. pts[] is a floating point
  91  * array which holds the vertex positions. ed[] holds the table of edges, and
  92  * also a relation table that determines how two vertices are connected to one
  93  * another. The relation table is redundant, but helps speed up the
  94  * computation. The function check_relations() checks that the relational table
  95  * is valid. */
  96 template <class n_option>
  97 class voronoicell_base {
  98         public:
  99                 /** This is a two dimensional array that holds information about
 100                  * the edge connections of the vertices that make up the cell.
 101                  * The two dimensional array is not allocated in the usual method.
 102                  * To account for the fact the different vertices have different
 103                  * orders, and thus require different amounts of storage, the
 104                  * elements of ed[i] point to one-dimensional arrays in the mep[]
 105                  * array of different sizes.
 106                  *
 107                  * More specifically, if vertex i has order m, then ed[i]
 108                  * points to a one-dimensional array in mep[m] that has 2*m+1
 109                  * entries. The first m elements hold the neighboring edges, so
 110                  * that the jth edge of vertex i is held in ed[i][j]. The next
 111                  * m elements hold a table of relations which is redundant but
 112                  * helps speed up the computation. It satisfies the relation
 113                  * ed[ed[i][j]][ed[i][m+j]]=i. The final entry holds a back
 114                  * pointer, so that ed[i+2*m]=i. These are used when
 115                  * rearranging the memory. */
 116                 int **ed;
 117                 /** This array holds the order of the vertices in the Voronoi cell.
 118                  * This array is dynamically allocated, with its current size
 119                  * held by current_vertices. */
 120                 int *nu;
 121                 /** This holds the current size of the arrays ed and nu, which
 122                  * hold the vertex information. If more vertices are created
 123                  * than can fit in this array, then it is dynamically extended
 124                  * using the add_memory_vertices routine. */
 125                 int current_vertices;
 126                 /** This holds the current maximum allowed order of a vertex,
 127                  * which sets the size of the mem, mep, and mec arrays. If a
 128                  * vertex is created with more vertices than this, the arrays
 129                  * are dynamically extended using the add_memory_vorder routine.
 130                  */
 131                 int current_vertex_order;
 132                 /** This sets the size of the main delete stack. */
 133                 int current_delete_size;
 134                 /** This sets the size of the auxiliary delete stack. */
 135                 int current_delete2_size;
 136                 /** This in an array with size 3*current_vertices for holding
 137                  * the positions of the vertices. */
 138                 fpoint *pts;
 139                 /** This sets the total number of vertices in the current cell.
 140                  */
 141                 int p;
 142                 /** This is the index of particular point in the cell, which is used to start
 143                  * the tracing routines for plane intersection and cutting. These
 144                  * routines will work starting from any point, but it's often most
 145                  * efficient to start from the last point considered, since in many cases,
 146                  * the cell construction algorithm may consider many planes with similar
 147                  * vectors concurrently. */
 148                 int up;
 149                 /** This is a class used in the plane routine for carrying out
 150                  * reliable comparisons of whether points in the cell are
 151                  * inside, outside, or on the current cutting plane. */
 152                 suretest sure;
 153                 voronoicell_base();
 154                 ~voronoicell_base();
 155                 void init(fpoint xmin,fpoint xmax,fpoint ymin,fpoint ymax,fpoint zmin,fpoint zmax);
 156                 inline void init_octahedron(fpoint l);
 157                 inline void init_tetrahedron(fpoint x0,fpoint y0,fpoint z0,fpoint x1,fpoint y1,fpoint z1,fpoint x2,fpoint y2,fpoint z2,fpoint x3,fpoint y3,fpoint z3);
 158                 inline void init_test(int n);
 159                 inline void add_vertex(fpoint x,fpoint y,fpoint z,int a);
 160                 inline void add_vertex(fpoint x,fpoint y,fpoint z,int a,int b);
 161                 inline void add_vertex(fpoint x,fpoint y,fpoint z,int a,int b,int c);
 162                 inline void add_vertex(fpoint x,fpoint y,fpoint z,int a,int b,int c,int d);
 163                 inline void add_vertex(fpoint x,fpoint y,fpoint z,int a,int b,int c,int d,int e);
 164                 void draw_pov(ostream &os,fpoint x,fpoint y,fpoint z);
 165                 inline void draw_pov(const char *filename,fpoint x,fpoint y,fpoint z);
 166                 inline void draw_pov(fpoint x,fpoint y,fpoint z);
 167                 void draw_pov_mesh(ostream &os,fpoint x,fpoint y,fpoint z);
 168                 inline void draw_pov_mesh(const char *filename,fpoint x,fpoint y,fpoint z);
 169                 inline void draw_pov_mesh(fpoint x,fpoint y,fpoint z);
 170                 void draw_gnuplot(ostream &os,fpoint x,fpoint y,fpoint z);
 171                 inline void draw_gnuplot(const char *filename,fpoint x,fpoint y,fpoint z);
 172                 inline void draw_gnuplot(fpoint x,fpoint y,fpoint z);
 173                 inline void check_relations();
 174                 inline void check_duplicates();
 175                 inline void construct_relations();
 176                 fpoint volume();
 177                 fpoint maxradsq();
 178                 int number_of_faces();
 179                 void print_edges();
 180                 inline void perturb(fpoint r);
 181                 void facets(ostream &os);
 182                 inline void facets();
 183                 inline void facets(const char *filename);
 184                 void facet_statistics(ostream &os);
 185                 inline void facet_statistics();
 186                 inline void facet_statistics(const char *filename);
 187                 bool nplane(fpoint x,fpoint y,fpoint z,fpoint rs,int p_id);
 188                 inline bool nplane(fpoint x,fpoint y,fpoint z,int p_id);
 189                 inline bool plane(fpoint x,fpoint y,fpoint z,fpoint rs);
 190                 inline bool plane(fpoint x,fpoint y,fpoint z);
 191                 bool plane_intersects(fpoint x,fpoint y,fpoint z,fpoint rs);
 192                 bool plane_intersects_guess(fpoint x,fpoint y,fpoint z,fpoint rs);
 193                 void label_facets();
 194                 void neighbors(ostream &os);
 195                 void check_facets();
 196         private:
 197                 /** This a one dimensional array that holds the current sizes
 198                  * of the memory allocations for them mep array.*/
 199                 int *mem;
 200                 /** This is a two dimensional array for holding the information
 201                  * about the edges of the Voronoi cell. mep[p] is a
 202                  * one-dimensional array for holding the edge information about
 203                  * all vertices of order p, with each vertex holding 2*p+1
 204                  * integers of information. The total number of vertices held
 205                  * on mep[p] is stored in mem[p]. If the space runs out, the
 206                  * code allocates more using the add_memory() routine. */
 207                 int **mep;
 208                 /** This is a one dimensional array that holds the current
 209                  * number of vertices of order p that are stored in the mep[p]
 210                  * array. */
 211                 int *mec;
 212                 /** This is the delete stack, used to store the vertices which
 213                  * are going to be deleted during the plane cutting procedure.
 214                  */
 215                 int *ds;
 216                 /** This is the auxiliary delete stack, which has size set by
 217                  * current_delete2_size. */
 218                 int *ds2;
 219                 /** This holds the number of points currently on the auxiliary
 220                  * delete stack. */
 221                 int stack2;
 222                 /** This object contains all the functions required to carry
 223                  * out the neighbor computation. If the neighbor_none class is
 224                  * used for n_option, then all these functions are blank. If
 225                  * the neighbor_track class is used, then the neighbor tracking
 226                  * is enabled. All the functions for the n_option classes are
 227                  * declared inline, so that they should all be completely
 228                  * integrated into the routine during compilation. */
 229                 n_option neighbor;
 230                 inline int cycle_up(int a,int p);
 231                 inline int cycle_down(int a,int p);
 232                 void add_memory(int i);
 233                 void add_memory_vertices();
 234                 void add_memory_vorder();
 235                 void add_memory_ds();
 236                 void add_memory_ds2();
 237                 inline bool collapse_order1();
 238                 inline bool collapse_order2();
 239                 inline bool delete_connection(int j,int k,bool hand);
 240                 inline bool plane_intersects_track(fpoint x,fpoint y,fpoint z,fpoint rs,fpoint g);
 241                 friend class neighbor_track;
 242 };
 243
 244 /** \brief A class passed to the voronoicell_base template to switch off
 245  * neighbor computation.
 246  *
 247  * This is a class full of empty routines for neighbor computation. If the
 248  * voronoicell_base template is instantiated with this class, then it
 249  * has the effect of switching off all neighbor computation. Since all these
 250  * routines are declared inline, it should have the effect of a zero speed
 251  * overhead in the resulting code. */
 252 class neighbor_none {
 253         public:
 254                 /** This is a blank constructor. */
 255                 neighbor_none(voronoicell_base<neighbor_none> *ivc) {};
 256                 /** This is a blank placeholder function that does nothing. */
 257                 inline void allocate(int i,int m) {};
 258                 /** This is a blank placeholder function that does nothing. */
 259                 inline void add_memory_vertices(int i) {};
 260                 /** This is a blank placeholder function that does nothing. */
 261                 inline void add_memory_vorder(int i) {};
 262                 /** This is a blank placeholder function that does nothing. */
 263                 inline void init() {};
 264                 /** This is a blank placeholder function that does nothing. */
 265                 inline void init_octahedron() {};
 266                 /** This is a blank placeholder function that does nothing. */
 267                 inline void init_tetrahedron() {};
 268                 /** This is a blank placeholder function that does nothing. */
 269                 inline void set_pointer(int p,int n) {};
 270                 /** This is a blank placeholder function that does nothing. */
 271                 inline void copy(int a,int b,int c,int d) {};
 272                 /** This is a blank placeholder function that does nothing. */
 273                 inline void set(int a,int b,int c) {};
 274                 /** This is a blank placeholder function that does nothing. */
 275                 inline void set_aux1(int k) {};
 276                 /** This is a blank placeholder function that does nothing. */
 277                 inline void copy_aux1(int a,int b) {};
 278                 /** This is a blank placeholder function that does nothing. */
 279                 inline void copy_aux1_shift(int a,int b) {};
 280                 /** This is a blank placeholder function that does nothing. */
 281                 inline void set_aux2_copy(int a,int b) {};
 282                 /** This is a blank placeholder function that does nothing. */
 283                 inline void copy_pointer(int a,int b) {};
 284                 /** This is a blank placeholder function that does nothing. */
 285                 inline void set_to_aux1(int j) {};
 286                 /** This is a blank placeholder function that does nothing. */
 287                 inline void set_to_aux2(int j) {};
 288                 /** This is a blank placeholder function that does nothing. */
 289                 inline void print_edges(int i) {};
 290                 /** This is a blank placeholder function that does nothing. */
 291                 inline void allocate_aux1(int i) {};
 292                 /** This is a blank placeholder function that does nothing. */
 293                 inline void switch_to_aux1(int i) {};
 294                 /** This is a blank placeholder function that does nothing. */
 295                 inline void copy_to_aux1(int i,int m) {};
 296                 /** This is a blank placeholder function that does nothing. */
 297                 inline void set_to_aux1_offset(int k,int m) {};
 298                 /** This is a blank placeholder function that does nothing. */
 299                 inline void print(ostream &os,int i,int j);
 300                 /** This is a blank placeholder function that does nothing. */
 301                 inline void label_facets() {};
 302                 /** This is a blank placeholder function that does nothing. */
 303                 inline void neighbors(ostream &os) {};
 304                 /** This is a blank placeholder function that does nothing. */
 305                 inline void check_facets() {};
 306                 inline bool internal_wall(int i) {return true;}
 307                 inline bool internal_wall(int i,int j) {return true;}
 308 };
 309
 310 /** \brief A class passed to the voronoicell_base template to switch on the
 311  * neighbor computation.
 312  *
 313  * This class encapsulates all the routines which are required to carry out the
 314  * neighbor tracking. If the voronoicell_base template is instantiated with
 315  * this class, then the neighbor computation is enabled. All these routines are
 316  * simple and declared inline, so they should be directly integrated into the
 317  * functions in the voronoicell class during compilation, without zero function
 318  * call overhead. */
 319 class neighbor_track {
 320         public:
 321                 /** This two dimensional array holds the neighbor information
 322                  * associated with each vertex. mne[p] is a one dimensional
 323                  * array which holds all of the neighbor information for
 324                  * vertices of order p. */
 325                 int **mne;
 326                 /** This is a two dimensional array that holds the neighbor
 327                  * information associated with each vertex. ne[i] points to a
 328                  * one-dimensional array in mne[nu[i]]. ne[i][j] holds the
 329                  * neighbor information associated with the jth edge of vertex
 330                  * i. It is set to the ID number of the plane that made the
 331                  * face that is clockwise from the jth edge. */
 332                 int **ne;
 333                 neighbor_track(voronoicell_base<neighbor_track> *ivc);
 334                 ~neighbor_track();
 335                 /** This is a pointer back to the voronoicell class which
 336                  * created this class. It is used to reference the members of
 337                  * that class in computations. */
 338                 voronoicell_base<neighbor_track> *vc;
 339                 inline void allocate(int i,int m);
 340                 inline void add_memory_vertices(int i);
 341                 inline void add_memory_vorder(int i);
 342                 inline void init();
 343                 inline void init_octahedron();
 344                 inline void init_tetrahedron();
 345                 inline void set_pointer(int p,int n);
 346                 inline void copy(int a,int b,int c,int d);
 347                 inline void set(int a,int b,int c);
 348                 inline void set_aux1(int k);
 349                 inline void copy_aux1(int a,int b);
 350                 inline void copy_aux1_shift(int a,int b);
 351                 inline void set_aux2_copy(int a,int b);
 352                 inline void copy_pointer(int a,int b);
 353                 inline void set_to_aux1(int j);
 354                 inline void set_to_aux2(int j);
 355                 inline void print_edges(int i);
 356                 inline void allocate_aux1(int i);
 357                 inline void switch_to_aux1(int i);
 358                 inline void copy_to_aux1(int i,int m);
 359                 inline void set_to_aux1_offset(int k,int m);
 360                 inline void print(ostream &os,int i,int j);
 361                 inline void label_facets();
 362                 inline void neighbors(ostream &os);
 363                 inline void check_facets();
 364                 inline bool internal_wall(int i);
 365                 inline bool internal_wall(int i,int j);
 366         private:
 367                 /** This is an auxiliary pointer which is used in some of the
 368                  * low level neighbor operations. */
 369                 int *paux1;
 370                 /** This is a second auxiliary pointer which is used in some
 371                  * of the low level neighbor operations. */
 372                 int *paux2;
 373 };
 374
 375 /** The basic voronoicell class. */
 376 typedef voronoicell_base<neighbor_none> voronoicell;
 377
 378 /** A neighbor-tracking version of the voronoicell. */
 379 typedef voronoicell_base<neighbor_track> voronoicell_neighbor;
 380 #endif