| blob | dfa3d92d7e536a5b8a1814693e4eb66620716c33 |
1 // Voro++, a 3D cell-based Voronoi library
2 //
3 // Author : Chris H. Rycroft (LBL / UC Berkeley)
4 // Email : chr@alum.mit.edu
5 // Date : August 30th 2011
7 /** \file cell.cc
8 * \brief Function implementations for the voronoicell and related classes. */
10 #include <cmath>
11 #include <cstring>
19 /** Constructs a Voronoi cell and sets up the initial memory. */
38 }
44 }
45 }
47 /** The voronoicell destructor deallocates all the dynamic memory. */
56 }
58 /** Ensures that enough memory is allocated prior to carrying out a copy.
59 * \param[in] vc a reference to the specialized version of the calling class.
60 * \param[in] vb a pointered to the class to be copied. */
66 }
68 /** Copies the vertex and edge information from another class. The routine
69 * assumes that enough memory is available for the copy.
70 * \param[in] vb a pointer to the class to copy. */
78 }
81 }
83 /** Copies the information from another voronoicell class into this
84 * class, extending memory allocation if necessary.
85 * \param[in] c the class to copy. */
93 }
94 }
96 /** Copies the information from another voronoicell_neighbor class into this
97 * class, extending memory allocation if necessary.
98 * \param[in] c the class to copy. */
106 }
107 }
109 /** Translates the vertices of the Voronoi cell by a given vector.
110 * \param[in] (x,y,z) the coordinates of the vector. */
116 }
117 }
119 /** Increases the memory storage for a particular vertex order, by increasing
120 * the size of the of the corresponding mep array. If the arrays already exist,
121 * their size is doubled; if they don't exist, then new ones of size
122 * init_n_vertices are allocated. The routine also ensures that the pointers in
123 * the ed array are updated, by making use of the back pointers. For the cases
124 * where the back pointer has been temporarily overwritten in the marginal
125 * vertex code, the auxiliary delete stack is scanned to find out how to update
126 * the ed value. If the template has been instantiated with the neighbor
127 * tracking turned on, then the routine also reallocates the corresponding mne
128 * array.
129 * \param[in] i the order of the vertex memory to be increased. */
137 #if VOROPP_VERBOSE >=2
139 #endif
143 if(mem[i]>max_n_vertices) voro_fatal_error("Point memory allocation exceeded absolute maximum",VOROPP_MEMORY_ERROR);
144 #if VOROPP_VERBOSE >=2
146 #endif
162 }
163 }
170 }
171 }
173 }
174 #if VOROPP_VERBOSE >=3
176 #endif
177 }
180 }
184 }
185 }
187 /** Doubles the maximum number of vertices allowed, by reallocating the ed, nu,
188 * and pts arrays. If the allocation exceeds the absolute maximum set in
189 * max_vertices, then the routine exits with a fatal error. If the template has
190 * been instantiated with the neighbor tracking turned on, then the routine
191 * also reallocates the ne array. */
196 if(i>max_vertices) voro_fatal_error("Vertex memory allocation exceeded absolute maximum",VOROPP_MEMORY_ERROR);
197 #if VOROPP_VERBOSE >=2
199 #endif
216 }
218 /** Doubles the maximum allowed vertex order, by reallocating mem, mep, and mec
219 * arrays. If the allocation exceeds the absolute maximum set in
220 * max_vertex_order, then the routine causes a fatal error. If the template has
221 * been instantiated with the neighbor tracking turned on, then the routine
222 * also reallocates the mne array. */
226 if(i>max_vertex_order) voro_fatal_error("Vertex order memory allocation exceeded absolute maximum",VOROPP_MEMORY_ERROR);
227 #if VOROPP_VERBOSE >=2
229 #endif
241 }
243 /** Doubles the size allocation of the main delete stack. If the allocation
244 * exceeds the absolute maximum set in max_delete_size, then routine causes a
245 * fatal error. */
248 if(current_delete_size>max_delete_size) voro_fatal_error("Delete stack 1 memory allocation exceeded absolute maximum",VOROPP_MEMORY_ERROR);
249 #if VOROPP_VERBOSE >=2
251 #endif
256 }
258 /** Doubles the size allocation of the auxiliary delete stack. If the
259 * allocation exceeds the absolute maximum set in max_delete2_size, then the
260 * routine causes a fatal error. */
263 if(current_delete2_size>max_delete2_size) voro_fatal_error("Delete stack 2 memory allocation exceeded absolute maximum",VOROPP_MEMORY_ERROR);
264 #if VOROPP_VERBOSE >=2
266 #endif
271 }
273 /** Doubles the size allocation of the auxiliary delete stack. If the
274 * allocation exceeds the absolute maximum set in max_delete2_size, then the
275 * routine causes a fatal error. */
278 if(current_xsearch_size>max_xsearch_size) voro_fatal_error("Extra search stack memory allocation exceeded absolute maximum",VOROPP_MEMORY_ERROR);
279 #if VOROPP_VERBOSE >=2
281 #endif
286 }
289 /** Initializes a Voronoi cell as a rectangular box with the given dimensions.
290 * \param[in] (xmin,xmax) the minimum and maximum x coordinates.
291 * \param[in] (ymin,ymax) the minimum and maximum y coordinates.
292 * \param[in] (zmin,zmax) the minimum and maximum z coordinates. */
293 void voronoicell_base::init_base(double xmin,double xmax,double ymin,double ymax,double zmin,double zmax) {
316 }
318 /** Initializes an L-shaped Voronoi cell of a fixed size for testing the
319 * convexity robustness. */
353 }
355 /** Initializes a Voronoi cell as a regular octahedron.
356 * \param[in] l The distance from the octahedron center to a vertex. Six
357 * vertices are initialized at (-l,0,0), (l,0,0), (0,-l,0),
358 * (0,l,0), (0,0,-l), and (0,0,l). */
377 }
379 /** Initializes a Voronoi cell as a tetrahedron. It assumes that the normal to
380 * the face for the first three vertices points inside.
381 * \param (x0,y0,z0) a position vector for the first vertex.
382 * \param (x1,y1,z1) a position vector for the second vertex.
383 * \param (x2,y2,z2) a position vector for the third vertex.
384 * \param (x3,y3,z3) a position vector for the fourth vertex. */
385 void voronoicell_base::init_tetrahedron_base(double x0,double y0,double z0,double x1,double y1,double z1,double x2,double y2,double z2,double x3,double y3,double z3) {
399 }
401 /** Checks that the relational table of the Voronoi cell is accurate, and
402 * prints out any errors. This algorithm is O(p), so running it every time the
403 * plane routine is called will result in a significant slowdown. */
408 }
410 /** This routine checks for any two vertices that are connected by more than
411 * one edge. The plane algorithm is designed so that this should not happen, so
412 * any occurrences are most likely errors. Note that the routine is O(p), so
413 * running it every time the plane routine is called will result in a
414 * significant slowdown. */
419 }
421 /** Constructs the relational table if the edges have been specified. */
428 l++;
430 }
432 }
433 }
435 /** Starting from a point within the current cutting plane, this routine attempts
436 * to find an edge to a point outside the cutting plane. This prevents the plane
437 * routine from .
438 * \param[in,out] up */
451 }
453 }
454 }
457 }
459 /** Adds a point to the auxiliary delete stack if it is not already there.
460 * \param[in] vc a reference to the specialized version of the calling class.
461 * \param[in] lp the index of the point to add.
462 * \param[in,out] stackp2 a pointer to the end of the stack entries. */
467 }
469 /** Assuming that the point up is outside the cutting plane, this routine
470 * searches upwards along edges trying to find an edge that intersects the
471 * cutting plane.
472 * \param[in] rsq the distance along this vector of the plane.
473 * \param[in,out] u the dot product of point up with the normal.
474 * \return True if the cutting plane was reached, false otherwise. */
475 inline bool voronoicell_base::search_upward(unsigned int &uw,int &lp,int &ls,int &us,double &l,double &u) {
479 // The test point is outside of the cutting space
484 }
488 }
491 //if(++count>=p) failsafe_find(lp,ls,us,l,u);
493 // Test all the neighbors of the current point
494 // and find the one which is closest to the
495 // plane
502 }
506 }
507 }
510 }
512 /** Checks whether a particular point lp is a definite maximum, searching
513 * through any possible minor non-convexities, for a better maximum.
514 * \param[in] (x,y,z) the normal vector to the plane. */
520 // Check to see whether point up is a well-defined maximum. Otherwise
521 // any neighboring vertices of up that are marginal need to be
522 // followed, to see if they lead to a better maximum.
527 }
530 // The point tp is marginal, so it will be necessary to do the
531 // flood-fill search. Mark the point tp and the point qp, and search
532 // any remaining neighbors of the point tp.
537 ts++;
545 }
546 ts++;
547 }
549 // Consider additional marginal points, starting with the original
550 // point qp
557 // Skip the point if it's already marked
561 // This point is a better maximum. Reset markers and
562 // return true.
573 }
575 // The point is marginal and therefore must also be
576 // considered
582 }
585 }
586 }
587 }
589 // Reset markers and return false
593 }
595 inline bool voronoicell_base::search_downward(unsigned int &lw,int &lp,int &ls,int &us,double &l,double &u) {
598 // The test point is outside of the cutting space
603 }
607 //if(++count>=p) failsafe_find(lp,ls,us,l,u);
609 // Test all the neighbors of the current point
610 // and find the one which is closest to the
611 // plane
618 }
620 }
623 }
630 // Check to see whether point up is a well-defined maximum. Otherwise
631 // any neighboring vertices of up that are marginal need to be
632 // followed, to see if they lead to a better maximum.
637 }
640 // The point tp is marginal, so it will be necessary to do the
641 // flood-fill search. Mark the point tp and the point qp, and search
642 // any remaining neighbors of the point tp.
647 ts++;
655 }
656 ts++;
657 }
659 // Consider additional marginal points, starting with the original
660 // point qp
667 // Skip the point if it's already marked
671 // This point is a better minimum. Reset markers and
672 // return true.
683 }
685 // The point is marginal and therefore must also be
686 // considered
692 }
695 }
696 }
697 }
699 // Reset markers and return false
703 }
705 /** Cuts the Voronoi cell by a particle whose center is at a separation of
706 * (x,y,z) from the cell center. The value of rsq should be initially set to
707 * \f$x^2+y^2+z^2\f$.
708 * \param[in] vc a reference to the specialized version of the calling class.
709 * \param[in] (x,y,z) the normal vector to the plane.
710 * \param[in] rsq the distance along this vector of the plane.
711 * \param[in] p_id the plane ID (for neighbor tracking only).
712 * \return False if the plane cut deleted the cell entirely, true otherwise. */
721 // Initialize the safe testing routine
735 }
737 // Set stack pointers
740 // Store initial number of vertices
751 // Skip if this is a new vertex
760 }
763 // This is a possible facet starting
764 // from a vertex on the cutting plane
768 // This is a new facet
772 }
773 }
774 xtra++;
775 }
777 // Reset back pointers on extra search stack
781 }
783 // Delete points: first, remove any duplicates
789 dsp++;
791 }
793 // Add the points in the auxiliary delete stack,
794 // and reset their back pointers
802 }
803 }
805 // Scan connections and add in extras
815 }
818 }
819 }
820 }
823 // Delete them from the array structure
834 }
838 // Vertex management
843 // Memory management
852 // Edge management
858 }
860 // Check for any vertices of zero order
863 // Collapse any order 2 vertices and exit
865 }
867 /** Creates a new facet.
868 * \return True if cell deleted, false otherwise. */
870 bool voronoicell_base::create_facet(vc_class &vc,int lp,int ls,double l,int us,double u,int p_id) {
876 // We're about to add the first point of the new facet. In either
877 // routine, we have to add a point, so first check there's space for
878 // it.
883 // We want to be strict about reaching the conclusion that the
884 // cell is entirely within the cutting plane. It's not enough
885 // to find a vertex that has edges which are all inside or on
886 // the plane. If the vertex has neighbors that are also on the
887 // plane, we should check those too.
890 // The search algorithm found a point which is on the cutting
891 // plane. We leave that point in place, and create a new one at
892 // the same location.
897 // Search for a collection of edges of the test vertex which
898 // are outside of the cutting space. Begin by testing the
899 // zeroth edge.
905 // The first edge is either inside the cutting space,
906 // or lies within the cutting plane. Test the edges
907 // sequentially until we find one that is outside.
910 i++;
912 // If we reached the last edge with no luck
913 // then all of the vertices are inside
914 // or on the plane, so the cell is completely
915 // deleted
922 // We found an edge outside the cutting space. Keep
923 // moving through these edges until we find one that's
924 // inside or on the plane.
929 j++;
930 }
932 // Compute the number of edges for the new vertex. In
933 // general it will be the number of outside edges
934 // found, plus two. But we need to recognize the
935 // special case when all but one edge is outside, and
936 // the remaining one is on the plane. For that case we
937 // have to reduce the edge count by one to prevent
938 // doubling up.
945 // Add memory for the new vertex if needed, and
946 // initialize
953 // Copy the edges of the original vertex into the new
954 // one. Delete the edges of the original vertex, and
955 // update the relational table.
967 }
971 // In this case, the zeroth edge is outside the cutting
972 // plane. Begin by searching backwards from the last
973 // edge until we find an edge which isn't outside.
978 i--;
980 // If i reaches zero, then we have a point in
981 // the plane all of whose edges are outside
982 // the cutting space, so we just exit
986 }
988 // Now search forwards from zero
993 j++;
996 }
998 // Compute the number of edges for the new vertex. In
999 // general it will be the number of outside edges
1000 // found, plus two. But we need to recognize the
1001 // special case when all but one edge is outside, and
1002 // the remaining one is on the plane. For that case we
1003 // have to reduce the edge count by one to prevent
1004 // doubling up.
1010 }
1012 // Add memory to store the vertex if it doesn't exist
1013 // already
1018 // Copy the edges of the original vertex into the new
1019 // one. Delete the edges of the original vertex, and
1020 // update the relational table.
1035 }
1047 }
1049 }
1055 // Add this point to the auxiliary delete stack
1059 // Look at the edges on either side of the group that was
1060 // detected. We're going to commence facet computation by
1061 // moving along one of them. We are going to end up coming back
1062 // along the other one.
1072 // The search algorithm found an intersected edge between the
1073 // points lp and up. Create a new vertex between them which
1074 // lies on the cutting plane. Since u and l differ by at least
1075 // the tolerance, this division should never screw up.
1083 // This point will always have three edges. Connect one of them
1084 // to lp.
1100 // Set the direction to move in
1103 }
1105 // When the code reaches here, we have initialized the first point, and
1106 // we have a direction for moving it to construct the rest of the facet
1110 // We're currently tracing round an intersected facet. Keep
1111 // moving around it until we find a point or edge which
1112 // intersects the plane.
1118 // The point is still in the cutting space. Just add it
1119 // to the delete stack and keep moving.
1128 // The point is outside of the cutting space, so we've
1129 // found an intersected edge. Introduce a regular point
1130 // at the point of intersection. Connect it to the
1131 // point we just tested. Also connect it to the previous
1132 // new point in the facet we're constructing.
1161 // We've found a point which is on the cutting plane.
1162 // We're going to introduce a new point right here, but
1163 // first we need to figure out the number of edges it
1164 // has.
1167 // If the previous vertex detected a double edge, our
1168 // new vertex will have one less edge.
1174 // Start testing the edges of the current point until
1175 // we find one which isn't outside the cutting space
1177 k++;
1183 // Now we need to find out whether this marginal vertex
1184 // we are on has been visited before, because if that's
1185 // the case, we need to add vertices to the existing
1186 // new vertex, rather than creating a fresh one. We also
1187 // need to figure out whether we're in a case where we
1188 // might be creating a duplicate edge.
1192 // If we're heading into the final part of the
1193 // new facet, then we never worry about the
1194 // duplicate edge calculation.
1200 // This vertex was visited before, so
1201 // count those vertices to the ones we
1202 // already have.
1205 // The only time when we might make a
1206 // duplicate edge is if the point we're
1207 // going to move to next is also a
1208 // marginal point, so test for that
1209 // first.
1212 // Now see whether this marginal point
1213 // has been visited before.
1217 // Now see if the last edge of that other
1218 // marginal point actually ends up here.
1225 // That marginal point hasn't been visited
1226 // before, so we probably don't have to worry
1227 // about duplicate edges, except in the
1228 // case when that's the way into the end
1229 // of the facet, because that way always creates
1230 // an edge.
1235 }
1239 // The vertex hasn't been visited
1240 // before, but let's see if it's
1241 // marginal
1244 // If it is, we need to check
1245 // for the case that it's a
1246 // small branch, and that we're
1247 // heading right back to where
1248 // we came from
1255 }
1256 }
1258 // k now holds the number of edges of the new vertex
1259 // we are forming. Add memory for it if it doesn't exist
1260 // already.
1264 // Now create a new vertex with order k, or augment
1265 // the existing one
1268 // If we're augmenting a vertex but we don't
1269 // actually need any more edges, just skip this
1270 // routine to avoid memory confusion
1273 // Allocate memory and copy the edges
1274 // of the previous instance into it
1282 i++;
1283 }
1286 // Remove the previous instance with
1287 // fewer vertices from the memory
1288 // structure
1295 }
1301 // Allocate a new vertex of order k
1313 }
1316 // Unless the previous case was a double edge, connect
1317 // the first available edge of the new vertex to the
1318 // last one in the facet
1325 i++;
1326 }
1328 // Copy in the edges of the underlying vertex,
1329 // and do one less if this was a double edge
1340 i++;
1341 }
1347 // Update the double_edge flag, to pass it
1348 // to the next instance of this routine
1350 }
1351 }
1353 // Connect the final created vertex to the initial one
1359 }
1361 /** During the creation of a new facet in the plane routine, it is possible
1362 * that some order two vertices may arise. This routine removes them.
1363 * Suppose an order two vertex joins c and d. If there's a edge between
1364 * c and d already, then the order two vertex is just removed; otherwise,
1365 * the order two vertex is removed and c and d are joined together directly.
1366 * It is possible this process will create order two or order one vertices,
1367 * and the routine is continually run until all of them are removed.
1368 * \return False if the vertex removal was unsuccessful, indicative of the cell
1369 * reducing to zero volume and disappearing; true if the vertex removal
1370 * was successful. */
1377 // Pick a order 2 vertex and read in its edges
1381 #if VOROPP_VERBOSE >=1
1383 #endif
1385 }
1387 // Scan the edges of j to see if joins k
1390 }
1392 // If j doesn't already join k, join them together.
1393 // Otherwise delete the connection to the current
1394 // vertex from j and k.
1404 }
1406 // Compact the memory
1419 }
1421 // Collapse any order 1 vertices if they were created
1423 }
1425 }
1427 /** Order one vertices can potentially be created during the order two collapse
1428 * routine. This routine keeps removing them until there are none left.
1429 * \return False if the vertex removal was unsuccessful, indicative of the cell
1430 * having zero volume and disappearing; true if the vertex removal was
1431 * successful. */
1437 #if VOROPP_VERBOSE >=1
1439 #endif
1456 }
1457 }
1459 }
1461 /** This routine deletes the kth edge of vertex j and reorganizes the memory.
1462 * If the neighbor computation is enabled, we also have to supply an handedness
1463 * flag to decide whether to preserve the plane on the left or right of the
1464 * connection.
1465 * \return False if a zero order vertex was formed, indicative of the cell
1466 * disappearing; true if the vertex removal was successful. */
1471 #if VOROPP_VERBOSE >=1
1475 }
1476 #endif
1482 l++;
1483 }
1489 }
1496 l++;
1497 }
1508 }
1510 /** This routine is a fall-back, in case floating point errors caused the usual
1511 * search routine to fail. In the fall-back routine, we just test every edge to
1512 * find one straddling the plane. */
1516 /* qw=1;lw=0;
1517 for(qp=0;qp<p;qp++) {
1518 qw=m_test(qp,q);
1519 if(qw==1) {
1521 // The point is inside the cutting space. Now
1522 // see if we can find a neighbor which isn't.
1523 for(us=0;us<nu[qp];us++) {
1524 lp=ed[qp][us];
1525 if(lp<qp) {
1526 lw=m_test(lp,l);
1527 if(lw!=1) break;
1528 }
1529 }
1530 if(us<nu[qp]) {
1531 up=qp;
1532 if(lw==0) {
1533 complicated_setup=true;
1534 } else {
1535 complicated_setup=false;
1536 u=q;
1537 ls=ed[up][nu[up]+us];
1538 }
1539 break;
1540 }
1541 } else if(qw==-1) {
1543 // The point is outside the cutting space. See
1544 // if we can find a neighbor which isn't.
1545 for(ls=0;ls<nu[qp];ls++) {
1546 up=ed[qp][ls];
1547 if(up<qp) {
1548 uw=m_test(up,u);
1549 if(uw!=-1) break;
1550 }
1551 }
1552 if(ls<nu[qp]) {
1553 if(uw==0) {
1554 up=qp;
1555 complicated_setup=true;
1556 } else {
1557 complicated_setup=false;
1558 lp=qp;l=q;
1559 us=ed[lp][nu[lp]+ls];
1560 }
1561 break;
1562 }
1563 } else {
1565 // The point is in the plane, so we just
1566 // proceed with the complicated setup routine
1567 up=qp;
1568 complicated_setup=true;
1569 break;
1570 }
1571 }
1572 if(qp==p) return qw==-1?true:false;*/
1573 }
1575 /** Calculates the volume of the Voronoi cell, by decomposing the cell into
1576 * tetrahedra extending outward from the zeroth vertex, whose volumes are
1577 * evaluated using a scalar triple product.
1578 * \return A floating point number holding the calculated volume. */
1605 }
1606 }
1607 }
1608 }
1611 }
1613 /** Calculates the contributions to the Minkowski functionals for this Voronoi cell.
1614 * \param[in] r the radius to consider.
1615 * \param[out] ar the area functional.
1616 * \param[out] vo the volume functional. */
1631 }
1632 }
1633 }
1637 }
1639 inline void voronoicell_base::minkowski_contrib(int i,int k,int m,double r,double &ar,double &vo) {
1650 // Compute second orthonormal vector
1657 }
1661 // Compute third orthonormal vector
1679 }
1681 void voronoicell_base::minkowski_edge(double x0,double r1,double s1,double r2,double s2,double r,double &ar,double &vo) {
1691 }
1693 void voronoicell_base::minkowski_formula(double x0,double y0,double z0,double r,double &ar,double &vo) {
1700 double xs=x0*x0,ys=y0*y0,zs=z0*z0,res=xs+ys,rvs=res+zs,theta=atan(z0/y0),rs=r*r,rc=rs*r,temp,voc,arc;
1718 arc=x0*r*temp-0.5*temp2*y0*r/((rs-xs)*temp5)+x0*y0/6.*r/temp5+rs*0.5*temp3+rs*rs/3.*2*xs*ys/(res*(rs-xs)*sqrt((rs-xs)*(rs-xs)-(x2s-y2s-xs)*(x2s-y2s-xs)))-rs*0.5*temp4;
1722 }
1725 }
1728 /** Calculates the areas of each face of the Voronoi cell and prints the
1729 * results to an output stream.
1730 * \param[out] v the vector to store the results in. */
1757 }
1759 }
1760 }
1762 }
1765 /** Calculates the total surface area of the Voronoi cell.
1766 * \return The computed area. */
1791 }
1792 }
1793 }
1796 }
1799 /** Calculates the centroid of the Voronoi cell, by decomposing the cell into
1800 * tetrahedra extending outward from the zeroth vertex.
1801 * \param[out] (cx,cy,cz) references to floating point numbers in which to
1802 * pass back the centroid vector. */
1832 }
1833 }
1834 }
1835 }
1843 }
1845 /** Computes the maximum radius squared of a vertex from the center of the
1846 * cell. It can be used to determine when enough particles have been testing an
1847 * all planes that could cut the cell have been considered.
1848 * \return The maximum radius squared of a vertex.*/
1857 }
1859 }
1861 /** Calculates the total edge distance of the Voronoi cell.
1862 * \return A floating point number holding the calculated distance. */
1873 }
1874 }
1876 }
1878 /** Outputs the edges of the Voronoi cell in POV-Ray format to an open file
1879 * stream, displacing the cell by given vector.
1880 * \param[in] (x,y,z) a displacement vector to be added to the cell's position.
1881 * \param[in] fp a file handle to write to. */
1894 }
1895 }
1896 }
1897 }
1899 /** Outputs the edges of the Voronoi cell in gnuplot format to an output stream.
1900 * \param[in] (x,y,z) a displacement vector to be added to the cell's position.
1901 * \param[in] fp a file handle to write to. */
1916 }
1917 }
1919 }
1925 }
1927 }
1929 /** Outputs the Voronoi cell in the POV mesh2 format, described in section
1930 * 1.3.2.2 of the POV-Ray documentation. The mesh2 output consists of a list of
1931 * vertex vectors, followed by a list of triangular faces. The routine also
1932 * makes use of the optional inside_vector specification, which makes the mesh
1933 * object solid, so that the POV-Ray Constructive Solid Geometry (CSG) can be
1934 * applied.
1935 * \param[in] (x,y,z) a displacement vector to be added to the cell's position.
1936 * \param[in] fp a file handle to write to. */
1954 }
1955 }
1956 }
1959 }
1961 /** Several routines in the class that gather cell-based statistics internally
1962 * track their progress by flipping edges to negative so that they know what
1963 * parts of the cell have already been tested. This function resets them back
1964 * to positive. When it is called, it assumes that every edge in the routine
1965 * should have already been flipped to negative, and it bails out with an
1966 * internal error if it encounters a positive edge. */
1970 if(ed[i][j]>=0) voro_fatal_error("Edge reset routine found a previously untested edge",VOROPP_INTERNAL_ERROR);
1972 }
1973 }
1975 /** Checks to see if a given vertex is inside, outside or within the test
1976 * plane. If the point is far away from the test plane, the routine immediately
1977 * returns whether it is inside or outside. If the routine is close the the
1978 * plane and within the specified tolerance, then the special check_marginal()
1979 * routine is called.
1980 * \param[in] n the vertex to test.
1981 * \param[out] ans the result of the scalar product used in evaluating the
1982 * location of the point.
1983 * \return -1 if the point is inside the plane, 1 if the point is outside the
1984 * plane, or 0 if the point is within the plane. */
1990 }
2001 }
2004 /** Checks to see if a given vertex is inside, outside or within the test
2005 * plane. If the point is far away from the test plane, the routine immediately
2006 * returns whether it is inside or outside. If the routine is close the the
2007 * plane and within the specified tolerance, then the special check_marginal()
2008 * routine is called.
2009 * \param[in] n the vertex to test.
2010 * \param[out] ans the result of the scalar product used in evaluating the
2011 * location of the point.
2012 * \return -1 if the point is inside the plane, 1 if the point is outside the
2013 * plane, or 0 if the point is within the plane. */
2024 }
2026 }
2028 /** This routine calculates the unit normal vectors for every face.
2029 * \param[out] v the vector to store the results in. */
2036 }
2038 }
2040 /** This inline routine is called by normals(). It attempts to construct a
2041 * single normal vector that is associated with a particular face. It first
2042 * traces around the face, trying to find two vectors along the face edges
2043 * whose vector product is above the numerical tolerance. It then constructs
2044 * the normal vector using this product. If the face is too small, and none of
2045 * the vector products are large enough, the routine may return (0,0,0) as the
2046 * normal vector.
2047 * \param[in] v the vector to store the results in.
2048 * \param[in] i the initial vertex of the face to test.
2049 * \param[in] j the index of an edge of the vertex.
2050 * \param[in] k the neighboring vertex of i, set to ed[i][j]. */
2061 // Test to see if the length of this edge is above the tolerance
2070 // Construct the vector product of this edge with
2071 // the previous one
2077 // Test to see if this vector product of the
2078 // two edges is above the tolerance
2081 // Construct the normal vector and print it
2087 // Mark all of the remaining edges of this
2088 // face and exit
2092 }
2094 }
2095 }
2100 }
2107 }
2110 /** Returns the number of faces of a computed Voronoi cell.
2111 * \return The number of faces. */
2117 s++;
2127 }
2128 }
2131 }
2133 /** Returns a vector of the vertex orders.
2134 * \param[out] v the vector to store the results in. */
2138 }
2140 /** Outputs the vertex orders.
2141 * \param[out] fp the file handle to write to. */
2146 }
2147 }
2149 /** Returns a vector of the vertex vectors using the local coordinate system.
2150 * \param[out] v the vector to store the results in. */
2158 }
2159 }
2161 /** Outputs the vertex vectors using the local coordinate system.
2162 * \param[out] fp the file handle to write to. */
2166 for(double *ptsp=pts+4;ptsp<pts+(p<<2);ptsp+=4) fprintf(fp," (%g,%g,%g)",*ptsp*0.5,ptsp[1]*0.5,ptsp[2]*0.5);
2167 }
2168 }
2171 /** Returns a vector of the vertex vectors in the global coordinate system.
2172 * \param[out] v the vector to store the results in.
2173 * \param[in] (x,y,z) the position vector of the particle in the global
2174 * coordinate system. */
2182 }
2183 }
2185 /** Outputs the vertex vectors using the global coordinate system.
2186 * \param[out] fp the file handle to write to.
2187 * \param[in] (x,y,z) the position vector of the particle in the global
2188 * coordinate system. */
2192 for(double *ptsp=pts+4;ptsp<pts+(p<<2);ptsp+=4) fprintf(fp," (%g,%g,%g)",x+*ptsp*0.5,y+ptsp[1]*0.5,z+ptsp[2]*0.5);
2193 }
2194 }
2196 /** This routine returns the perimeters of each face.
2197 * \param[out] v the vector to store the results in. */
2222 }
2223 }
2225 }
2227 /** For each face, this routine outputs a bracketed sequence of numbers
2228 * containing a list of all the vertices that make up that face.
2229 * \param[out] v the vector to store the results in. */
2250 }
2251 }
2253 }
2255 /** Outputs a list of the number of edges in each face.
2256 * \param[out] v the vector to store the results in. */
2267 q++;
2274 }
2275 }
2277 }
2279 /** Computes the number of edges that each face has and outputs a frequency
2280 * table of the results.
2281 * \param[out] v the vector to store the results in. */
2292 q++;
2300 }
2301 }
2303 }
2305 /** This routine tests to see whether the cell intersects a plane by starting
2306 * from the guess point up. If up intersects, then it immediately returns true.
2307 * Otherwise, it calls the plane_intersects_track() routine.
2308 * \param[in] (x,y,z) the normal vector to the plane.
2309 * \param[in] rsq the distance along this vector of the plane.
2310 * \return False if the plane does not intersect the plane, true if it does. */
2315 }
2317 /** This routine tests to see if a cell intersects a plane. It first tests a
2318 * random sample of approximately sqrt(p)/4 points. If any of those are
2319 * intersect, then it immediately returns true. Otherwise, it takes the closest
2320 * point and passes that to plane_intersect_track() routine.
2321 * \param[in] (x,y,z) the normal vector to the plane.
2322 * \param[in] rsq the distance along this vector of the plane.
2323 * \return False if the plane does not intersect the plane, true if it does. */
2335 }
2337 }
2339 }
2341 }
2343 /* This routine tests to see if a cell intersects a plane, by tracing over the
2344 * cell from vertex to vertex, starting at up. It is meant to be called either
2345 * by plane_intersects() or plane_intersects_track(), when those routines
2346 * cannot immediately resolve the case.
2347 * \param[in] (x,y,z) the normal vector to the plane.
2348 * \param[in] rsq the distance along this vector of the plane.
2349 * \param[in] g the distance of up from the plane.
2350 * \return False if the plane does not intersect the plane, true if it does. */
2351 inline bool voronoicell_base::plane_intersects_track(double x,double y,double z,double rsq,double g) {
2355 /*
2356 int ls,us,lp;
2357 double l,u;
2358 unsigned int uw;
2360 // Initialize the safe testing routine
2361 px=x;py=y;pz=z;prsq=rsq;
2362 maskc+=4;
2363 if(maskc<4) reset_mask();
2365 return search_upward(uw,lp,ls,us,l,u);
2366 }*/
2367 /*
2368 int count=0,ls,us,tp;
2369 double t;
2370 // The test point is outside of the cutting space
2371 for(us=0;us<nu[up];us++) {
2372 tp=ed[up][us];
2373 t=x*pts[tp<<2]+y*pts[(tp<<2)+1]+z*pts[(tp<<2)+2];
2374 if(t>g) {
2375 ls=ed[up][nu[up]+us];
2376 up=tp;
2377 while (t<rsq) {
2378 if(++count>=p) {
2379 #if VOROPP_VERBOSE >=1
2380 fputs("Bailed out of convex calculation",stderr);
2381 #endif
2382 for(tp=0;tp<p;tp++) if(x*pts[tp<<2]+y*pts[(tp<<2)+1]+z*pts[(tp<<2)+2]>rsq) return true;
2383 return false;
2384 }
2386 // Test all the neighbors of the current point
2387 // and find the one which is closest to the
2388 // plane
2389 for(us=0;us<ls;us++) {
2390 tp=ed[up][us];double *pp=pts+(tp<<2);
2391 g=x*(*pp)+y*pp[1]+z*pp[2];
2392 if(g>t) break;
2393 }
2394 if(us==ls) {
2395 us++;
2396 while(us<nu[up]) {
2397 tp=ed[up][us];double *pp=pts+(tp<<2);
2398 g=x*(*pp)+y*pp[1]+z*pp[2];
2399 if(g>t) break;
2400 us++;
2401 }
2402 if(us==nu[up]) return false;
2403 }
2404 ls=ed[up][nu[up]+us];up=tp;t=g;
2405 }
2406 return true;
2407 }
2408 }
2409 return false;*/
2410 }
2412 /** Counts the number of edges of the Voronoi cell.
2413 * \return the number of edges. */
2418 }
2420 /** Outputs a custom string of information about the Voronoi cell. The string
2421 * of information follows a similar style as the C printf command, and detailed
2422 * information about its format is available at
2423 * http://math.lbl.gov/voro++/doc/custom.html.
2424 * \param[in] format the custom string to print.
2425 * \param[in] i the ID of the particle associated with this Voronoi cell.
2426 * \param[in] (x,y,z) the position of the particle associated with this Voronoi
2427 * cell.
2428 * \param[in] r a radius associated with the particle.
2429 * \param[in] fp the file handle to write to. */
2430 void voronoicell_base::output_custom(const char *format,int i,double x,double y,double z,double r,FILE *fp) {
2436 fmp++;
2439 // Particle-related output
2447 // Vertex-related output
2454 // Edge-related output
2459 // Face-related output
2479 // Volume-related output
2492 // End-of-string reached
2495 // The percent sign is not part of a
2496 // control sequence
2498 }
2500 fmp++;
2501 }
2503 }
2505 /** This initializes the class to be a rectangular box. It calls the base class
2506 * initialization routine to set up the edge and vertex information, and then
2507 * sets up the neighbor information, with initial faces being assigned ID
2508 * numbers from -1 to -6.
2509 * \param[in] (xmin,xmax) the minimum and maximum x coordinates.
2510 * \param[in] (ymin,ymax) the minimum and maximum y coordinates.
2511 * \param[in] (zmin,zmax) the minimum and maximum z coordinates. */
2512 void voronoicell_neighbor::init(double xmin,double xmax,double ymin,double ymax,double zmin,double zmax) {
2525 }
2527 /** This initializes the class to be an octahedron. It calls the base class
2528 * initialization routine to set up the edge and vertex information, and then
2529 * sets up the neighbor information, with the initial faces being assigned ID
2530 * numbers from -1 to -8.
2531 * \param[in] l The distance from the octahedron center to a vertex. Six
2532 * vertices are initialized at (-l,0,0), (l,0,0), (0,-l,0),
2533 * (0,l,0), (0,0,-l), and (0,0,l). */
2544 }
2546 /** This initializes the class to be a tetrahedron. It calls the base class
2547 * initialization routine to set up the edge and vertex information, and then
2548 * sets up the neighbor information, with the initial faces being assigned ID
2549 * numbers from -1 to -4.
2550 * \param (x0,y0,z0) a position vector for the first vertex.
2551 * \param (x1,y1,z1) a position vector for the second vertex.
2552 * \param (x2,y2,z2) a position vector for the third vertex.
2553 * \param (x3,y3,z3) a position vector for the fourth vertex. */
2554 void voronoicell_neighbor::init_tetrahedron(double x0,double y0,double z0,double x1,double y1,double z1,double x2,double y2,double z2,double x3,double y3,double z3) {
2562 }
2564 /** This routine checks to make sure the neighbor information of each face is
2565 * consistent. */
2577 if(ne[k][l]!=q) fprintf(stderr,"Facet error at (%d,%d)=%d, started from (%d,%d)=%d\n",k,l,ne[k][l],i,j,q);
2581 }
2582 }
2584 }
2586 /** The class constructor allocates memory for storing neighbor information. */
2594 }
2596 /** The class destructor frees the dynamically allocated memory for storing
2597 * neighbor information. */
2602 }
2604 /** Computes a vector list of neighbors. */
2620 }
2621 }
2623 }
2625 /** Prints the vertices, their edges, the relation table, and also notifies if
2626 * any memory errors are visible. */
2640 }
2641 }
2643 /** This prints out the neighbor information for vertex i. */
2651 }
2653 // Explicit instantiation
2657 template void voronoicell_base::check_memory_for_copy(voronoicell_neighbor&,voronoicell_base*);
2659 }