| blob | d0340fa411ce8447ba0a82ca4b2cb9b366a35ab1 |
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 areas of each face of the Voronoi cell and prints the
1614 * results to an output stream.
1615 * \param[out] v the vector to store the results in. */
1642 }
1644 }
1645 }
1647 }
1650 /** Calculates the total surface area of the Voronoi cell.
1651 * \return The computed area. */
1676 }
1677 }
1678 }
1681 }
1684 /** Calculates the centroid of the Voronoi cell, by decomposing the cell into
1685 * tetrahedra extending outward from the zeroth vertex.
1686 * \param[out] (cx,cy,cz) references to floating point numbers in which to
1687 * pass back the centroid vector. */
1717 }
1718 }
1719 }
1720 }
1728 }
1730 /** Computes the maximum radius squared of a vertex from the center of the
1731 * cell. It can be used to determine when enough particles have been testing an
1732 * all planes that could cut the cell have been considered.
1733 * \return The maximum radius squared of a vertex.*/
1742 }
1744 }
1746 /** Calculates the total edge distance of the Voronoi cell.
1747 * \return A floating point number holding the calculated distance. */
1758 }
1759 }
1761 }
1763 /** Outputs the edges of the Voronoi cell in POV-Ray format to an open file
1764 * stream, displacing the cell by given vector.
1765 * \param[in] (x,y,z) a displacement vector to be added to the cell's position.
1766 * \param[in] fp a file handle to write to. */
1779 }
1780 }
1781 }
1782 }
1784 /** Outputs the edges of the Voronoi cell in gnuplot format to an output stream.
1785 * \param[in] (x,y,z) a displacement vector to be added to the cell's position.
1786 * \param[in] fp a file handle to write to. */
1801 }
1802 }
1804 }
1810 }
1812 }
1814 /** Outputs the Voronoi cell in the POV mesh2 format, described in section
1815 * 1.3.2.2 of the POV-Ray documentation. The mesh2 output consists of a list of
1816 * vertex vectors, followed by a list of triangular faces. The routine also
1817 * makes use of the optional inside_vector specification, which makes the mesh
1818 * object solid, so that the POV-Ray Constructive Solid Geometry (CSG) can be
1819 * applied.
1820 * \param[in] (x,y,z) a displacement vector to be added to the cell's position.
1821 * \param[in] fp a file handle to write to. */
1839 }
1840 }
1841 }
1844 }
1846 /** Several routines in the class that gather cell-based statistics internally
1847 * track their progress by flipping edges to negative so that they know what
1848 * parts of the cell have already been tested. This function resets them back
1849 * to positive. When it is called, it assumes that every edge in the routine
1850 * should have already been flipped to negative, and it bails out with an
1851 * internal error if it encounters a positive edge. */
1855 if(ed[i][j]>=0) voro_fatal_error("Edge reset routine found a previously untested edge",VOROPP_INTERNAL_ERROR);
1857 }
1858 }
1860 /** Checks to see if a given vertex is inside, outside or within the test
1861 * plane. If the point is far away from the test plane, the routine immediately
1862 * returns whether it is inside or outside. If the routine is close the the
1863 * plane and within the specified tolerance, then the special check_marginal()
1864 * routine is called.
1865 * \param[in] n the vertex to test.
1866 * \param[out] ans the result of the scalar product used in evaluating the
1867 * location of the point.
1868 * \return -1 if the point is inside the plane, 1 if the point is outside the
1869 * plane, or 0 if the point is within the plane. */
1875 }
1886 }
1889 /** Checks to see if a given vertex is inside, outside or within the test
1890 * plane. If the point is far away from the test plane, the routine immediately
1891 * returns whether it is inside or outside. If the routine is close the the
1892 * plane and within the specified tolerance, then the special check_marginal()
1893 * routine is called.
1894 * \param[in] n the vertex to test.
1895 * \param[out] ans the result of the scalar product used in evaluating the
1896 * location of the point.
1897 * \return -1 if the point is inside the plane, 1 if the point is outside the
1898 * plane, or 0 if the point is within the plane. */
1909 }
1911 }
1913 /** This routine calculates the unit normal vectors for every face.
1914 * \param[out] v the vector to store the results in. */
1921 }
1923 }
1925 /** This inline routine is called by normals(). It attempts to construct a
1926 * single normal vector that is associated with a particular face. It first
1927 * traces around the face, trying to find two vectors along the face edges
1928 * whose vector product is above the numerical tolerance. It then constructs
1929 * the normal vector using this product. If the face is too small, and none of
1930 * the vector products are large enough, the routine may return (0,0,0) as the
1931 * normal vector.
1932 * \param[in] v the vector to store the results in.
1933 * \param[in] i the initial vertex of the face to test.
1934 * \param[in] j the index of an edge of the vertex.
1935 * \param[in] k the neighboring vertex of i, set to ed[i][j]. */
1946 // Test to see if the length of this edge is above the tolerance
1955 // Construct the vector product of this edge with
1956 // the previous one
1962 // Test to see if this vector product of the
1963 // two edges is above the tolerance
1966 // Construct the normal vector and print it
1972 // Mark all of the remaining edges of this
1973 // face and exit
1977 }
1979 }
1980 }
1985 }
1992 }
1995 /** Returns the number of faces of a computed Voronoi cell.
1996 * \return The number of faces. */
2002 s++;
2012 }
2013 }
2016 }
2018 /** Returns a vector of the vertex orders.
2019 * \param[out] v the vector to store the results in. */
2023 }
2025 /** Outputs the vertex orders.
2026 * \param[out] fp the file handle to write to. */
2031 }
2032 }
2034 /** Returns a vector of the vertex vectors using the local coordinate system.
2035 * \param[out] v the vector to store the results in. */
2043 }
2044 }
2046 /** Outputs the vertex vectors using the local coordinate system.
2047 * \param[out] fp the file handle to write to. */
2051 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);
2052 }
2053 }
2056 /** Returns a vector of the vertex vectors in the global coordinate system.
2057 * \param[out] v the vector to store the results in.
2058 * \param[in] (x,y,z) the position vector of the particle in the global
2059 * coordinate system. */
2067 }
2068 }
2070 /** Outputs the vertex vectors using the global coordinate system.
2071 * \param[out] fp the file handle to write to.
2072 * \param[in] (x,y,z) the position vector of the particle in the global
2073 * coordinate system. */
2077 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);
2078 }
2079 }
2081 /** This routine returns the perimeters of each face.
2082 * \param[out] v the vector to store the results in. */
2107 }
2108 }
2110 }
2112 /** For each face, this routine outputs a bracketed sequence of numbers
2113 * containing a list of all the vertices that make up that face.
2114 * \param[out] v the vector to store the results in. */
2135 }
2136 }
2138 }
2140 /** Outputs a list of the number of edges in each face.
2141 * \param[out] v the vector to store the results in. */
2152 q++;
2159 }
2160 }
2162 }
2164 /** Computes the number of edges that each face has and outputs a frequency
2165 * table of the results.
2166 * \param[out] v the vector to store the results in. */
2177 q++;
2185 }
2186 }
2188 }
2190 /** This routine tests to see whether the cell intersects a plane by starting
2191 * from the guess point up. If up intersects, then it immediately returns true.
2192 * Otherwise, it calls the plane_intersects_track() routine.
2193 * \param[in] (x,y,z) the normal vector to the plane.
2194 * \param[in] rsq the distance along this vector of the plane.
2195 * \return False if the plane does not intersect the plane, true if it does. */
2200 }
2202 /** This routine tests to see if a cell intersects a plane. It first tests a
2203 * random sample of approximately sqrt(p)/4 points. If any of those are
2204 * intersect, then it immediately returns true. Otherwise, it takes the closest
2205 * point and passes that to plane_intersect_track() routine.
2206 * \param[in] (x,y,z) the normal vector to the plane.
2207 * \param[in] rsq the distance along this vector of the plane.
2208 * \return False if the plane does not intersect the plane, true if it does. */
2220 }
2222 }
2224 }
2226 }
2228 /* This routine tests to see if a cell intersects a plane, by tracing over the
2229 * cell from vertex to vertex, starting at up. It is meant to be called either
2230 * by plane_intersects() or plane_intersects_track(), when those routines
2231 * cannot immediately resolve the case.
2232 * \param[in] (x,y,z) the normal vector to the plane.
2233 * \param[in] rsq the distance along this vector of the plane.
2234 * \param[in] g the distance of up from the plane.
2235 * \return False if the plane does not intersect the plane, true if it does. */
2236 inline bool voronoicell_base::plane_intersects_track(double x,double y,double z,double rsq,double g) {
2240 /*
2241 int ls,us,lp;
2242 double l,u;
2243 unsigned int uw;
2245 // Initialize the safe testing routine
2246 px=x;py=y;pz=z;prsq=rsq;
2247 maskc+=4;
2248 if(maskc<4) reset_mask();
2250 return search_upward(uw,lp,ls,us,l,u);
2251 }*/
2252 /*
2253 int count=0,ls,us,tp;
2254 double t;
2255 // The test point is outside of the cutting space
2256 for(us=0;us<nu[up];us++) {
2257 tp=ed[up][us];
2258 t=x*pts[tp<<2]+y*pts[(tp<<2)+1]+z*pts[(tp<<2)+2];
2259 if(t>g) {
2260 ls=ed[up][nu[up]+us];
2261 up=tp;
2262 while (t<rsq) {
2263 if(++count>=p) {
2264 #if VOROPP_VERBOSE >=1
2265 fputs("Bailed out of convex calculation",stderr);
2266 #endif
2267 for(tp=0;tp<p;tp++) if(x*pts[tp<<2]+y*pts[(tp<<2)+1]+z*pts[(tp<<2)+2]>rsq) return true;
2268 return false;
2269 }
2271 // Test all the neighbors of the current point
2272 // and find the one which is closest to the
2273 // plane
2274 for(us=0;us<ls;us++) {
2275 tp=ed[up][us];double *pp=pts+(tp<<2);
2276 g=x*(*pp)+y*pp[1]+z*pp[2];
2277 if(g>t) break;
2278 }
2279 if(us==ls) {
2280 us++;
2281 while(us<nu[up]) {
2282 tp=ed[up][us];double *pp=pts+(tp<<2);
2283 g=x*(*pp)+y*pp[1]+z*pp[2];
2284 if(g>t) break;
2285 us++;
2286 }
2287 if(us==nu[up]) return false;
2288 }
2289 ls=ed[up][nu[up]+us];up=tp;t=g;
2290 }
2291 return true;
2292 }
2293 }
2294 return false;*/
2295 }
2297 /** Counts the number of edges of the Voronoi cell.
2298 * \return the number of edges. */
2303 }
2305 /** Outputs a custom string of information about the Voronoi cell. The string
2306 * of information follows a similar style as the C printf command, and detailed
2307 * information about its format is available at
2308 * http://math.lbl.gov/voro++/doc/custom.html.
2309 * \param[in] format the custom string to print.
2310 * \param[in] i the ID of the particle associated with this Voronoi cell.
2311 * \param[in] (x,y,z) the position of the particle associated with this Voronoi
2312 * cell.
2313 * \param[in] r a radius associated with the particle.
2314 * \param[in] fp the file handle to write to. */
2315 void voronoicell_base::output_custom(const char *format,int i,double x,double y,double z,double r,FILE *fp) {
2321 fmp++;
2324 // Particle-related output
2332 // Vertex-related output
2339 // Edge-related output
2344 // Face-related output
2364 // Volume-related output
2377 // End-of-string reached
2380 // The percent sign is not part of a
2381 // control sequence
2383 }
2385 fmp++;
2386 }
2388 }
2390 /** This initializes the class to be a rectangular box. It calls the base class
2391 * initialization routine to set up the edge and vertex information, and then
2392 * sets up the neighbor information, with initial faces being assigned ID
2393 * numbers from -1 to -6.
2394 * \param[in] (xmin,xmax) the minimum and maximum x coordinates.
2395 * \param[in] (ymin,ymax) the minimum and maximum y coordinates.
2396 * \param[in] (zmin,zmax) the minimum and maximum z coordinates. */
2397 void voronoicell_neighbor::init(double xmin,double xmax,double ymin,double ymax,double zmin,double zmax) {
2410 }
2412 /** This initializes the class to be an octahedron. It calls the base class
2413 * initialization routine to set up the edge and vertex information, and then
2414 * sets up the neighbor information, with the initial faces being assigned ID
2415 * numbers from -1 to -8.
2416 * \param[in] l The distance from the octahedron center to a vertex. Six
2417 * vertices are initialized at (-l,0,0), (l,0,0), (0,-l,0),
2418 * (0,l,0), (0,0,-l), and (0,0,l). */
2429 }
2431 /** This initializes the class to be a tetrahedron. It calls the base class
2432 * initialization routine to set up the edge and vertex information, and then
2433 * sets up the neighbor information, with the initial faces being assigned ID
2434 * numbers from -1 to -4.
2435 * \param (x0,y0,z0) a position vector for the first vertex.
2436 * \param (x1,y1,z1) a position vector for the second vertex.
2437 * \param (x2,y2,z2) a position vector for the third vertex.
2438 * \param (x3,y3,z3) a position vector for the fourth vertex. */
2439 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) {
2447 }
2449 /** This routine checks to make sure the neighbor information of each face is
2450 * consistent. */
2462 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);
2466 }
2467 }
2469 }
2471 /** The class constructor allocates memory for storing neighbor information. */
2479 }
2481 /** The class destructor frees the dynamically allocated memory for storing
2482 * neighbor information. */
2487 }
2489 /** Computes a vector list of neighbors. */
2505 }
2506 }
2508 }
2510 /** Prints the vertices, their edges, the relation table, and also notifies if
2511 * any memory errors are visible. */
2525 }
2526 }
2528 /** This prints out the neighbor information for vertex i. */
2536 }
2538 // Explicit instantiation
2542 template void voronoicell_base::check_memory_for_copy(voronoicell_neighbor&,voronoicell_base*);
2544 }