branches/dynamic/src/cell.cc

   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.cc
   8  * \brief Function implementations for the voronoicell_base template and
   9  * related classes. */
  10
  11 #include "cell.hh"
  12
  13 /** Constructs a Voronoi cell and sets up the initial memory. */
  14 template<class n_option>
  15 voronoicell_base<n_option>::voronoicell_base() :
  16         current_vertices(init_vertices), current_vertex_order(init_vertex_order),
  17         current_delete_size(init_delete_size), current_delete2_size(init_delete2_size),
  18         neighbor(this) {
  19         int i;
  20         ds=new int[current_delete_size];
  21         ds2=new int[current_delete2_size];
  22         mem=new int[current_vertex_order];
  23         mec=new int[current_vertex_order];
  24         mep=new int*[current_vertex_order];
  25         ed=new int*[current_vertices];
  26         nu=new int[current_vertices];
  27         pts=new fpoint[3*current_vertices];
  28         sure.p=pts;
  29         for(i=0;i<3;i++) {
  30                 mem[i]=init_n_vertices;
  31                 mep[i]=new int[init_n_vertices*(2*i+1)];
  32                 mec[i]=0;
  33         }
  34         mem[3]=init_3_vertices;
  35         mep[3]=new int[init_3_vertices*7];
  36         mec[3]=0;
  37         for(i=4;i<current_vertex_order;i++) {
  38                 mem[i]=init_n_vertices;
  39                 mep[i]=new int[init_n_vertices*(2*i+1)];
  40                 mec[i]=0;
  41         }
  42 }
  43
  44 /** The voronoicell destructor deallocates all the dynamic memory. */
  45 template<class n_option>
  46 voronoicell_base<n_option>::~voronoicell_base() {
  47         delete [] ds;
  48         delete [] ds2;
  49         for(int i=0;i<current_vertex_order;i++) if (mem[i]>0) delete [] mep[i];
  50         delete [] mem;
  51         delete [] mec;
  52         delete [] mep;
  53         delete [] ed;
  54         delete [] nu;
  55         delete [] pts;
  56 }
  57
  58 /** Increases the memory storage for a particular vertex order, by increasing
  59  * the size of the of the mep and mne <em>(neighbor option only)</em> arrays.
  60  * If the arrays already exist, their size is doubled; if they don't exist,
  61  * then new ones of size init_n_vertices are allocated. The routine also ensures
  62  * that the pointers in the ed array are updated, by making use of the back
  63  * pointers. For the cases where the back pointer has been temporarily
  64  * overwritten in the marginal vertex code, the auxiliary delete stack is
  65  * scanned to find out how to update the ed value.
  66  * \param[in] i The order of the vertex memory to be increased.*/
  67 template<class n_option>
  68 void voronoicell_base<n_option>::add_memory(int i) {
  69         int s=2*i+1;
  70         if(mem[i]==0) {
  71                 neighbor.allocate(i,init_n_vertices);
  72                 mep[i]=new int[init_n_vertices*s];
  73                 mem[i]=init_n_vertices;
  74 #if VOROPP_VERBOSE >=2
  75                 cerr << "Order " << i << " vertex memory created " << endl;
  76 #endif
  77         } else {
  78                 int j=0,k,*l;
  79                 mem[i]*=2;
  80                 if (mem[i]>max_n_vertices) throw fatal_error("Point memory allocation exceeded absolute maximum");
  81 #if VOROPP_VERBOSE >=2
  82                 cerr << "Order " << i << " vertex memory scaled up to " << mem[i] << endl;
  83 #endif
  84                 l=new int[s*mem[i]];
  85                 int m=0;
  86                 neighbor.allocate_aux1(i);
  87                 while(j<s*mec[i]) {
  88                         k=mep[i][j+2*i];
  89                         if(k>=0) {
  90                                 ed[k]=l+j;
  91                                 neighbor.set_to_aux1_offset(k,m);
  92                         } else {
  93                                 int o;
  94                                 for(o=0;o<stack2;o++) {
  95                                         if(ed[ds2[o]]==mep[i]+j) {
  96                                                 ed[ds2[o]]=l+j;
  97                                                 neighbor.set_to_aux1_offset(ds2[o],m);
  98                                                 break;
  99                                         }
 100                                 }
 101                                 if(o==stack2) throw fatal_error("Couldn't relocate dangling pointer");
 102 #if VOROPP_VERBOSE >=3
 103                                 cerr << "Relocated dangling pointer" << endl;
 104 #endif
 105                         }
 106                         for(k=0;k<s;k++,j++) l[j]=mep[i][j];
 107                         for(k=0;k<i;k++,m++) neighbor.copy_to_aux1(i,m);
 108                 }
 109                 delete [] mep[i];
 110                 mep[i]=l;
 111                 neighbor.switch_to_aux1(i);
 112         }
 113 }
 114
 115 /** Doubles the maximum number of vertices allowed, by reallocating the ed, nu,
 116  * pts, and ne <em>(neighbor option only)</em> arrays. If the allocation
 117  * exceeds the absolute maximum set in max_vertices, then the routine throws a
 118  * fatal error. */
 119 template<class n_option>
 120 void voronoicell_base<n_option>::add_memory_vertices() {
 121         int i=2*current_vertices,j,**pp,*pnu;
 122         if (i>max_vertices) throw fatal_error("Vertex memory allocation exceeded absolute maximum");
 123 #if VOROPP_VERBOSE >=2
 124         cerr << "Vertex memory scaled up to " << i << endl;
 125 #endif
 126         fpoint *ppts;
 127         pp=new int*[i];
 128         for(j=0;j<current_vertices;j++) pp[j]=ed[j];
 129         delete [] ed;ed=pp;
 130         neighbor.add_memory_vertices(i);
 131         pnu=new int[i];
 132         for(j=0;j<current_vertices;j++) pnu[j]=nu[j];
 133         delete [] nu;nu=pnu;
 134         ppts=new fpoint[3*i];
 135         for(j=0;j<3*current_vertices;j++) ppts[j]=pts[j];
 136         delete [] pts;sure.p=pts=ppts;
 137         current_vertices=i;
 138 }
 139
 140 /** Doubles the maximum allowed vertex order, by reallocating mem, mep, mec,
 141  * and mnu <em>(neighbor option only)</em> arrays. If the allocation exceeds
 142  * the absolute maximum set in max_vertex_order, then the routine causes a fatal
 143  * error. */
 144 template<class n_option>
 145 void voronoicell_base<n_option>::add_memory_vorder() {
 146         int i=2*current_vertex_order,j,*p1,**p2;
 147         if (i>max_vertex_order) throw fatal_error("Vertex order memory allocation exceeded absolute maximum");
 148 #if VOROPP_VERBOSE >=2
 149         cerr << "Vertex order memory scaled up to " << i << endl;
 150 #endif
 151         p1=new int[i];
 152         for(j=0;j<current_vertex_order;j++) p1[j]=mem[j];while(j<i) p1[j++]=0;
 153         delete [] mem;mem=p1;
 154         p2=new int*[i];
 155         for(j=0;j<current_vertex_order;j++) p2[j]=mep[j];
 156         delete [] mep;mep=p2;
 157         p1=new int[i];
 158         for(j=0;j<current_vertex_order;j++) p1[j]=mec[j];while(j<i) p1[j++]=0;
 159         delete [] mec;mec=p1;
 160         neighbor.add_memory_vorder(i);
 161         current_vertex_order=i;
 162 }
 163
 164 /** Doubles the size allocation of the main delete stack. If the allocation
 165  * exceeds the absolute maximum set in max_delete_size, then routine causes a
 166  * fatal error. */
 167 template<class n_option>
 168 void voronoicell_base<n_option>::add_memory_ds() {
 169         int i=2*current_delete_size,j,*pds;
 170         if (i>max_delete_size) throw fatal_error("Delete stack 1 memory allocation exceeded absolute maximum");
 171 #if VOROPP_VERBOSE >=2
 172         cerr << "Delete stack 1 memory scaled up to " << i << endl;
 173 #endif
 174         pds=new int[i];
 175         for(j=0;j<current_delete_size;j++) pds[j]=ds[j];
 176         delete [] ds;ds=pds;
 177         current_delete_size=i;
 178 }
 179
 180 /** Doubles the size allocation of the auxiliary delete stack. If the
 181  * allocation exceeds the absolute maximum set in max_delete2_size, then the
 182  * routine causes a fatal error. */
 183 template<class n_option>
 184 void voronoicell_base<n_option>::add_memory_ds2() {
 185         int i=2*current_delete2_size,j,*pds2;
 186         if (i>max_delete2_size) throw fatal_error("Delete stack 2 memory allocation exceeded absolute maximum");
 187 #if VOROPP_VERBOSE >=2
 188         cerr << "Delete stack 2 memory scaled up to " << i << endl;
 189 #endif
 190         pds2=new int[i];
 191         for(j=0;j<current_delete2_size;j++) pds2[j]=ds2[j];
 192         delete [] ds2;ds2=pds2;
 193         current_delete2_size=i;
 194 }
 195
 196 /** Initializes a Voronoi cell as a rectangular box with the given dimensions */
 197 template<class n_option>
 198 void voronoicell_base<n_option>::init(fpoint xmin,fpoint xmax,fpoint ymin,fpoint ymax,fpoint zmin,fpoint zmax) {
 199         for(int i=0;i<current_vertex_order;i++) mec[i]=0;up=0;
 200         mec[3]=p=8;xmin*=2;xmax*=2;ymin*=2;ymax*=2;zmin*=2;zmax*=2;
 201         pts[0]=xmin;pts[1]=ymin;pts[2]=zmin;
 202         pts[3]=xmax;pts[4]=ymin;pts[5]=zmin;
 203         pts[6]=xmin;pts[7]=ymax;pts[8]=zmin;
 204         pts[9]=xmax;pts[10]=ymax;pts[11]=zmin;
 205         pts[12]=xmin;pts[13]=ymin;pts[14]=zmax;
 206         pts[15]=xmax;pts[16]=ymin;pts[17]=zmax;
 207         pts[18]=xmin;pts[19]=ymax;pts[20]=zmax;
 208         pts[21]=xmax;pts[22]=ymax;pts[23]=zmax;
 209         int *q=mep[3];
 210         q[0]=1;q[1]=4;q[2]=2;q[3]=2;q[4]=1;q[5]=0;q[6]=0;
 211         q[7]=3;q[8]=5;q[9]=0;q[10]=2;q[11]=1;q[12]=0;q[13]=1;
 212         q[14]=0;q[15]=6;q[16]=3;q[17]=2;q[18]=1;q[19]=0;q[20]=2;
 213         q[21]=2;q[22]=7;q[23]=1;q[24]=2;q[25]=1;q[26]=0;q[27]=3;
 214         q[28]=6;q[29]=0;q[30]=5;q[31]=2;q[32]=1;q[33]=0;q[34]=4;
 215         q[35]=4;q[36]=1;q[37]=7;q[38]=2;q[39]=1;q[40]=0;q[41]=5;
 216         q[42]=7;q[43]=2;q[44]=4;q[45]=2;q[46]=1;q[47]=0;q[48]=6;
 217         q[49]=5;q[50]=3;q[51]=6;q[52]=2;q[53]=1;q[54]=0;q[55]=7;
 218         ed[0]=q;ed[1]=q+7;ed[2]=q+14;ed[3]=q+21;
 219         ed[4]=q+28;ed[5]=q+35;ed[6]=q+42;ed[7]=q+49;
 220         neighbor.init();
 221         nu[0]=nu[1]=nu[2]=nu[3]=nu[4]=nu[5]=nu[6]=nu[7]=3;
 222 }
 223
 224 /** Initializes a Voronoi cell as a regular octahedron.
 225  * \param[in] l The distance from the octahedron center to a vertex. Six vertices
 226  * are initialized at (-l,0,0), (l,0,0), (0,-l,0), (0,l,0), (0,0,-l), and (0,0,l).*/
 227 template<class n_option>
 228 inline void voronoicell_base<n_option>::init_octahedron(fpoint l) {
 229         for(int i=0;i<current_vertex_order;i++) mec[i]=0;up=0;
 230         mec[4]=p=6;l*=2;
 231         pts[0]=-l;pts[1]=0;pts[2]=0;
 232         pts[3]=l;pts[4]=0;pts[5]=0;
 233         pts[6]=0;pts[7]=-l;pts[8]=0;
 234         pts[9]=0;pts[10]=l;pts[11]=0;
 235         pts[12]=0;pts[13]=0;pts[14]=-l;
 236         pts[15]=0;pts[16]=0;pts[17]=l;
 237         int *q=mep[4];
 238         q[0]=2;q[1]=5;q[2]=3;q[3]=4;q[4]=0;q[5]=0;q[6]=0;q[7]=0;q[8]=0;
 239         q[9]=2;q[10]=4;q[11]=3;q[12]=5;q[13]=2;q[14]=2;q[15]=2;q[16]=2;q[17]=1;
 240         q[18]=0;q[19]=4;q[20]=1;q[21]=5;q[22]=0;q[23]=3;q[24]=0;q[25]=1;q[26]=2;
 241         q[27]=0;q[28]=5;q[29]=1;q[30]=4;q[31]=2;q[32]=3;q[33]=2;q[34]=1;q[35]=3;
 242         q[36]=0;q[37]=3;q[38]=1;q[39]=2;q[40]=3;q[41]=3;q[42]=1;q[43]=1;q[44]=4;
 243         q[45]=0;q[46]=2;q[47]=1;q[48]=3;q[49]=1;q[50]=3;q[51]=3;q[52]=1;q[53]=5;
 244         ed[0]=q;ed[1]=q+9;ed[2]=q+18;ed[3]=q+27;ed[4]=q+36;ed[5]=q+45;
 245         neighbor.init_octahedron();
 246         nu[0]=nu[1]=nu[2]=nu[3]=nu[4]=nu[5]=4;
 247 }
 248
 249 /** Initializes a Voronoi cell as a tetrahedron. It assumes that the normal to
 250  * the face for the first three vertices points inside.
 251  * \param (x0,y0,z0) A position vector for the first vertex.
 252  * \param (x1,y1,z1) A position vector for the second vertex.
 253  * \param (x2,y2,z2) A position vector for the third vertex.
 254  * \param (x3,y3,z3) A position vector for the fourth vertex. */
 255 template<class n_option>
 256 inline void voronoicell_base<n_option>::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) {
 257         for(int i=0;i<current_vertex_order;i++) mec[i]=0;up=0;
 258         mec[3]=p=4;
 259         pts[0]=x0*2;pts[1]=y0*2;pts[2]=z0*2;
 260         pts[3]=x1*2;pts[4]=y1*2;pts[5]=z1*2;
 261         pts[6]=x2*2;pts[7]=y2*2;pts[8]=z2*2;
 262         pts[9]=x3*2;pts[10]=y3*2;pts[11]=z3*2;
 263         int *q=mep[3];
 264         q[0]=1;q[1]=3;q[2]=2;q[3]=0;q[4]=0;q[5]=0;q[6]=0;
 265         q[7]=0;q[8]=2;q[9]=3;q[10]=0;q[11]=2;q[12]=1;q[13]=1;
 266         q[14]=0;q[15]=3;q[16]=1;q[17]=2;q[18]=2;q[19]=1;q[20]=2;
 267         q[21]=0;q[22]=1;q[23]=2;q[24]=1;q[25]=2;q[26]=1;q[27]=3;
 268         ed[0]=q;ed[1]=q+7;ed[2]=q+14;ed[3]=q+21;
 269         neighbor.init_tetrahedron();
 270         nu[0]=nu[1]=nu[2]=nu[3]=3;
 271 }
 272
 273 /** Initializes an arbitrary test object using the add_vertex() and
 274  * construct_relations() routines. See the source code for information about
 275  * the specific objects.
 276  * \param[in] n the number of the test object (from 0 to 9)
 277  */
 278 template<class n_option>
 279 inline void voronoicell_base<n_option>::init_test(int n) {
 280         for(int i=0;i<current_vertex_order;i++) mec[i]=0;up=p=0;
 281         switch(n) {
 282                 case 0:
 283                         // A peaked object, with a high vertex 6, and a ridge
 284                         // at z=0 from vertex 1 to 2. This can be used to test
 285                         // the order 1 vertex collapse routine.
 286                         add_vertex(1,-2,-1,3,1,5);
 287                         add_vertex(0,-1,1,5,0,2);
 288                         add_vertex(0,1,0,1,3,6,4);
 289                         add_vertex(1,4,-1,4,6,2,0);
 290                         add_vertex(-1,4,-1,3,5,2,6);
 291                         add_vertex(-1,-2,-1,4,0,1);
 292                         add_vertex(0,3,0,4,2,3);
 293                         break;
 294                 case 1:
 295                         // A truncated pyramid shape, with vertex 4 in the z=0
 296                         // plane. This can be used to test order 4 vertex
 297                         // generation.
 298                         add_vertex(-2,2,-1,3,4,1);
 299                         add_vertex(2,2,-1,0,5,2);
 300                         add_vertex(2,-2,-1,1,6,3);
 301                         add_vertex(-2,-2,-1,2,7,0);
 302                         add_vertex(-1,1,0,0,7,5);
 303                         add_vertex(1,1,1,1,4,6);
 304                         add_vertex(1,-1,1,7,2,5);
 305                         add_vertex(-1,-1,1,4,3,6);
 306                         break;
 307                 case 2:
 308                         // An object with two peaks at vertices 4 and 6,
 309                         // connected with a trough at vertex 5. It can be used
 310                         // to test the part of the routine that deals with
 311                         // augmenting existing vertices.
 312                         add_vertex(1,-2,-1,1,3,4);
 313                         add_vertex(-1,-2,-1,2,0,4,5);
 314                         add_vertex(-1,2,-1,1,6,3);
 315                         add_vertex(1,2,-1,2,6,5,0);
 316                         add_vertex(0,-1,1,5,1,0);
 317                         add_vertex(0,0,0,1,4,3,6);
 318                         add_vertex(0,1,1,2,5,3);
 319                         break;
 320                 case 3:
 321                         // A box with a pyramid on top of it. This a good
 322                         // test object to make sure that the code can handle
 323                         // object with vertices of different orders.
 324                         add_vertex(-1,-1,-1,4,3,1);
 325                         add_vertex(1,-1,-1,0,2,5);
 326                         add_vertex(1,1,-1,6,1,3);
 327                         add_vertex(-1,1,-1,2,0,7);
 328                         add_vertex(-1,-1,1,7,0,5,8);
 329                         add_vertex(1,-1,1,4,1,6,8);
 330                         add_vertex(1,1,1,5,2,7,8);
 331                         add_vertex(-1,1,1,6,3,4,8);
 332                         add_vertex(0,0,2,7,4,5,6);
 333                         break;
 334                 case 4:
 335                         // A shape with two peaks (vertices 6 and 9) connected
 336                         // with a trough (vertices 7 and 8). It can be used as
 337                         // a basic test of the double edge skipping in the plane
 338                         // generation routine, whereby an edge is omitted if
 339                         // the routine is tracing along a part it has already
 340                         // been down.
 341                         add_vertex(1,-3,-1,5,6,1);
 342                         add_vertex(-1,-3,-1,0,6,2);
 343                         add_vertex(-3,0,-1,1,7,8,3);
 344                         add_vertex(-1,3,-1,2,9,4);
 345                         add_vertex(1,3,-1,3,9,5);
 346                         add_vertex(3,0,-1,4,8,7,0);
 347                         add_vertex(0,-2,1,7,1,0);
 348                         add_vertex(0,-1,0,8,2,6,5);
 349                         add_vertex(0,1,0,9,2,7,5);
 350                         add_vertex(0,2,1,3,8,4);
 351                         break;
 352                 case 5:
 353                         // An object with four peaks (vertices 8 to 11)
 354                         // connected by a trough (vertex 16). It can be used to
 355                         // the test multiple augmentation of edges of a
 356                         // particular vertex.
 357                         add_vertex(-1,-3,-1,7,1,8,12);
 358                         add_vertex(1,-3,-1,2,12,8,0);
 359                         add_vertex(3,-1,-1,3,9,13,1);
 360                         add_vertex(3,1,-1,4,13,9,2);
 361                         add_vertex(1,3,-1,5,10,14,3);
 362                         add_vertex(-1,3,-1,6,14,10,4);
 363                         add_vertex(-3,1,-1,7,11,15,5);
 364                         add_vertex(-3,-1,-1,0,15,11,6);
 365                         add_vertex(0,-2,1,0,1,12);
 366                         add_vertex(2,0,1,2,3,13);
 367                         add_vertex(0,2,1,14,4,5);
 368                         add_vertex(-2,0,1,7,15,6);
 369                         add_vertex(0,-1,0.5,0,8,1,16);
 370                         add_vertex(1,0,0.5,2,9,3,16);
 371                         add_vertex(0,1,0.5,4,10,5,16);
 372                         add_vertex(-1,0,0.5,6,11,7,16);
 373                         add_vertex(0,0,0,15,12,13,14);
 374                         break;
 375                 case 6:
 376                         // An object with four peaks (vertices 8 to 11)
 377                         // connected by a sequence of ridges. It can be used to
 378                         // the test multiple cases of double edge skipping
 379                         // on the same vertex.
 380                         add_vertex(-1,-3,-1,7,1,8,12);
 381                         add_vertex(1,-3,-1,2,12,8,0);
 382                         add_vertex(3,-1,-1,3,9,13,1);
 383                         add_vertex(3,1,-1,4,13,9,2);
 384                         add_vertex(1,3,-1,5,10,14,3);
 385                         add_vertex(-1,3,-1,6,14,10,4);
 386                         add_vertex(-3,1,-1,7,11,15,5);
 387                         add_vertex(-3,-1,-1,0,15,11,6);
 388                         add_vertex(0,-2,1,0,1,12);
 389                         add_vertex(2,0,1,2,3,13);
 390                         add_vertex(0,2,1,14,4,5);
 391                         add_vertex(-2,0,1,7,15,6);
 392                         add_vertex(0,-1,0,0,8,1,16);
 393                         add_vertex(1,0,0,2,9,3,16);
 394                         add_vertex(0,1,0,4,10,5,16);
 395                         add_vertex(-1,0,0,6,11,7,16);
 396                         add_vertex(0,0,0,15,12,13,14);
 397                         break;
 398                 case 7:
 399                         // A variation on the zeroth test shape, with a peak
 400                         // (vertices 7 and 8) connected to a ridge at z=0
 401                         // (vertices 4 and 5)
 402                         add_vertex(2,-3,-1,3,4,1);
 403                         add_vertex(-2,-3,-1,0,4,2);
 404                         add_vertex(-2,3,-1,1,7,3);
 405                         add_vertex(2,3,-1,2,6,0);
 406                         add_vertex(0,-2,0,0,5,1);
 407                         add_vertex(0,1,0,4,6,7);
 408                         add_vertex(1,2,1,7,5,3);
 409                         add_vertex(-1,2,1,5,6,2);
 410                         break;
 411                 case 8:
 412                         // A triangular object with a skewed peak, that can be used to
 413                         // test the order two vertex removal routine
 414                         add_vertex(3,-2,-1,2,3,1);
 415                         add_vertex(-3,-2,-1,0,4,2);
 416                         add_vertex(0,4,-1,1,5,0);
 417                         add_vertex(1.5,-1,0,7,0,6);
 418                         add_vertex(-1.5,-1,0,1,7,8);
 419                         add_vertex(0,2,0,8,6,2);
 420                         add_vertex(0.75,0.5,0,9,3,5);
 421                         add_vertex(0,-1,0,4,3,9);
 422                         add_vertex(-0.75,0.5,0,4,9,5);
 423                         add_vertex(0,0,1,8,7,6);
 424                         break;
 425                 case 9:
 426                         // This a tetrahedron with some low-order extraneous edges, and can
 427                         // be used to test the order 1 and order 2 removal routines
 428                         add_vertex(0,0,0,3,1,2);
 429                         add_vertex(1,0,1,3,2,0);
 430                         add_vertex(1,1,0,3,0,1);
 431                         add_vertex(2,0,0,6,4,2,1,0);
 432                         add_vertex(3,1,0,3,6,8,5);
 433                         add_vertex(3,2,0,4);
 434                         add_vertex(4,0,0,4,3,7,8);
 435                         add_vertex(5,0,0,6);
 436                         add_vertex(4,1,0,6,4);
 437         }
 438
 439         construct_relations();
 440         neighbor.label_facets();
 441 }
 442
 443 /** Adds an order one vertex to the memory structure, and specifies its edge.
 444  * \param[in] (x,y,z) are the coordinates of the vertex
 445  * \param[in] a is the first and only edge of this vertex
 446  */
 447 template<class n_option>
 448 void voronoicell_base<n_option>::add_vertex(fpoint x,fpoint y,fpoint z,int a) {
 449         pts[3*p]=x;pts[3*p+1]=y;pts[3*p+2]=z;nu[p]=1;
 450         if (mem[1]==mec[1]) add_memory(1);
 451         neighbor.set_pointer(p,1);
 452         int *q=mep[1]+3*mec[1]++;ed[p]=q;
 453         q[0]=a;q[2]=p++;
 454 }
 455
 456 /** Adds an order 2 vertex to the memory structure, and specifies its edges. */
 457 template<class n_option>
 458 void voronoicell_base<n_option>::add_vertex(fpoint x,fpoint y,fpoint z,int a,int b) {
 459         pts[3*p]=x;pts[3*p+1]=y;pts[3*p+2]=z;nu[p]=2;
 460         if (mem[2]==mec[2]) add_memory(2);
 461         neighbor.set_pointer(p,2);
 462         int *q=mep[2]+5*mec[2]++;ed[p]=q;
 463         q[0]=a;q[1]=b;q[4]=p++;
 464 }
 465
 466 /** Adds an order 3 vertex to the memory structure, and specifies its edges. */
 467 template<class n_option>
 468 void voronoicell_base<n_option>::add_vertex(fpoint x,fpoint y,fpoint z,int a,int b,int c) {
 469         pts[3*p]=x;pts[3*p+1]=y;pts[3*p+2]=z;nu[p]=3;
 470         if (mem[3]==mec[3]) add_memory(3);
 471         neighbor.set_pointer(p,3);
 472         int *q=mep[3]+7*mec[3]++;ed[p]=q;
 473         q[0]=a;q[1]=b;q[2]=c;q[6]=p++;
 474 }
 475
 476 /** Adds an order 4 vertex to the memory structure, and specifies its edges. */
 477 template<class n_option>
 478 void voronoicell_base<n_option>::add_vertex(fpoint x,fpoint y,fpoint z,int a,int b,int c,int d) {
 479         pts[3*p]=x;pts[3*p+1]=y;pts[3*p+2]=z;nu[p]=4;
 480         if (mem[4]==mec[4]) add_memory(4);
 481         neighbor.set_pointer(p,4);
 482         int *q=mep[4]+9*mec[4]++;ed[p]=q;
 483         q[0]=a;q[1]=b;q[2]=c;q[3]=d;q[8]=p++;
 484 }
 485
 486 /** Adds an order 5 vertex to the memory structure, and specifies its edges. */
 487 template<class n_option>
 488 void voronoicell_base<n_option>::add_vertex(fpoint x,fpoint y,fpoint z,int a,int b,int c,int d,int e) {
 489         pts[3*p]=x;pts[3*p+1]=y;pts[3*p+2]=z;nu[p]=5;
 490         if (mem[5]==mec[5]) add_memory(5);
 491         neighbor.set_pointer(p,5);
 492         int *q=mep[5]+11*mec[5]++;ed[p]=q;
 493         q[0]=a;q[1]=b;q[2]=c;q[3]=d;q[4]=e;q[10]=p++;
 494 }
 495
 496 /** Checks that the relational table of the Voronoi cell is accurate, and prints
 497  * out any errors. This algorithm is O(p), so running it every time the plane
 498  * routine is called will result in a significant slowdown. */
 499 template<class n_option>
 500 inline void voronoicell_base<n_option>::check_relations() {
 501         int i,j;
 502         for(i=0;i<p;i++) {
 503                 for(j=0;j<nu[i];j++) {
 504                         if (ed[ed[i][j]][ed[i][nu[i]+j]]!=i) cout << "Relational error at point " << i << ", edge " << j << "." << endl;
 505                 }
 506         }
 507 }
 508
 509 /** This routine checks for any two vertices that are connected by more than one
 510  * edge. The plane algorithm is designed so that this should not happen, so any
 511  * occurrences are most likely errors. Note that the routine is O(p), so
 512  * running it every time the plane routine is called will result in a significant
 513  * slowdown. */
 514 template<class n_option>
 515 inline void voronoicell_base<n_option>::check_duplicates() {
 516         int i,j,k;
 517         for(i=0;i<p;i++) {
 518                 for(j=1;j<nu[i];j++) {
 519                         for(k=0;k<j;k++) {
 520                                 if (ed[i][j]==ed[i][k]) cout << "Duplicate edges: (" << i << "," << j << ") and (" << i << "," << k << ") [" << ed[i][j] << "]" << endl;
 521                         }
 522                 }
 523         }
 524 }
 525
 526 /** Constructs the relational table if the edges have been specified. */
 527 template<class n_option>
 528 inline void voronoicell_base<n_option>::construct_relations() {
 529         int i,j,k,l;
 530         for(i=0;i<p;i++) for(j=0;j<nu[i];j++) {
 531                 k=ed[i][j];
 532                 l=0;
 533                 while(ed[k][l]!=i) {
 534                         l++;
 535                         if (l==nu[k]) throw fatal_error("Relation table construction failed");
 536                 }
 537                 ed[i][nu[i]+j]=l;
 538         }
 539 }
 540
 541 /** Cuts the Voronoi cell by a particle whose center is at a separation of
 542  * (x,y,z) from the cell center. The value of rsq should be initially set to
 543  * \f$x^2+y^2+z^2\f$. */
 544 template<class n_option>
 545 bool voronoicell_base<n_option>::nplane(fpoint x,fpoint y,fpoint z,fpoint rsq,int p_id) {
 546         int count=0,i,j,k,lp=up,cp,qp,rp,stack=0;stack2=0;
 547         int us=0,ls=0,qs,iqs,cs,uw,qw,lw;
 548         int *edp,*edd;
 549         fpoint u,l,r,q;bool complicated_setup=false,new_double_edge=false,double_edge=false;
 550         //Initialize the safe testing routine
 551         sure.init(x,y,z,rsq);
 552
 553         //Test approximately sqrt(n)/4 points for their proximity to the plane
 554         //and keep the one which is closest
 555         uw=sure.test(up,u);
 556
 557         // Starting from an initial guess, we now move from vertex to vertex,
 558         // to try and find an edge which intersects the cutting plane,
 559         // or a vertex which is on the plane
 560         try {
 561                 if(uw==1) {
 562                         // The test point is inside the cutting plane.
 563                         us=0;
 564                         do {
 565                                 lp=ed[up][us];
 566                                 lw=sure.test(lp,l);
 567                                 if(l<u) break;
 568                                 us++;
 569                         } while (us<nu[up]);
 570
 571                         if(us==nu[up]) {
 572                                 return false;
 573                         }
 574
 575                         ls=ed[up][nu[up]+us];
 576                         while(lw==1) {
 577                                 if(++count>=p) throw true;
 578                                 u=l;up=lp;
 579                                 for(us=0;us<ls;us++) {
 580                                         lp=ed[up][us];
 581                                         lw=sure.test(lp,l);
 582                                         if(l<u) break;
 583                                 }
 584                                 if(us==ls) {
 585                                         us++;
 586                                         while(us<nu[up]) {
 587                                                 lp=ed[up][us];
 588                                                 lw=sure.test(lp,l);
 589                                                 if(l<u) break;
 590                                                 us++;
 591                                         }
 592                                         if(us==nu[up]) {
 593                                                 return false;
 594                                         }
 595                                 }
 596                                 ls=ed[up][nu[up]+us];
 597                         }
 598
 599                         // If the last point in the iteration is within the
 600                         // plane, we need to do the complicated setup
 601                         // routine. Otherwise, we use the regular iteration.
 602                         if (lw==0) {
 603                                 up=lp;
 604                                 complicated_setup=true;
 605                         } else complicated_setup=false;
 606                 } else if (uw==-1) {
 607                         us=0;
 608                         do {
 609                                 qp=ed[up][us];
 610                                 qw=sure.test(qp,q);
 611                                 if(u<q) break;
 612                                 us++;
 613                         } while (us<nu[up]);
 614                         if(us==nu[up]) return true;
 615
 616                         while(qw==-1) {
 617                                 qs=ed[up][nu[up]+us];
 618                                 if(++count>=p) throw true;
 619                                 u=q;up=qp;
 620                                 for(us=0;us<qs;us++) {
 621                                         qp=ed[up][us];
 622                                         qw=sure.test(qp,q);
 623                                         if(u<q) break;
 624                                 }
 625                                 if(us==qs) {
 626                                         us++;
 627                                         while(us<nu[up]) {
 628                                                 qp=ed[up][us];
 629                                                 qw=sure.test(qp,q);
 630                                                 if(u<q) break;
 631                                                 us++;
 632                                         }
 633                                         if(us==nu[up]) return true;
 634                                 }
 635                         }
 636                         if (qw==1) {
 637                                 lp=up;ls=us;l=u;
 638                                 up=qp;us=ed[lp][nu[lp]+ls];u=q;
 639                                 complicated_setup=false;
 640                         } else {
 641                                 up=qp;
 642                                 complicated_setup=true;
 643                         }
 644                 } else {
 645
 646                         // Our original test point was on the plane, so we
 647                         // automatically head for the complicated setup
 648                         // routine
 649                         complicated_setup=true;
 650                 }
 651         }
 652         catch(bool except) {
 653                 // This routine is a fall-back, in case floating point errors
 654                 // cause the usual search routine to fail. In the fall-back
 655                 // routine, we just test every edge to find one straddling
 656                 // the plane.
 657 #if VOROPP_VERBOSE >=1
 658                 cerr << "Bailed out of convex calculation\n";
 659 #endif
 660                 qw=1;lw=0;
 661                 for(qp=0;qp<p;qp++) {
 662                         qw=sure.test(qp,q);
 663                         if (qw==1) {
 664
 665                                 // The point is inside the cutting space. Now
 666                                 // see if we can find a neighbor which isn't.
 667                                 for(us=0;us<nu[qp];us++) {
 668                                         lp=ed[qp][us];
 669                                         if(lp<qp) {
 670                                                 lw=sure.test(lp,l);
 671                                                 if (lw!=1) break;
 672                                         }
 673                                 }
 674                                 if(us<nu[qp]) {
 675                                         up=qp;
 676                                         if(lw==0) {
 677                                                 complicated_setup=true;
 678                                         } else {
 679                                                 complicated_setup=false;
 680                                                 u=q;
 681                                                 ls=ed[up][nu[up]+us];
 682                                         }
 683                                         break;
 684                                 }
 685                         } else if (qw==-1) {
 686
 687                                 // The point is outside the cutting space. See
 688                                 // if we can find a neighbor which isn't.
 689                                 for(ls=0;ls<nu[qp];ls++) {
 690                                         up=ed[qp][ls];
 691                                         if(up<qp) {
 692                                                 uw=sure.test(up,u);
 693                                                 if (uw!=-1) break;
 694                                         }
 695                                 }
 696                                 if(ls<nu[qp]) {
 697                                         if(uw==0) {
 698                                                 up=qp;
 699                                                 complicated_setup=true;
 700                                         } else {
 701                                                 complicated_setup=false;
 702                                                 lp=qp;l=q;
 703                                                 us=ed[lp][nu[lp]+ls];
 704                                         }
 705                                         break;
 706                                 }
 707                         } else {
 708
 709                                 // The point is in the plane, so we just
 710                                 // proceed with the complicated setup routine
 711                                 up=qp;
 712                                 complicated_setup=true;
 713                                 break;
 714                         }
 715                 }
 716                 if(qp==p) return qw==-1?true:false;
 717         }
 718
 719         // We're about to add the first point of the new facet. In either
 720         // routine, we have to add a point, so first check there's space for
 721         // it.
 722         if(p==current_vertices) add_memory_vertices();
 723
 724         if (complicated_setup) {
 725                 // The search algorithm found a point which is on the cutting
 726                 // plane. We leave that point in place, and create a new one at
 727                 // the same location.
 728                 pts[3*p]=pts[3*up];
 729                 pts[3*p+1]=pts[3*up+1];
 730                 pts[3*p+2]=pts[3*up+2];
 731
 732                 // Search for a collection of edges of the test vertex which
 733                 // are outside of the cutting space. Begin by testing the
 734                 // zeroth edge.
 735                 i=0;
 736                 lp=ed[up][0];
 737                 lw=sure.test(lp,l);
 738                 if(lw!=-1) {
 739
 740                         // The first edge is either inside the cutting space,
 741                         // or lies within the cutting plane. Test the edges
 742                         // sequentially until we find one that is outside.
 743                         rp=lw;
 744                         do {
 745                                 i++;
 746
 747                                 // If we reached the last edge with no luck
 748                                 // then all of the vertices are inside
 749                                 // or on the plane, so the cell is completely
 750                                 // deleted
 751                                 if (i==nu[up]) return false;
 752                                 lp=ed[up][i];
 753                                 lw=sure.test(lp,l);
 754                         } while (lw!=-1);
 755                         j=i+1;
 756
 757                         // We found an edge outside the cutting space. Keep
 758                         // moving through these edges until we find one that's
 759                         // inside or on the plane.
 760                         while(j<nu[up]) {
 761                                 lp=ed[up][j];
 762                                 lw=sure.test(lp,l);
 763                                 if (lw!=-1) break;
 764                                 j++;
 765                         }
 766
 767                         // Compute the number of edges for the new vertex. In
 768                         // general it will be the number of outside edges
 769                         // found, plus two. But we need to recognize the
 770                         // special case when all but one edge is outside, and
 771                         // the remaining one is on the plane. For that case we
 772                         // have to reduce the edge count by one to prevent
 773                         // doubling up.
 774                         if(j==nu[up]&&i==1&&rp==0) {
 775                                 nu[p]=nu[up];
 776                                 double_edge=true;
 777                         } else nu[p]=j-i+2;
 778                         k=1;
 779
 780                         // Add memory for the new vertex if needed, and
 781                         // initialize
 782                         while (nu[p]>=current_vertex_order) add_memory_vorder();
 783                         if (mec[nu[p]]==mem[nu[p]]) add_memory(nu[p]);
 784                         neighbor.set_pointer(p,nu[p]);
 785                         ed[p]=mep[nu[p]]+(2*nu[p]+1)*mec[nu[p]]++;
 786                         ed[p][2*nu[p]]=p;
 787
 788                         // Copy the edges of the original vertex into the new
 789                         // one. Delete the edges of the original vertex, and
 790                         // update the relational table.
 791                         us=cycle_down(i,up);
 792                         while(i<j) {
 793                                 qp=ed[up][i];
 794                                 qs=ed[up][nu[up]+i];
 795                                 neighbor.copy(p,k,up,i);
 796                                 ed[p][k]=qp;
 797                                 ed[p][nu[p]+k]=qs;
 798                                 ed[qp][qs]=p;
 799                                 ed[qp][nu[qp]+qs]=k;
 800                                 ed[up][i]=-1;
 801                                 i++;k++;
 802                         }
 803                         qs=i==nu[up]?0:i;
 804                 } else {
 805
 806                         // In this case, the zeroth edge is outside the cutting
 807                         // plane. Begin by searching backwards from the last
 808                         // edge until we find an edge which isn't outside.
 809                         i=nu[up]-1;
 810                         lp=ed[up][i];
 811                         lw=sure.test(lp,l);
 812                         while(lw==-1) {
 813                                 i--;
 814
 815                                 // If i reaches zero, then we have a point in
 816                                 // the plane all of whose edges are outside
 817                                 // the cutting space, so we just exit
 818                                 if (i==0) return true;
 819                                 lp=ed[up][i];
 820                                 lw=sure.test(lp,l);
 821                         }
 822
 823                         // Now search forwards from zero
 824                         j=1;
 825                         qp=ed[up][j];
 826                         qw=sure.test(qp,q);
 827                         while(qw==-1) {
 828                                 j++;
 829                                 qp=ed[up][j];
 830                                 qw=sure.test(qp,l);
 831                         }
 832
 833                         // Compute the number of edges for the new vertex. In
 834                         // general it will be the number of outside edges
 835                         // found, plus two. But we need to recognize the
 836                         // special case when all but one edge is outside, and
 837                         // the remaining one is on the plane. For that case we
 838                         // have to reduce the edge count by one to prevent
 839                         // doubling up.
 840                         if (i==j&&qw==0) {
 841                                 double_edge=true;
 842                                 nu[p]=nu[up];
 843                         } else {
 844                                 nu[p]=nu[up]-i+j+1;
 845                         }
 846
 847                         // Add memory to store the vertex if it doesn't exist
 848                         // already
 849                         k=1;
 850                         while(nu[p]>=current_vertex_order) add_memory_vorder();
 851                         if (mec[nu[p]]==mem[nu[p]]) add_memory(nu[p]);
 852
 853                         // Copy the edges of the original vertex into the new
 854                         // one. Delete the edges of the original vertex, and
 855                         // update the relational table.
 856                         neighbor.set_pointer(p,nu[p]);
 857                         ed[p]=mep[nu[p]]+(2*nu[p]+1)*mec[nu[p]]++;
 858                         ed[p][2*nu[p]]=p;
 859                         us=i++;
 860                         while(i<nu[up]) {
 861                                 qp=ed[up][i];
 862                                 qs=ed[up][nu[up]+i];
 863                                 neighbor.copy(p,k,up,i);
 864                                 ed[p][k]=qp;
 865                                 ed[p][nu[p]+k]=qs;
 866                                 ed[qp][qs]=p;
 867                                 ed[qp][nu[qp]+qs]=k;
 868                                 ed[up][i]=-1;
 869                                 i++;k++;
 870                         }
 871                         i=0;
 872                         while(i<j) {
 873                                 qp=ed[up][i];
 874                                 qs=ed[up][nu[up]+i];
 875                                 neighbor.copy(p,k,up,i);
 876                                 ed[p][k]=qp;
 877                                 ed[p][nu[p]+k]=qs;
 878                                 ed[qp][qs]=p;
 879                                 ed[qp][nu[qp]+qs]=k;
 880                                 ed[up][i]=-1;
 881                                 i++;k++;
 882                         }
 883                         qs=j;
 884                 }
 885                 if (!double_edge) {
 886                         neighbor.copy(p,k,up,qs);
 887                         neighbor.set(p,0,p_id);
 888                 } else neighbor.copy(p,0,up,qs);
 889
 890                 // Add this point to the auxiliary delete stack
 891                 if (stack2==current_delete2_size) add_memory_ds2();
 892                 ds2[stack2++]=up;
 893
 894                 // Look at the edges on either side of the group that was
 895                 // detected. We're going to commence facet computation by
 896                 // moving along one of them. We are going to end up coming back
 897                 // along the other one.
 898                 cs=k;
 899                 qp=up;q=u;
 900                 i=ed[up][us];
 901                 us=ed[up][nu[up]+us];
 902                 up=i;
 903                 ed[qp][2*nu[qp]]=-p;
 904
 905         } else {
 906
 907                 // The search algorithm found an intersected edge between the
 908                 // points lp and up. Create a new vertex between them which
 909                 // lies on the cutting plane. Since u and l differ by at least
 910                 // the tolerance, this division should never screw up.
 911                 if (stack==current_delete_size) add_memory_ds();
 912                 ds[stack++]=up;
 913                 r=1/(u-l);
 914                 pts[3*p]=(pts[3*lp]*u-pts[3*up]*l)*r;
 915                 pts[3*p+1]=(pts[3*lp+1]*u-pts[3*up+1]*l)*r;
 916                 pts[3*p+2]=(pts[3*lp+2]*u-pts[3*up+2]*l)*r;
 917
 918                 // This point will always have three edges. Connect one of them
 919                 // to lp.
 920                 nu[p]=3;
 921                 if (mec[3]==mem[3]) add_memory(3);
 922                 neighbor.set_pointer(p,3);
 923                 neighbor.set(p,0,p_id);
 924                 neighbor.copy(p,1,up,us);
 925                 neighbor.copy(p,2,lp,ls);
 926                 ed[p]=mep[3]+7*mec[3]++;
 927                 ed[p][6]=p;
 928                 ed[up][us]=-1;
 929                 ed[lp][ls]=p;
 930                 ed[lp][nu[lp]+ls]=1;
 931                 ed[p][1]=lp;
 932                 ed[p][nu[p]+1]=ls;
 933                 cs=2;
 934
 935                 // Set the direction to move in
 936                 qs=cycle_up(us,up);
 937                 qp=up;q=u;
 938         }
 939
 940         // When the code reaches here, we have initialized the first point, and
 941         // we have a direction for moving it to construct the rest of the facet
 942         cp=p;rp=p;p++;
 943         while(qp!=up||qs!=us) {
 944
 945                 // We're currently tracing round an intersected facet. Keep
 946                 // moving around it until we find a point or edge which
 947                 // intersects the plane.
 948                 lp=ed[qp][qs];
 949                 lw=sure.test(lp,l);
 950
 951                 if (lw==1) {
 952
 953                         // The point is still in the cutting space. Just add it
 954                         // to the delete stack and keep moving.
 955                         if (stack==current_delete_size) add_memory_ds();
 956                         qs=cycle_up(ed[qp][nu[qp]+qs],lp);
 957                         qp=lp;
 958                         q=l;
 959                         ds[stack++]=qp;
 960
 961                 } else if (lw==-1) {
 962
 963                         // The point is outside of the cutting space, so we've
 964                         // found an intersected edge. Introduce a regular point
 965                         // at the point of intersection. Connect it to the
 966                         // point we just tested. Also connect it to the previous
 967                         // new point in the facet we're constructing.
 968                         if(p==current_vertices) add_memory_vertices();
 969                         r=1/(q-l);
 970                         pts[3*p]=(pts[3*lp]*q-pts[3*qp]*l)*r;
 971                         pts[3*p+1]=(pts[3*lp+1]*q-pts[3*qp+1]*l)*r;
 972                         pts[3*p+2]=(pts[3*lp+2]*q-pts[3*qp+2]*l)*r;
 973                         nu[p]=3;
 974                         if (mec[3]==mem[3]) add_memory(3);
 975                         ls=ed[qp][qs+nu[qp]];
 976                         neighbor.set_pointer(p,3);
 977                         neighbor.set(p,0,p_id);
 978                         neighbor.copy(p,1,qp,qs);
 979                         neighbor.copy(p,2,lp,ls);
 980                         ed[p]=mep[3]+7*mec[3]++;
 981                         ed[p][6]=p;
 982                         ed[lp][ls]=p;
 983                         ed[lp][nu[lp]+ls]=1;
 984                         ed[p][1]=lp;
 985                         ed[p][0]=cp;
 986                         ed[p][nu[p]+1]=ls;
 987                         ed[p][nu[p]]=cs;
 988                         ed[cp][cs]=p;
 989                         ed[cp][nu[cp]+cs]=0;
 990                         ed[qp][qs]=-1;
 991                         qs=cycle_up(qs,qp);
 992                         cp=p++;
 993                         cs=2;
 994                 } else {
 995
 996                         // We've found a point which is on the cutting plane.
 997                         // We're going to introduce a new point right here, but
 998                         // first we need to figure out the number of edges it
 999                         // has.
1000                         if(p==current_vertices) add_memory_vertices();
1001
1002                         // If the previous vertex detected a double edge, our
1003                         // new vertex will have one less edge.
1004                         k=double_edge?0:1;
1005                         qs=ed[qp][nu[qp]+qs];
1006                         qp=lp;
1007                         iqs=qs;
1008
1009                         // Start testing the edges of the current point until
1010                         // we find one which isn't outside the cutting space
1011                         do {
1012                                 k++;
1013                                 qs=cycle_up(qs,qp);
1014                                 lp=ed[qp][qs];
1015                                 lw=sure.test(lp,l);
1016                         } while (lw==-1);
1017
1018                         // Now we need to find out whether this marginal vertex
1019                         // we are on has been visited before, because if that's
1020                         // the case, we need to add vertices to the existing
1021                         // new vertex, rather than creating a fresh one. We also
1022                         // need to figure out whether we're in a case where we
1023                         // might be creating a duplicate edge.
1024                         j=-ed[qp][2*nu[qp]];
1025                         if(qp==up&&qs==us) {
1026
1027                                 // If we're heading into the final part of the
1028                                 // new facet, then we never worry about the
1029                                 // duplicate edge calculation.
1030                                 new_double_edge=false;
1031                                 if(j>0) k+=nu[j];
1032                         } else {
1033                                 if (j>0) {
1034
1035                                         // This vertex was visited before, so
1036                                         // count those vertices to the ones we
1037                                         // already have.
1038                                         k+=nu[j];
1039
1040                                         // The only time when we might make a
1041                                         // duplicate edge is if the point we're
1042                                         // going to move to next is also a
1043                                         // marginal point, so test for that
1044                                         // first.
1045                                         if(lw==0) {
1046
1047                                                 // Now see whether this marginal point
1048                                                 // has been visited before.
1049                                                 i=-ed[lp][2*nu[lp]];
1050                                                 if(i>0) {
1051
1052                                                         // Now see if the last edge of that other
1053                                                         // marginal point actually ends up here.
1054                                                         if(ed[i][nu[i]-1]==j) {
1055                                                                 new_double_edge=true;
1056                                                                 k-=1;
1057                                                         } else new_double_edge=false;
1058                                                 } else {
1059
1060                                                         // That marginal point hasn't been visited
1061                                                         // before, so we probably don't have to worry
1062                                                         // about duplicate edges, except in the
1063                                                         // case when that's the way into the end
1064                                                         // of the facet, because that way always creates
1065                                                         // an edge.
1066                                                         if (j==rp&&lp==up&&ed[qp][nu[qp]+qs]==us) {
1067                                                                 new_double_edge=true;
1068                                                                 k-=1;
1069                                                         } else new_double_edge=false;
1070                                                 }
1071                                         } else new_double_edge=false;
1072                                 } else {
1073
1074                                         // The vertex hasn't been visited
1075                                         // before, but let's see if it's
1076                                         // marginal
1077                                         if(lw==0) {
1078
1079                                                 // If it is, we need to check
1080                                                 // for the case that it's a
1081                                                 // small branch, and that we're
1082                                                 // heading right back to where
1083                                                 // we came from
1084                                                 i=-ed[lp][2*nu[lp]];
1085                                                 if(i==cp) {
1086                                                         new_double_edge=true;
1087                                                         k-=1;
1088                                                 } else new_double_edge=false;
1089                                         } else new_double_edge=false;
1090                                 }
1091                         }
1092
1093                         // k now holds the number of edges of the new vertex
1094                         // we are forming. Add memory for it if it doesn't exist
1095                         // already.
1096                         while(k>=current_vertex_order) add_memory_vorder();
1097                         if (mec[k]==mem[k]) add_memory(k);
1098
1099                         // Now create a new vertex with order k, or augment
1100                         // the existing one
1101                         if(j>0) {
1102
1103                                 // If we're augmenting a vertex but we don't
1104                                 // actually need any more edges, just skip this
1105                                 // routine to avoid memory confusion
1106                                 if(nu[j]!=k) {
1107                                         // Allocate memory and copy the edges
1108                                         // of the previous instance into it
1109                                         neighbor.set_aux1(k);
1110                                         edp=mep[k]+(2*k+1)*mec[k]++;
1111                                         i=0;
1112                                         while(i<nu[j]) {
1113                                                 neighbor.copy_aux1(j,i);
1114                                                 edp[i]=ed[j][i];
1115                                                 edp[k+i]=ed[j][nu[j]+i];
1116                                                 i++;
1117                                         }
1118                                         edp[2*k]=j;
1119
1120                                         // Remove the previous instance with
1121                                         // fewer vertices from the memory
1122                                         // structure
1123                                         edd=mep[nu[j]]+(2*nu[j]+1)*--mec[nu[j]];
1124                                         if(edd!=ed[j]) {
1125                                                 for(lw=0;lw<=2*nu[j];lw++) ed[j][lw]=edd[lw];
1126                                                 neighbor.set_aux2_copy(j,nu[j]);
1127                                                 neighbor.copy_pointer(edd[2*nu[j]],j);
1128                                                 ed[edd[2*nu[j]]]=ed[j];
1129                                         }
1130                                         neighbor.set_to_aux1(j);
1131                                         ed[j]=edp;
1132                                 } else i=nu[j];
1133                         } else {
1134
1135                                 // Allocate a new vertex of order k
1136                                 neighbor.set_pointer(p,k);
1137                                 ed[p]=mep[k]+(2*k+1)*mec[k]++;
1138                                 ed[p][2*k]=p;
1139                                 if (stack2==current_delete2_size) add_memory_ds2();
1140                                 ds2[stack2++]=qp;
1141                                 pts[3*p]=pts[3*qp];
1142                                 pts[3*p+1]=pts[3*qp+1];
1143                                 pts[3*p+2]=pts[3*qp+2];
1144                                 ed[qp][2*nu[qp]]=-p;
1145                                 j=p++;
1146                                 i=0;
1147                         }
1148                         nu[j]=k;
1149
1150                         // Unless the previous case was a double edge, connect
1151                         // the first available edge of the new vertex to the
1152                         // last one in the facet
1153                         if (!double_edge) {
1154                                 ed[j][i]=cp;
1155                                 ed[j][nu[j]+i]=cs;
1156                                 neighbor.set(j,i,p_id);
1157                                 ed[cp][cs]=j;
1158                                 ed[cp][nu[cp]+cs]=i;
1159                                 i++;
1160                         }
1161
1162                         // Copy in the edges of the underlying vertex,
1163                         // and do one less if this was a double edge
1164                         qs=iqs;
1165                         while(i<(new_double_edge?k:k-1)) {
1166                                 qs=cycle_up(qs,qp);
1167                                 lp=ed[qp][qs];ls=ed[qp][nu[qp]+qs];
1168                                 neighbor.copy(j,i,qp,qs);
1169                                 ed[j][i]=lp;
1170                                 ed[j][nu[j]+i]=ls;
1171                                 ed[lp][ls]=j;
1172                                 ed[lp][nu[lp]+ls]=i;
1173                                 ed[qp][qs]=-1;
1174                                 i++;
1175                         }
1176                         qs=cycle_up(qs,qp);
1177                         cs=i;
1178                         cp=j;
1179                         neighbor.copy(j,new_double_edge?0:cs,qp,qs);
1180
1181                         // Update the double_edge flag, to pass it
1182                         // to the next instance of this routine
1183                         double_edge=new_double_edge;
1184                 }
1185         }
1186
1187         // Connect the final created vertex to the initial one
1188         ed[cp][cs]=rp;
1189         ed[rp][0]=cp;
1190         ed[cp][nu[cp]+cs]=0;
1191         ed[rp][nu[rp]+0]=cs;
1192
1193         // Delete points: first, remove any duplicates
1194         i=0;
1195         while(i<stack) {
1196                 j=ds[i];
1197                 if(ed[j][nu[j]]!=-1) {
1198                         ed[j][nu[j]]=-1;
1199                         i++;
1200                 } else ds[i]=ds[--stack];
1201         }
1202
1203         // Add the points in the auxiliary delete stack,
1204         // and reset their back pointers
1205         for(i=0;i<stack2;i++) {
1206                 j=ds2[i];
1207                 ed[j][2*nu[j]]=j;
1208                 if(ed[j][nu[j]]!=-1) {
1209                         ed[j][nu[j]]=-1;
1210                         if (stack==current_delete_size) add_memory_ds();
1211                         ds[stack++]=j;
1212                 }
1213         }
1214
1215         // Scan connections and add in extras
1216         for(i=0;i<stack;i++) {
1217                 cp=ds[i];
1218                 for(j=0;j<nu[cp];j++) {
1219                         qp=ed[cp][j];
1220                         if(qp!=-1) {
1221                                 if (ed[qp][nu[qp]]!=-1) {
1222                                         if (stack==current_delete_size) add_memory_ds();
1223                                         ds[stack++]=qp;
1224                                         ed[qp][nu[qp]]=-1;
1225                                 }
1226                         }
1227                 }
1228         }
1229         up=0;
1230
1231         // Delete them from the array structure
1232         while(stack>0) {
1233                 while(ed[--p][nu[p]]==-1) {
1234                         j=nu[p];
1235                         mec[j]--;
1236                         for(i=0;i<=2*j;i++) ed[p][i]=(mep[j]+(2*j+1)*mec[j])[i];
1237                         neighbor.set_aux2_copy(p,j);
1238                         neighbor.copy_pointer(ed[p][2*j],p);
1239                         ed[ed[p][2*j]]=ed[p];
1240                 }
1241                 up=ds[--stack];
1242                 if (up<p) {
1243
1244                         // Vertex management
1245                         pts[3*up]=pts[3*p];
1246                         pts[3*up+1]=pts[3*p+1];
1247                         pts[3*up+2]=pts[3*p+2];
1248
1249                         // Memory management
1250                         j=nu[up];
1251                         mec[j]--;
1252                         for(i=0;i<=2*j;i++) ed[up][i]=(mep[j]+(2*j+1)*mec[j])[i];
1253                         neighbor.set_aux2_copy(up,j);
1254                         neighbor.copy_pointer(ed[up][2*j],up);
1255                         neighbor.copy_pointer(up,p);
1256                         ed[ed[up][2*j]]=ed[up];
1257
1258                         // Edge management
1259                         ed[up]=ed[p];
1260                         nu[up]=nu[p];
1261                         for(i=0;i<nu[up];i++) {
1262                                 if (ed[up][i]==-1) throw fatal_error("fishy");
1263                                 ed[ed[up][i]][ed[up][nu[up]+i]]=up;
1264                         }
1265                         ed[up][2*nu[up]]=up;
1266                 } else up=p++;
1267         }
1268
1269         // Check for any vertices of zero order
1270         if (mec[0]>0) throw fatal_error("Zero order vertex formed");
1271
1272         // Collapse any order 2 vertices and exit
1273         return collapse_order2();
1274 }
1275
1276 /** During the creation of a new facet in the plane routine, it is possible
1277  * that some order two vertices may arise. This routine removes them.
1278  * Suppose an order two vertex joins c and d. If there's a edge between
1279  * c and d already, then the order two vertex is just removed; otherwise,
1280  * the order two vertex is removed and c and d are joined together directly.
1281  * It is possible this process will create order two or order one vertices,
1282  * and the routine is continually run until all of them are removed.
1283  * \return false if the vertex removal was unsuccessful, indicative of
1284  * the cell reducing to zero volume and disappearing; true if the vertex
1285  * removal was successful. */
1286 template<class n_option>
1287 inline bool voronoicell_base<n_option>::collapse_order2() {
1288         if(!collapse_order1()) return false;
1289         int a,b,i,j,k,l;
1290         while(mec[2]>0) {
1291
1292                 // Pick a order 2 vertex and read in its edges
1293                 i=--mec[2];
1294                 j=mep[2][5*i];k=mep[2][5*i+1];
1295                 if (j==k) {
1296 #if VOROPP_VERBOSE >=1
1297                         cerr << "Order two vertex joins itself" << endl;
1298 #endif
1299                         return false;
1300                 }
1301
1302                 // Scan the edges of j to see if joins k
1303                 for(l=0;l<nu[j];l++) {
1304                         if(ed[j][l]==k) break;
1305                 }
1306
1307                 // If j doesn't already join k, join them together.
1308                 // Otherwise delete the connection to the current
1309                 // vertex from j and k.
1310                 a=mep[2][5*i+2];b=mep[2][5*i+3];i=mep[2][5*i+4];
1311                 if(l==nu[j]) {
1312                         ed[j][a]=k;
1313                         ed[k][b]=j;
1314                         ed[j][nu[j]+a]=b;
1315                         ed[k][nu[k]+b]=a;
1316                 } else {
1317                         if (!delete_connection(j,a,false)) return false;
1318                         if (!delete_connection(k,b,true)) return false;
1319                 }
1320
1321                 // Compact the memory
1322                 --p;
1323                 if(up==i) up=0;
1324                 if(p!=i) {
1325                         if (up==p) up=i;
1326                         pts[3*i]=pts[3*p];
1327                         pts[3*i+1]=pts[3*p+1];
1328                         pts[3*i+2]=pts[3*p+2];
1329                         for(k=0;k<nu[p];k++) ed[ed[p][k]][ed[p][nu[p]+k]]=i;
1330                         neighbor.copy_pointer(i,p);
1331                         ed[i]=ed[p];
1332                         nu[i]=nu[p];
1333                         ed[i][2*nu[i]]=i;
1334                 }
1335
1336                 // Collapse any order 1 vertices if they were created
1337                 if(!collapse_order1()) return false;
1338         }
1339         return true;
1340 }
1341
1342 /** Order one vertices can potentially be created during the order two collapse
1343  * routine. This routine keeps removing them until there are none left.
1344  * \return false if the vertex removal was unsuccessful, indicative of
1345  * the cell having zero volume and disappearing; true if the vertex removal
1346  * was successful. */
1347 template<class n_option>
1348 inline bool voronoicell_base<n_option>::collapse_order1() {
1349         int i,j,k;
1350         while(mec[1]>0) {
1351                 up=0;
1352 #if VOROPP_VERBOSE >=1
1353                 cerr << "Order one collapse" << endl;
1354 #endif
1355                 i=--mec[1];
1356                 j=mep[1][3*i];k=mep[1][3*i+1];
1357                 i=mep[1][3*i+2];
1358                 if(!delete_connection(j,k,false)) return false;
1359                 --p;
1360                 if(up==i) up=0;
1361                 if(p!=i) {
1362                         if (up==p) up=i;
1363                         pts[3*i]=pts[3*p];
1364                         pts[3*i+1]=pts[3*p+1];
1365                         pts[3*i+2]=pts[3*p+2];
1366                         for(k=0;k<nu[p];k++) ed[ed[p][k]][ed[p][nu[p]+k]]=i;
1367                         neighbor.copy_pointer(i,p);
1368                         ed[i]=ed[p];
1369                         nu[i]=nu[p];
1370                         ed[i][2*nu[i]]=i;
1371                 }
1372         }
1373         return true;
1374 }
1375
1376 /** This routine deletes the kth edge of vertex j and reorganizes the memory.
1377  * If the neighbor computation is enabled, we also have to supply an
1378  * handedness flag to decide whether to preserve the plane on the left
1379  * or right of the connection.
1380  * \return false if a zero order vertex was formed, indicative of the cell
1381  * disappearing; true if the vertex removal was successful. */
1382 template<class n_option>
1383 inline bool voronoicell_base<n_option>::delete_connection(int j,int k,bool hand) {
1384         int q=hand?k:cycle_up(k,j);
1385         int i=nu[j]-1,l,*edp,*edd,m;
1386         if(i<1) {
1387                 cout << "Zero order vertex formed" << endl;
1388                 return false;
1389         }
1390         if(mec[i]==mem[i]) add_memory(i);
1391         neighbor.set_aux1(i);
1392         for(l=0;l<q;l++) neighbor.copy_aux1(j,l);
1393         while(l<i) {
1394                 neighbor.copy_aux1_shift(j,l);
1395                 l++;
1396         }
1397         edp=mep[i]+(2*i+1)*mec[i]++;
1398         edp[2*i]=j;
1399         for(l=0;l<k;l++) {
1400                 edp[l]=ed[j][l];
1401                 edp[l+i]=ed[j][l+nu[j]];
1402         }
1403         while(l<i) {
1404                 m=ed[j][l+1];
1405                 edp[l]=m;
1406                 k=ed[j][l+nu[j]+1];
1407                 edp[l+i]=k;
1408                 ed[m][nu[m]+k]--;
1409                 l++;
1410         }
1411
1412         edd=mep[nu[j]]+(2*nu[j]+1)*--mec[nu[j]];
1413         for(l=0;l<=2*nu[j];l++) ed[j][l]=edd[l];
1414         neighbor.set_aux2_copy(j,nu[j]);
1415         neighbor.set_to_aux2(edd[2*nu[j]]);
1416         neighbor.set_to_aux1(j);
1417         ed[edd[2*nu[j]]]=edd;
1418         ed[j]=edp;
1419         nu[j]=i;
1420         return true;
1421 }
1422
1423 /** Cuts a Voronoi cell using the influence of a particle at (x,y,z), first
1424  * calculating the modulus squared of this vector before passing it to the
1425  * main nplane() routine. Zero is supplied as the plane ID, which will be
1426  * ignored unless neighbor tracking is enabled.
1427  * \param[in] (x,y,z) The vector to cut the cell by. */
1428 template<class n_option>
1429 inline bool voronoicell_base<n_option>::plane(fpoint x,fpoint y,fpoint z) {
1430         fpoint rsq=x*x+y*y+z*z;
1431         return nplane(x,y,z,rsq,0);
1432 }
1433
1434 /** This version of the plane routine just makes up the plane ID to be zero. It
1435  * will only be referenced if neighbor tracking is enabled.
1436  * \param[in] (x,y,z) The vector to cut the cell by.
1437  * \param[in] rsq The modulus squared of the vector. */
1438 template<class n_option>
1439 inline bool voronoicell_base<n_option>::plane(fpoint x,fpoint y,fpoint z,fpoint rsq) {
1440         return nplane(x,y,z,rsq,0);
1441 }
1442
1443 /** This routine calculates the modulus squared of the vector before passing it
1444  * to the main nplane() routine with full arguments.
1445  * \param[in] (x,y,z) The vector to cut the cell by.
1446  * \param[in] p_id The plane ID (for neighbor tracking only).*/
1447 template<class n_option>
1448 inline bool voronoicell_base<n_option>::nplane(fpoint x,fpoint y,fpoint z,int p_id) {
1449         fpoint rsq=x*x+y*y+z*z;
1450         return nplane(x,y,z,rsq,p_id);
1451 }
1452
1453 /** This is a simple inline function for picking out the index of the next edge
1454  * counterclockwise at the current vertex.
1455  * \param[in] a The index of an edge of the current vertex.
1456  * \param[in] p The number of the vertex.
1457  * \return 0 if a=nu[p]-1, or a+1 otherwise. */
1458 template<class n_option>
1459 inline int voronoicell_base<n_option>::cycle_up(int a,int p) {
1460         return a==nu[p]-1?0:a+1;
1461 }
1462
1463 /** This is a simple inline function for picking out the index of the next edge
1464  * clockwise from the current vertex.
1465  * \param[in] a The index of an edge of the current vertex.
1466  * \param[in] p The number of the vertex.
1467  * \return nu[p]-1 if a=0, or a-1 otherwise. */
1468 template<class n_option>
1469 inline int voronoicell_base<n_option>::cycle_down(int a,int p) {
1470         return a==0?nu[p]-1:a-1;
1471 }
1472
1473 /** Calculates the volume of the Voronoi cell, by decomposing the cell into
1474  * tetrahedra extending outward from the zeroth vertex, which are evaluated
1475  * using a scalar triple product.
1476  * \return A floating point number holding the calculated volume. */
1477 template<class n_option>
1478 fpoint voronoicell_base<n_option>::volume() {
1479         const fpoint fe=1/48.0;
1480         fpoint vol=0;
1481         int i,j,k,l,m,n;
1482         fpoint ux,uy,uz,vx,vy,vz,wx,wy,wz;
1483         for(i=1;i<p;i++) {
1484                 ux=pts[0]-pts[3*i];
1485                 uy=pts[1]-pts[3*i+1];
1486                 uz=pts[2]-pts[3*i+2];
1487                 for(j=0;j<nu[i];j++) {
1488                         k=ed[i][j];
1489                         if (k>=0) {
1490                                 ed[i][j]=-1-k;
1491                                 l=cycle_up(ed[i][nu[i]+j],k);
1492                                 vx=pts[3*k]-pts[0];
1493                                 vy=pts[3*k+1]-pts[1];
1494                                 vz=pts[3*k+2]-pts[2];
1495                                 m=ed[k][l];ed[k][l]=-1-m;
1496                                 while(m!=i) {
1497                                         n=cycle_up(ed[k][nu[k]+l],m);
1498                                         wx=pts[3*m]-pts[0];
1499                                         wy=pts[3*m+1]-pts[1];
1500                                         wz=pts[3*m+2]-pts[2];
1501                                         vol+=ux*vy*wz+uy*vz*wx+uz*vx*wy-uz*vy*wx-uy*vx*wz-ux*vz*wy;
1502                                         k=m;l=n;vx=wx;vy=wy;vz=wz;
1503                                         m=ed[k][l];ed[k][l]=-1-m;
1504                                 }
1505                         }
1506                 }
1507         }
1508         for(i=0;i<p;i++) {
1509                 for(j=0;j<nu[i];j++) {
1510                         if(ed[i][j]>=0) throw fatal_error("Volume routine didn't look everywhere");
1511                         ed[i][j]=-1-ed[i][j];
1512                 }
1513         }
1514
1515         return vol*fe;
1516 }
1517
1518 /** Computes the maximum radius squared of a vertex from the center
1519  * of the cell. It can be used to determine when enough particles have
1520  * been testing an all planes that could cut the cell have been considered.
1521  * \return The maximum radius squared of a vertex.*/
1522 template<class n_option>
1523 fpoint voronoicell_base<n_option>::maxradsq() {
1524         int i;fpoint r,s;
1525         r=pts[0]*pts[0]+pts[1]*pts[1]+pts[2]*pts[2];
1526         for(i=3;i<3*p;i+=3) {
1527                 s=pts[i]*pts[i]+pts[i+1]*pts[i+1]+pts[i+2]*pts[i+2];
1528                 if(s>r) r=s;
1529         }
1530         return r;
1531 }
1532
1533 /** Outputs the edges of the Voronoi cell (in POV-Ray format) to an open file
1534  * stream, displacing the cell by an amount (x,y,z).
1535  * \param[in] &os A output stream to write to.
1536  * \param[in] (x,y,z) A displacement vector to be added to the cell's position. */
1537 template<class n_option>
1538 void voronoicell_base<n_option>::draw_pov(ostream &os,fpoint x,fpoint y,fpoint z) {
1539         int i,j,k;fpoint ux,uy,uz;
1540 //os << "merge{";
1541         for(i=0;i<p;i++) {
1542                 ux=x+0.5*pts[3*i];uy=y+0.5*pts[3*i+1];uz=z+0.5*pts[3*i+2];
1543                 if (neighbor.internal_wall(i)) os << "sphere{<" << ux << "," << uy << "," << uz << ">,r}\n";
1544                 for(j=0;j<nu[i];j++) {
1545                         k=ed[i][j];
1546                         if (k<i&&!neighbor.internal_wall(i,j)) os << "cylinder{<" << ux << "," << uy << "," << uz << ">,<" << x+0.5*pts[3*k] << "," << y+0.5*pts[3*k+1] << "," << z+0.5*pts[3*k+2] << ">,r}\n";
1547                 }
1548         }
1549 //      os << "}\n";
1550 }
1551
1552 inline bool neighbor_track::internal_wall(int i) {
1553         for(int j=0;j<vc->nu[i];j++) {
1554                 if(ne[i][j]!=-99) return true;
1555         }
1556         return false;
1557 }
1558
1559 inline bool neighbor_track::internal_wall(int i,int j) {
1560         return ne[i][j]==-99&&ne[i][j==vc->nu[i]-1?0:j+1]==-99;
1561 }
1562
1563 /** An overloaded version of the draw_pov routine, that outputs the edges of
1564  * the Voronoi cell (in POV-Ray format) to a file.
1565  * \param[in] filename The file to write to.
1566  * \param[in] (x,y,z) A displacement vector to be added to the cell's position.
1567  */
1568 template<class n_option>
1569 inline void voronoicell_base<n_option>::draw_pov(const char *filename,fpoint x,fpoint y,fpoint z) {
1570         ofstream os;
1571         os.open(filename,ofstream::out|ofstream::trunc);
1572         draw_pov(os,x,y,z);
1573         os.close();
1574 }
1575
1576 /** An overloaded version of the draw_pov routine, that outputs the edges of the
1577  * Voronoi cell (in POV-Ray format) to standard output.
1578  * \param[in] (x,y,z) A displacement vector to be added to the cell's position.
1579  */
1580 template<class n_option>
1581 inline void voronoicell_base<n_option>::draw_pov(fpoint x,fpoint y,fpoint z) {
1582         draw_pov(cout,x,y,z);
1583 }
1584
1585 /** Outputs the edges of the Voronoi cell (in gnuplot format) to an output stream.
1586  * \param[in] &os A reference to an output stream to write to.
1587  * \param[in] (x,y,z) A displacement vector to be added to the cell's position.
1588  */
1589 template<class n_option>
1590 void voronoicell_base<n_option>::draw_gnuplot(ostream &os,fpoint x,fpoint y,fpoint z) {
1591         int i,j,k;fpoint ux,uy,uz;
1592         for(i=0;i<p;i++) {
1593                 ux=x+0.5*pts[3*i];uy=y+0.5*pts[3*i+1];uz=z+0.5*pts[3*i+2];
1594                 for(j=0;j<nu[i];j++) {
1595                         k=ed[i][j];
1596                         if (ed[i][j]<i) os << ux << " " << uy << " " << uz << "\n" << x+0.5*pts[3*k] << " " << y+0.5*pts[3*k+1] << " " << z+0.5*pts[3*k+2] << "\n\n\n";
1597                 }
1598         }
1599 }
1600
1601 /** An overloaded version of the draw_gnuplot routine that writes directly to
1602  * a file.
1603  * \param[in] filename The name of the file to write to.
1604  * \param[in] (x,y,z) A displacement vector to be added to the cell's position.
1605  */
1606 template<class n_option>
1607 inline void voronoicell_base<n_option>::draw_gnuplot(const char *filename,fpoint x,fpoint y,fpoint z) {
1608         ofstream os;
1609         os.open(filename,ofstream::out|ofstream::trunc);
1610         draw_gnuplot(os,x,y,z);
1611         os.close();
1612 }
1613
1614 /** An overloaded version of the draw_gnuplot routine, that prints to the
1615  * standard output.
1616  * \param[in] (x,y,z) A displacement vector to be added to the cell's position.
1617  */
1618 template<class n_option>
1619 inline void voronoicell_base<n_option>::draw_gnuplot(fpoint x,fpoint y,fpoint z) {
1620         draw_gnuplot(cout,x,y,z);
1621 }
1622
1623 /** Outputs the Voronoi cell in the POV mesh2 format, described in section
1624  * 1.3.2.2 of the POV-Ray documentation. The mesh2 output consists of a list of
1625  * vertex vectors, followed by a list of triangular faces. The routine also
1626  * makes use of the optional inside_vector specification, which makes the mesh
1627  * object solid, so the the POV-Ray Constructive Solid Geometry (CSG) can be
1628  * applied.
1629  * \param[in] &os An output stream to write to.
1630  * \param[in] (x,y,z) A displacement vector to be added to the cell's position.
1631  */
1632 template<class n_option>
1633 inline void voronoicell_base<n_option>::draw_pov_mesh(ostream &os,fpoint x,fpoint y,fpoint z) {
1634         int i,j,k,l,m,n;
1635         os << "mesh2 {\nvertex_vectors {\n" << p << ",\n";
1636         for(i=0;i<p;i++) {
1637                 os << "<" << x+0.5*pts[3*i] << "," << y+0.5*pts[3*i+1] << "," << z+0.5*pts[3*i+2] << ">,\n";
1638         }
1639         os << "}\nface_indices {\n" << 2*(p-2) << ",\n";
1640         for(i=1;i<p;i++) {
1641                 for(j=0;j<nu[i];j++) {
1642                         k=ed[i][j];
1643                         if (k>=0) {
1644                                 ed[i][j]=-1-k;
1645                                 l=cycle_up(ed[i][nu[i]+j],k);
1646                                 m=ed[k][l];ed[k][l]=-1-m;
1647                                 while(m!=i) {
1648                                         n=cycle_up(ed[k][nu[k]+l],m);
1649                                         os << "<" << i << "," << k << "," << m << ">\n";
1650                                         k=m;l=n;
1651                                         m=ed[k][l];ed[k][l]=-1-m;
1652                                 }
1653                         }
1654                 }
1655         }
1656         for(i=0;i<p;i++) {
1657                 for(j=0;j<nu[i];j++) {
1658                         if(ed[i][j]>=0) throw fatal_error("Mesh routine didn't look everywhere");
1659                         ed[i][j]=-1-ed[i][j];
1660                 }
1661         }
1662         os << "}\ninside_vector <0,0,1>\n}\n";
1663 }
1664
1665 /** An overloaded version of the draw_pov_mesh routine, that writes directly to a
1666  * file.
1667  * \param[in] filename A filename to write to.
1668  * \param[in] (x,y,z) A displacement vector to be added to the cell's position.
1669  */
1670 template<class n_option>
1671 inline void voronoicell_base<n_option>::draw_pov_mesh(const char *filename,fpoint x,fpoint y,fpoint z) {
1672         ofstream os;
1673         os.open(filename,ofstream::out|ofstream::trunc);
1674         draw_pov_mesh(os,x,y,z);
1675         os.close();
1676 }
1677
1678 /** An overloaded version of the draw_pov_mesh routine, that prints to the
1679  * standard output.
1680  * \param[in] (x,y,z) A displacement vector to be added to the cell's position.
1681  */
1682 template<class n_option>
1683 inline void voronoicell_base<n_option>::draw_pov_mesh(fpoint x,fpoint y,fpoint z) {
1684         draw_pov_mesh(cout,x,y,z);
1685 }
1686
1687 /** Randomly perturbs the points in the Voronoi cell by an amount r. */
1688 template<class n_option>
1689 inline void voronoicell_base<n_option>::perturb(fpoint r) {
1690         for(int i=0;i<3*p;i++) {
1691                 pts[i]+=(2*fpoint(rand())/RAND_MAX-1)*r;
1692         }
1693 }
1694
1695 /** Initializes the suretest class and creates a buffer for marginal points. */
1696 suretest::suretest() : current_marginal(init_marginal) {
1697         sn=new int[2*current_marginal];
1698 }
1699
1700 /** Suretest destructor to free memory allocation. */
1701 suretest::~suretest() {
1702         delete [] sn;
1703 }
1704
1705 /** Sets up the suretest class with a particular test plane, and removes
1706  * any special cases from the table. */
1707 inline void suretest::init(fpoint x,fpoint y,fpoint z,fpoint rsq) {
1708         sc=0;px=x;py=y;pz=z;prsq=rsq;
1709 }
1710
1711 /** */
1712 inline int suretest::test(int n,fpoint &ans) {
1713         ans=px*p[3*n]+py*p[3*n+1]+pz*p[3*n+2]-prsq;
1714         if(ans<-tolerance2) {
1715                 return -1;
1716         } else if(ans>tolerance2) {
1717                 return 1;
1718         }
1719         return check_marginal(n,ans);
1720 }
1721
1722 int suretest::check_marginal(int n,fpoint &ans) {
1723         int i;
1724         for(i=0;i<sc;i+=2) if(sn[i]==n) return sn[i+1];
1725         if (sc==2*current_marginal) {
1726                 i=2*current_marginal;
1727                 if (i>max_marginal) throw fatal_error("Marginal case buffer allocation exceeded absolute maximum");
1728 #if VOROPP_VERBOSE >=2
1729                 cerr << "Marginal cases buffer scaled up to " << i << endl;
1730 #endif
1731                 int *psn=new int[2*i];
1732                 for(int j=0;j<2*current_marginal;j++) psn[j]=sn[j];
1733                 delete [] sn;sn=psn;
1734         }
1735         sn[sc++]=n;
1736         sn[sc++]=ans>tolerance?1:(ans<-tolerance?-1:0);
1737         return sn[sc-1];
1738 }
1739
1740 /** Prints the vertices, their edges, the relation table, and also notifies if
1741  * any glaring memory errors are visible. */
1742 template<class n_option>
1743 void voronoicell_base<n_option>::print_edges() {
1744         int j;
1745         for(int i=0;i<p;i++) {
1746                 cout << i << " " << nu[i] << "  ";
1747                 for(j=0;j<nu[i];j++) cout << " " << ed[i][j];
1748                 cout << "   ";
1749                 while(j<2*nu[i]) cout << " " << ed[i][j++];
1750                 cout << "    " << ed[i][j];
1751                 neighbor.print_edges(i);
1752                 cout << "    " << pts[3*i] << " " << pts[3*i+1] << " " << pts[3*i+2];
1753                 cout << " " << ed[i];
1754                 if (ed[i]>=mep[nu[i]]+mec[nu[i]]*(2*nu[i]+1)) cout << " Memory error";
1755                 cout << endl;
1756         }
1757 }
1758
1759 /** Prints out a list of all the facets and their vertices. If the neighbor option
1760  * is defined, it lists each cutting plane. */
1761 template<class n_option>
1762 void voronoicell_base<n_option>::facets(ostream &os) {
1763         int i,j,k,l,m;
1764         for(i=0;i<p;i++) {
1765                 for(j=0;j<nu[i];j++) {
1766                         k=ed[i][j];
1767                         if (k>=0) {
1768                                 neighbor.print(os,i,j);
1769                                 ed[i][j]=-1-k;
1770                                 l=cycle_up(ed[i][nu[i]+j],k);
1771                                 do {
1772                                         os << " ";
1773                                         neighbor.print(os,k,l);
1774                                         m=ed[k][l];
1775                                         ed[k][l]=-1-m;
1776                                         l=cycle_up(ed[k][nu[k]+l],m);
1777                                         k=m;
1778                                 } while (k!=i);
1779                                 os << endl;
1780                         }
1781                 }
1782         }
1783         for(i=0;i<p;i++) {
1784                 for(j=0;j<nu[i];j++) {
1785                         if(ed[i][j]>=0) throw fatal_error("Facet evaluation routine didn't look everywhere");
1786                         ed[i][j]=-1-ed[i][j];
1787                 }
1788         }
1789 }
1790
1791 /** Returns the number of faces of a computed Voronoi cell.
1792  * \return The number of faces. */
1793 template<class n_option>
1794 int voronoicell_base<n_option>::number_of_faces() {
1795         int i,j,k,l,m,s=0;
1796         for(i=0;i<p;i++) {
1797                 for(j=0;j<nu[i];j++) {
1798                         k=ed[i][j];
1799                         if (k>=0) {
1800                                 s++;
1801                                 ed[i][j]=-1-k;
1802                                 l=cycle_up(ed[i][nu[i]+j],k);
1803                                 do {
1804                                         m=ed[k][l];
1805                                         ed[k][l]=-1-m;
1806                                         l=cycle_up(ed[k][nu[k]+l],m);
1807                                         k=m;
1808                                 } while (k!=i);
1809                         }
1810                 }
1811         }
1812         for(i=0;i<p;i++) {
1813                 for(j=0;j<nu[i];j++) {
1814                         if(ed[i][j]>=0) throw fatal_error("Face evaluation routine didn't look everywhere");
1815                         ed[i][j]=-1-ed[i][j];
1816                 }
1817         }
1818         return s;
1819 }
1820
1821 /** This routine is a placeholder which just prints the ID of a
1822  * vertex.
1823  * \param[in] &os The output stream to write to.
1824  * \param[in] i The ID of a vertex.
1825  * \param[in] j The particular plane of interest (ignored in this routine). */
1826 void neighbor_none::print(ostream &os,int i,int j) {
1827         os << i;
1828 }
1829
1830 /** An overloaded version of facets() which output the results to the standard
1831  * output. */
1832 template<class n_option>
1833 inline void voronoicell_base<n_option>::facets() {
1834         facets(cout);
1835 }
1836
1837 /** An overloaded version of facets(), which outputs the results to a file.
1838  * \param[in] filename The name of the file to write to. */
1839 template<class n_option>
1840 inline void voronoicell_base<n_option>::facets(const char *filename) {
1841         ofstream os;
1842         os.open(filename,ofstream::out|ofstream::trunc);
1843         facets(os);
1844         os.close();
1845 }
1846
1847 /** Examines all the facets, and evaluates them by the number of vertices that
1848  * they have.
1849  * \param[in] &os An open output stream to write to. */
1850 template<class n_option>
1851 void voronoicell_base<n_option>::facet_statistics(ostream &os) {
1852         int *stat,*pstat,current_facet_size=init_facet_size,newc,maxf=0;
1853         stat=new int[current_facet_size];
1854         int i,j,k,l,m,q;
1855         for(i=0;i<current_facet_size;i++) stat[i]=0;
1856         for(i=0;i<p;i++) {
1857                 for(j=0;j<nu[i];j++) {
1858                         k=ed[i][j];
1859                         if (k>=0) {
1860                                 q=1;
1861                                 ed[i][j]=-1-k;
1862                                 l=cycle_up(ed[i][nu[i]+j],k);
1863                                 do {
1864                                         q++;
1865                                         m=ed[k][l];
1866                                         ed[k][l]=-1-m;
1867                                         l=cycle_up(ed[k][nu[k]+l],m);
1868                                         k=m;
1869                                 } while (k!=i);
1870                                 if (q>=current_facet_size) {
1871                                         newc=current_facet_size*2;
1872                                         pstat=new int[newc];
1873                                         for(k=0;k<current_facet_size;k++) pstat[k]=stat[k];
1874                                         while(k<newc) pstat[k]=0;
1875                                         delete [] stat;
1876                                         current_facet_size=newc;
1877                                         stat=pstat;
1878                                 }
1879                                 stat[q]++;
1880                                 if (q>maxf) maxf=q;
1881                         }
1882                 }
1883         }
1884         for(i=0;i<p;i++) {
1885                 for(j=0;j<nu[i];j++) {
1886                         if(ed[i][j]>=0) throw fatal_error("Facet statistics routine didn't look everywhere");
1887                         ed[i][j]=-1-ed[i][j];
1888                 }
1889         }
1890         for(i=0;i<=maxf;i++) os << i << " " << stat[i] << endl;
1891         delete [] stat;
1892 }
1893
1894 /** An overloaded version of facet_statistics() which outputs the results to
1895  * standard output. */
1896 template<class n_option>
1897 inline void voronoicell_base<n_option>::facet_statistics() {
1898         facet_statistics(cout);
1899 }
1900
1901 /** An overloaded version of facet_statistics() which outputs the results to
1902  * a file.
1903  * \param[in] filename The name of the file to write to. */
1904 template<class n_option>
1905 inline void voronoicell_base<n_option>::facet_statistics(const char *filename) {
1906         ofstream os;
1907         os.open(filename,ofstream::out|ofstream::trunc);
1908         facet_statistics(os);
1909         os.close();
1910 }
1911
1912 /** If the template is instantiated with the neighbor tracking turned on,
1913  * then this routine will label all the facets of the current cell. Otherwise
1914  * this routine does nothing. */
1915 template<class n_option>
1916 void voronoicell_base<n_option>::label_facets() {
1917         neighbor.label_facets();
1918 }
1919
1920 /** If the template is instantiated with the neighbor tracking turned on, then
1921  * this routine will print out a list of all the neighbors of a given cell.
1922  * Otherwise, this routine does nothing.
1923  * \param[in] &os An open output stream to write to. */
1924 template<class n_option>
1925 void voronoicell_base<n_option>::neighbors(ostream &os) {
1926         neighbor.neighbors(os);
1927 }
1928
1929 /** If the template is instantiated with the neighbor tracking turned on, then
1930  * this routine will check that the neighbor information is consistent, by
1931  * tracing around every facet, and ensuring that all the neighbor information
1932  * for that facet refers to the same neighbor. If the neighbor tracking isn't
1933  * turned on, this routine does nothing. */
1934 template<class n_option>
1935 void voronoicell_base<n_option>::check_facets() {
1936         neighbor.check_facets();
1937 }
1938
1939 /** This routine tests to see whether the cell intersects a plane by starting
1940  * from the guess point up. If up intersects, then it immediately returns true.
1941  * Otherwise, it calls the plane_intersects_track() routine.
1942  * \param[in] (x,y,z) The normal vector to the plane.
1943  * \param[in] rsq The distance along this vector of the plane.
1944  * \return false if the plane does not intersect the plane, true if it does.*/
1945 template<class n_option>
1946 bool voronoicell_base<n_option>::plane_intersects(fpoint x,fpoint y,fpoint z,fpoint rsq) {
1947         fpoint g=x*pts[3*up]+y*pts[3*up+1]+z*pts[3*up+2];
1948         if (g<rsq) return plane_intersects_track(x,y,z,rsq,g);
1949         return true;
1950 }
1951
1952 /** This routine tests to see if a cell intersects a plane. It first tests a
1953  * random sample of approximately sqrt(p)/4 points. If any of those are
1954  * intersect, then it immediately returns true. Otherwise, it takes the closest
1955  * point and passes that to plane_intersect_track() routine.
1956  * \param[in] (x,y,z) The normal vector to the plane.
1957  * \param[in] rsq The distance along this vector of the plane.
1958  * \return false if the plane does not intersect the plane, true if it does.*/
1959 template<class n_option>
1960 bool voronoicell_base<n_option>::plane_intersects_guess(fpoint x,fpoint y,fpoint z,fpoint rsq) {
1961         up=0;
1962         fpoint g=x*pts[3*up]+y*pts[3*up+1]+z*pts[3*up+2];
1963         if (g<rsq) {
1964                 int ca=1,cc=p>>3,mp=1;
1965                 fpoint m;
1966                 while(ca<cc) {
1967                         m=x*pts[3*mp]+y*pts[3*mp+1]+z*pts[3*mp+2];
1968                         if (m>g) {
1969                                 if (m>rsq) return true;
1970                                 g=m;up=mp;
1971                         }
1972                         ca+=mp++;
1973                 }
1974                 return plane_intersects_track(x,y,z,rsq,g);
1975         }
1976         return true;
1977 }
1978
1979 /* This routine tests to see if a cell intersects a plane, by tracing over the cell from
1980  * vertex to vertex, starting at up. It is meant to be called either by plane_intersects()
1981  * or plane_intersects_track(), when those routines cannot immediately resolve the case.
1982  * \param[in] (x,y,z) The normal vector to the plane.
1983  * \param[in] rsq The distance along this vector of the plane.
1984  * \param[in] g The distance of up from the plane.
1985  * \return false if the plane does not intersect the plane, true if it does.*/
1986 template<class n_option>
1987 inline bool voronoicell_base<n_option>::plane_intersects_track(fpoint x,fpoint y,fpoint z,fpoint rsq,fpoint g) {
1988         int count=0,ls,us,tp;
1989         fpoint t;
1990         // The test point is outside of the cutting space
1991         for(us=0;us<nu[up];us++) {
1992                 tp=ed[up][us];
1993                 t=x*pts[3*tp]+y*pts[3*tp+1]+z*pts[3*tp+2];
1994                 if(t>g) {
1995                         ls=ed[up][nu[up]+us];
1996                         up=tp;
1997                         while (t<rsq) {
1998                                 if (++count>=p) {
1999 #if VOROPP_VERBOSE >=1
2000                                         cerr << "Bailed out of convex calculation" << endl;
2001 #endif
2002                                         for(tp=0;tp<p;tp++) if (x*pts[3*tp]+y*pts[3*tp+1]+z*pts[3*tp+2]>rsq) return true;
2003                                         return false;
2004                                 }
2005
2006                                 // Test all the neighbors of the current point
2007                                 // and find the one which is closest to the
2008                                 // plane
2009                                 for(us=0;us<ls;us++) {
2010                                         tp=ed[up][us];
2011                                         g=x*pts[3*tp]+y*pts[3*tp+1]+z*pts[3*tp+2];
2012                                         if(g>t) break;
2013                                 }
2014                                 if (us==ls) {
2015                                         us++;
2016                                         while(us<nu[up]) {
2017                                                 tp=ed[up][us];
2018                                                 g=x*pts[3*tp]+y*pts[3*tp+1]+z*pts[3*tp+2];
2019                                                 if(g>t) break;
2020                                                 us++;
2021                                         }
2022                                         if (us==nu[up]) return false;
2023                                 }
2024                                 ls=ed[up][nu[up]+us];up=tp;t=g;
2025                         }
2026                         return true;
2027                 }
2028         }
2029         return false;
2030 }
2031
2032 /** This constructs the neighbor_track class, within a current
2033  * voronoicell_neighbor class. It allocates memory for neighbor storage in a
2034  * similar way to the voronoicell constructor. */
2035 neighbor_track::neighbor_track(voronoicell_base<neighbor_track> *ivc) : vc(ivc) {
2036         int i;
2037         mne=new int*[vc->current_vertex_order];
2038         ne=new int*[vc->current_vertices];
2039         for(i=0;i<3;i++) mne[i]=new int[init_n_vertices*i];
2040         mne[3]=new int[init_3_vertices*3];
2041         for(i=4;i<vc->current_vertex_order;i++) mne[i]=new int[init_n_vertices*i];
2042 }
2043
2044 /** The destructor for the neighbor_track class deallocates the arrays
2045  * for neighbor tracking. */
2046 neighbor_track::~neighbor_track() {
2047         for(int i=0;i<vc->current_vertex_order;i++) if (vc->mem[i]>0) delete [] mne[i];
2048         delete [] mne;
2049         delete [] ne;
2050 }
2051
2052 /** This allocates a single array for neighbor tracking.
2053  * \param[in] i the vertex order of the array to be extended.
2054  * \param[in] m the size of the array to be extended. */
2055 inline void neighbor_track::allocate(int i,int m) {
2056         mne[i]=new int[m*i];
2057 }
2058
2059 /** This increases the size of the ne[] array.
2060  * \param[in] i the new size of the array. */
2061 inline void neighbor_track::add_memory_vertices(int i) {
2062         int **pp;
2063         pp=new int*[i];
2064         for(int j=0;j<vc->current_vertices;j++) pp[j]=ne[j];
2065         delete [] ne;ne=pp;
2066 }
2067
2068 /** This increases the size of the maximum allowable vertex
2069  * order in the neighbor tracking. */
2070 inline void neighbor_track::add_memory_vorder(int i) {
2071         int **p2;
2072         p2=new int*[i];
2073         for(int j=0;j<vc->current_vertex_order;j++) p2[j]=mne[j];
2074         delete [] mne;mne=p2;
2075 }
2076
2077 /** This initializes the neighbor information for a rectangular box
2078  * and is called during the initialization routine for the voronoicell
2079  * class. */
2080 inline void neighbor_track::init() {
2081         int *q;
2082         q=mne[3];
2083         q[0]=-5;q[1]=-3;q[2]=-1;
2084         q[3]=-5;q[4]=-2;q[5]=-3;
2085         q[6]=-5;q[7]=-1;q[8]=-4;
2086         q[9]=-5;q[10]=-4;q[11]=-2;
2087         q[12]=-6;q[13]=-1;q[14]=-3;
2088         q[15]=-6;q[16]=-3;q[17]=-2;
2089         q[18]=-6;q[19]=-4;q[20]=-1;
2090         q[21]=-6;q[22]=-2;q[23]=-4;
2091         ne[0]=q;ne[1]=q+3;ne[2]=q+6;ne[3]=q+9;
2092         ne[4]=q+12;ne[5]=q+15;ne[6]=q+18;ne[7]=q+21;
2093 }
2094
2095 /** This initializes the neighbor information for an octahedron. The eight
2096  * initial faces are assigned ID numbers from -1 to -8. */
2097 inline void neighbor_track::init_octahedron() {
2098         int *q;
2099         q=mne[4];
2100         q[0]=-5;q[1]=-6;q[2]=-7;q[3]=-8;
2101         q[4]=-1;q[5]=-2;q[6]=-3;q[7]=-4;
2102         q[8]=-6;q[9]=-5;q[10]=-2;q[11]=-1;
2103         q[12]=-8;q[13]=-7;q[14]=-4;q[15]=-3;
2104         q[16]=-5;q[17]=-8;q[18]=-3;q[19]=-2;
2105         q[20]=-7;q[21]=-6;q[22]=-1;q[23]=-4;
2106         ne[0]=q;ne[1]=q+4;ne[2]=q+8;ne[3]=q+12;ne[4]=q+16;ne[5]=q+20;
2107 }
2108
2109 /** This initializes the neighbor information for an tetrahedron. The four
2110  * initial faces are assigned ID numbers from -1 to -4.*/
2111 inline void neighbor_track::init_tetrahedron() {
2112         int *q;
2113         q=mne[3];
2114         q[0]=-4;q[1]=-3;q[2]=-2;
2115         q[3]=-3;q[4]=-4;q[5]=-1;
2116         q[6]=-4;q[7]=-2;q[8]=-1;
2117         q[9]=-2;q[10]=-3;q[11]=-1;
2118         ne[0]=q;ne[1]=q+3;ne[2]=q+6;ne[3]=q+9;
2119 }
2120
2121 /** This is a basic operation to set a new pointer in the ne[] array.
2122  * \param[in] p the index in the ne[] array to set.
2123  * \param[in] n the order of the vertex. */
2124 inline void neighbor_track::set_pointer(int p,int n) {
2125         ne[p]=mne[n]+n*vc->mec[n];
2126 }
2127
2128 /** This is a basic operation to copy ne[c][d] to ne[a][b]. */
2129 inline void neighbor_track::copy(int a,int b,int c,int d) {
2130         ne[a][b]=ne[c][d];
2131 }
2132
2133 /** This is a basic operation to carry out ne[a][b]=c. */
2134 inline void neighbor_track::set(int a,int b,int c) {
2135         ne[a][b]=c;
2136 }
2137
2138 /** This is a basic operation to set the auxiliary pointer
2139  * paux1.
2140  * \param[in] k the order of the vertex to point to. */
2141 inline void neighbor_track::set_aux1(int k) {
2142         paux1=mne[k]+k*vc->mec[k];
2143 }
2144
2145 /** This is a basic operation to copy a neighbor into paux1.*/
2146 inline void neighbor_track::copy_aux1(int a,int b) {
2147         paux1[b]=ne[a][b];
2148 }
2149
2150 /** This is a basic operation to copy a neighbor into paux1 with a shift. It is
2151  * used in the delete_connection() routine of the voronoicell class. */
2152 inline void neighbor_track::copy_aux1_shift(int a,int b) {
2153         paux1[b]=ne[a][b+1];
2154 }
2155
2156 /** This routine sets the second auxiliary pointer to a new section
2157  * of memory, and then copies existing neighbor information into it. */
2158 inline void neighbor_track::set_aux2_copy(int a,int b) {
2159         paux2=mne[b]+b*vc->mec[b];
2160         for(int i=0;i<b;i++) ne[a][i]=paux2[i];
2161 }
2162
2163 /** This is a basic routine to copy ne[b] into ne[a]. */
2164 inline void neighbor_track::copy_pointer(int a,int b) {
2165         ne[a]=ne[b];
2166 }
2167
2168 /** This sets ne[j] to the first auxiliary pointer. */
2169 inline void neighbor_track::set_to_aux1(int j) {
2170         ne[j]=paux1;
2171 }
2172
2173 /** This sets ne[j] to the second auxiliary pointer. */
2174 inline void neighbor_track::set_to_aux2(int j) {
2175         ne[j]=paux2;
2176 }
2177
2178 /** This prints out the neighbor information for vertex i. */
2179 inline void neighbor_track::print_edges(int i) {
2180         cout << "    (";
2181         for(int j=0;j<vc->nu[i];j++) {
2182                 cout << ne[i][j] << (j==vc->nu[i]-1?")":",");
2183         }
2184 }
2185
2186 /** This allocates a new array and sets the auxiliary pointer
2187  * to it. */
2188 inline void neighbor_track::allocate_aux1(int i) {
2189         paux1=new int[i*vc->mem[i]];
2190 }
2191
2192 /** This deletes a particular neighbor array and switches the
2193  * pointer to the auxiliary pointer. */
2194 inline void neighbor_track::switch_to_aux1(int i) {
2195         delete [] mne[i];
2196         mne[i]=paux1;
2197 }
2198
2199 /** This routine copies neighbor information into the
2200  * auxiliary pointer. */
2201 inline void neighbor_track::copy_to_aux1(int i,int m) {
2202         paux1[m]=mne[i][m];
2203 }
2204
2205 /** This sets ne[k] to the auxiliary pointer with an offset. */
2206 inline void neighbor_track::set_to_aux1_offset(int k,int m) {
2207         ne[k]=paux1+m;
2208 }
2209
2210 /** This routine checks to make sure the neighbor information of each facets is
2211  * consistent.*/
2212 void neighbor_track::check_facets() {
2213         int **edp,*nup;edp=vc->ed;nup=vc->nu;
2214         int i,j,k,l,m,q;
2215         for(i=0;i<vc->p;i++) {
2216                 for(j=0;j<nup[i];j++) {
2217                         k=edp[i][j];
2218                         if (k>=0) {
2219                                 edp[i][j]=-1-k;
2220                                 q=ne[i][j];
2221                                 l=vc->cycle_up(edp[i][nup[i]+j],k);
2222                                 do {
2223                                         m=edp[k][l];
2224                                         edp[k][l]=-1-m;
2225                                         if (ne[k][l]!=q) cerr << "Facet error at (" << k << "," << l << ")=" << ne[k][l] << ", started from (" << i << "," << j << ")=" << q << endl;
2226                                         l=vc->cycle_up(edp[k][nup[k]+l],m);
2227                                         k=m;
2228                                 } while (k!=i);
2229                         }
2230                 }
2231         }
2232         for(i=0;i<vc->p;i++) {
2233                 for(j=0;j<nup[i];j++) {
2234                         if(edp[i][j]>=0) throw fatal_error("Facet labeling routine didn't look everywhere");
2235                         edp[i][j]=-1-edp[i][j];
2236                 }
2237         }
2238 }
2239
2240 /** This routine provides a list of plane IDs.
2241  * \param[in] &os An output stream to write to. */
2242 void neighbor_track::neighbors(ostream &os) {
2243         int **edp,*nup;edp=vc->ed;nup=vc->nu;
2244         int i,j,k,l,m;
2245         for(i=0;i<vc->p;i++) {
2246                 for(j=0;j<nup[i];j++) {
2247                         k=edp[i][j];
2248                         if (k>=0) {
2249                                 os << " " << ne[i][j];
2250                                 edp[i][j]=-1-k;
2251                                 l=vc->cycle_up(edp[i][nup[i]+j],k);
2252                                 do {
2253                                         m=edp[k][l];
2254                                         edp[k][l]=-1-m;
2255                                         l=vc->cycle_up(edp[k][nup[k]+l],m);
2256                                         k=m;
2257                                 } while (k!=i);
2258                         }
2259                 }
2260         }
2261         for(i=0;i<vc->p;i++) {
2262                 for(j=0;j<nup[i];j++) {
2263                         if(edp[i][j]>=0) throw fatal_error("Neighbor routine didn't look everywhere");
2264                         edp[i][j]=-1-edp[i][j];
2265                 }
2266         }
2267 }
2268
2269 /** This routine labels the facets in an arbitrary order, starting from one. */
2270 void neighbor_track::label_facets() {
2271         int **edp,*nup;edp=vc->ed;nup=vc->nu;
2272         int i,j,k,l,m,q=1;
2273         for(i=0;i<vc->p;i++) {
2274                 for(j=0;j<nup[i];j++) {
2275                         k=edp[i][j];
2276                         if (k>=0) {
2277                                 edp[i][j]=-1-k;
2278                                 ne[i][j]=q;
2279                                 l=vc->cycle_up(edp[i][nup[i]+j],k);
2280                                 do {
2281                                         m=edp[k][l];
2282                                         edp[k][l]=-1-m;
2283                                         ne[k][l]=q;
2284                                         l=vc->cycle_up(edp[k][nup[k]+l],m);
2285                                         k=m;
2286                                 } while (k!=i);
2287                                 q++;
2288                         }
2289                 }
2290         }
2291         for(i=0;i<vc->p;i++) {
2292                 for(j=0;j<nup[i];j++) {
2293                         if(edp[i][j]>=0) throw fatal_error("Facet labeling routine didn't look everywhere");
2294                         edp[i][j]=-1-edp[i][j];
2295                 }
2296         }
2297 }
2298
2299 /** This routine prints out a bracketed pair showing a vertex number, and the
2300  * corresponding neighbor information.
2301  * \param[in] &os The output stream to write to.
2302  * \param[in] i The vertex number to print.
2303  * \param[in] j The index of the neighbor information to print. */
2304 void neighbor_track::print(ostream &os,int i,int j) {
2305         os << "(" << i << "," << ne[i][j] << ")";
2306 }