| blob | f480912ff590e1d8e1be2342df1c25c9e334e0e0 |
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. */
352 }
354 /** Initializes a Voronoi cell as a regular octahedron.
355 * \param[in] l The distance from the octahedron center to a vertex. Six
356 * vertices are initialized at (-l,0,0), (l,0,0), (0,-l,0),
357 * (0,l,0), (0,0,-l), and (0,0,l). */
376 }
378 /** Initializes a Voronoi cell as a tetrahedron. It assumes that the normal to
379 * the face for the first three vertices points inside.
380 * \param (x0,y0,z0) a position vector for the first vertex.
381 * \param (x1,y1,z1) a position vector for the second vertex.
382 * \param (x2,y2,z2) a position vector for the third vertex.
383 * \param (x3,y3,z3) a position vector for the fourth vertex. */
384 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) {
398 }
400 /** Checks that the relational table of the Voronoi cell is accurate, and
401 * prints out any errors. This algorithm is O(p), so running it every time the
402 * plane routine is called will result in a significant slowdown. */
407 }
409 /** This routine checks for any two vertices that are connected by more than
410 * one edge. The plane algorithm is designed so that this should not happen, so
411 * any occurrences are most likely errors. Note that the routine is O(p), so
412 * running it every time the plane routine is called will result in a
413 * significant slowdown. */
418 }
420 /** Constructs the relational table if the edges have been specified. */
427 l++;
429 }
431 }
432 }
434 /** Starting from a point within the current cutting plane, this routine attempts
435 * to find an edge to a point outside the cutting plane. This prevents the plane
436 * routine from .
437 * \param[in,out] up */
450 }
452 }
453 }
456 }
458 /** Adds a point to the auxiliary delete stack if it is not already there.
459 * \param[in] vc a reference to the specialized version of the calling class.
460 * \param[in] lp the index of the point to add.
461 * \param[in,out] stackp2 a pointer to the end of the stack entries. */
466 }
468 /** Assuming that the point up is outside the cutting plane, this routine
469 * searches upwards along edges trying to find an edge that intersects the
470 * cutting plane.
471 * \param[in] rsq the distance along this vector of the plane.
472 * \param[in,out] u the dot product of point up with the normal.
473 * \return True if the cutting plane was reached, false otherwise. */
474 inline bool voronoicell_base::search_upward(unsigned int &uw,int &lp,int &ls,int &us,double &l,double &u) {
478 // The test point is outside of the cutting space
483 }
487 }
490 //if(++count>=p) failsafe_find(lp,ls,us,l,u);
492 // Test all the neighbors of the current point
493 // and find the one which is closest to the
494 // plane
501 }
505 }
506 }
509 }
511 /** Checks whether a particular point lp is a definite maximum, searching
512 * through any possible minor non-convexities, for a better maximum.
513 * \param[in] (x,y,z) the normal vector to the plane. */
519 // Check to see whether point up is a well-defined maximum. Otherwise
520 // any neighboring vertices of up that are marginal need to be
521 // followed, to see if they lead to a better maximum.
526 }
529 // The point tp is marginal, so it will be necessary to do the
530 // flood-fill search. Mark the point tp and the point qp, and search
531 // any remaining neighbors of the point tp.
536 ts++;
544 }
545 ts++;
546 }
548 // Consider additional marginal points, starting with the original
549 // point qp
556 // Skip the point if it's already marked
560 // This point is a better maximum. Reset markers and
561 // return true.
572 }
574 // The point is marginal and therefore must also be
575 // considered
581 }
584 }
585 }
586 }
588 // Reset markers and return false
592 }
594 inline bool voronoicell_base::search_downward(unsigned int &lw,int &lp,int &ls,int &us,double &l,double &u) {
597 // The test point is outside of the cutting space
602 }
606 //if(++count>=p) failsafe_find(lp,ls,us,l,u);
608 // Test all the neighbors of the current point
609 // and find the one which is closest to the
610 // plane
617 }
619 }
622 }
629 // Check to see whether point up is a well-defined maximum. Otherwise
630 // any neighboring vertices of up that are marginal need to be
631 // followed, to see if they lead to a better maximum.
636 }
639 // The point tp is marginal, so it will be necessary to do the
640 // flood-fill search. Mark the point tp and the point qp, and search
641 // any remaining neighbors of the point tp.
646 ts++;
654 }
655 ts++;
656 }
658 // Consider additional marginal points, starting with the original
659 // point qp
666 // Skip the point if it's already marked
670 // This point is a better minimum. Reset markers and
671 // return true.
682 }
684 // The point is marginal and therefore must also be
685 // considered
691 }
694 }
695 }
696 }
698 // Reset markers and return false
702 }
704 /** Cuts the Voronoi cell by a particle whose center is at a separation of
705 * (x,y,z) from the cell center. The value of rsq should be initially set to
706 * \f$x^2+y^2+z^2\f$.
707 * \param[in] vc a reference to the specialized version of the calling class.
708 * \param[in] (x,y,z) the normal vector to the plane.
709 * \param[in] rsq the distance along this vector of the plane.
710 * \param[in] p_id the plane ID (for neighbor tracking only).
711 * \return False if the plane cut deleted the cell entirely, true otherwise. */
720 // Initialize the safe testing routine
734 }
736 // Set stack pointers
739 // Store initial number of vertices
750 // Skip if this is a new vertex
759 }
762 // This is a possible facet starting
763 // from a vertex on the cutting plane
767 // This is a new facet
771 }
772 }
773 xtra++;
774 }
776 // Reset back pointers on extra search stack
780 }
782 // Delete points: first, remove any duplicates
788 dsp++;
790 }
792 // Add the points in the auxiliary delete stack,
793 // and reset their back pointers
801 }
802 }
804 // Scan connections and add in extras
814 }
817 }
818 }
819 }
822 // Delete them from the array structure
833 }
837 // Vertex management
842 // Memory management
851 // Edge management
857 }
859 // Check for any vertices of zero order
862 // Collapse any order 2 vertices and exit
864 }
866 /** Creates a new facet.
867 * \return True if cell deleted, false otherwise. */
869 bool voronoicell_base::create_facet(vc_class &vc,int lp,int ls,double l,int us,double u,int p_id) {
875 // We're about to add the first point of the new facet. In either
876 // routine, we have to add a point, so first check there's space for
877 // it.
882 // We want to be strict about reaching the conclusion that the
883 // cell is entirely within the cutting plane. It's not enough
884 // to find a vertex that has edges which are all inside or on
885 // the plane. If the vertex has neighbors that are also on the
886 // plane, we should check those too.
889 // The search algorithm found a point which is on the cutting
890 // plane. We leave that point in place, and create a new one at
891 // the same location.
896 // Search for a collection of edges of the test vertex which
897 // are outside of the cutting space. Begin by testing the
898 // zeroth edge.
904 // The first edge is either inside the cutting space,
905 // or lies within the cutting plane. Test the edges
906 // sequentially until we find one that is outside.
909 i++;
911 // If we reached the last edge with no luck
912 // then all of the vertices are inside
913 // or on the plane, so the cell is completely
914 // deleted
921 // We found an edge outside the cutting space. Keep
922 // moving through these edges until we find one that's
923 // inside or on the plane.
928 j++;
929 }
931 // Compute the number of edges for the new vertex. In
932 // general it will be the number of outside edges
933 // found, plus two. But we need to recognize the
934 // special case when all but one edge is outside, and
935 // the remaining one is on the plane. For that case we
936 // have to reduce the edge count by one to prevent
937 // doubling up.
944 // Add memory for the new vertex if needed, and
945 // initialize
952 // Copy the edges of the original vertex into the new
953 // one. Delete the edges of the original vertex, and
954 // update the relational table.
966 }
970 // In this case, the zeroth edge is outside the cutting
971 // plane. Begin by searching backwards from the last
972 // edge until we find an edge which isn't outside.
977 i--;
979 // If i reaches zero, then we have a point in
980 // the plane all of whose edges are outside
981 // the cutting space, so we just exit
985 }
987 // Now search forwards from zero
992 j++;
995 }
997 // Compute the number of edges for the new vertex. In
998 // general it will be the number of outside edges
999 // found, plus two. But we need to recognize the
1000 // special case when all but one edge is outside, and
1001 // the remaining one is on the plane. For that case we
1002 // have to reduce the edge count by one to prevent
1003 // doubling up.
1009 }
1011 // Add memory to store the vertex if it doesn't exist
1012 // already
1017 // Copy the edges of the original vertex into the new
1018 // one. Delete the edges of the original vertex, and
1019 // update the relational table.
1034 }
1046 }
1048 }
1054 // Add this point to the auxiliary delete stack
1058 // Look at the edges on either side of the group that was
1059 // detected. We're going to commence facet computation by
1060 // moving along one of them. We are going to end up coming back
1061 // along the other one.
1071 // The search algorithm found an intersected edge between the
1072 // points lp and up. Create a new vertex between them which
1073 // lies on the cutting plane. Since u and l differ by at least
1074 // the tolerance, this division should never screw up.
1082 // This point will always have three edges. Connect one of them
1083 // to lp.
1099 // Set the direction to move in
1102 }
1104 // When the code reaches here, we have initialized the first point, and
1105 // we have a direction for moving it to construct the rest of the facet
1109 // We're currently tracing round an intersected facet. Keep
1110 // moving around it until we find a point or edge which
1111 // intersects the plane.
1117 // The point is still in the cutting space. Just add it
1118 // to the delete stack and keep moving.
1127 // The point is outside of the cutting space, so we've
1128 // found an intersected edge. Introduce a regular point
1129 // at the point of intersection. Connect it to the
1130 // point we just tested. Also connect it to the previous
1131 // new point in the facet we're constructing.
1160 // We've found a point which is on the cutting plane.
1161 // We're going to introduce a new point right here, but
1162 // first we need to figure out the number of edges it
1163 // has.
1166 // If the previous vertex detected a double edge, our
1167 // new vertex will have one less edge.
1173 // Start testing the edges of the current point until
1174 // we find one which isn't outside the cutting space
1176 k++;
1182 // Now we need to find out whether this marginal vertex
1183 // we are on has been visited before, because if that's
1184 // the case, we need to add vertices to the existing
1185 // new vertex, rather than creating a fresh one. We also
1186 // need to figure out whether we're in a case where we
1187 // might be creating a duplicate edge.
1191 // If we're heading into the final part of the
1192 // new facet, then we never worry about the
1193 // duplicate edge calculation.
1199 // This vertex was visited before, so
1200 // count those vertices to the ones we
1201 // already have.
1204 // The only time when we might make a
1205 // duplicate edge is if the point we're
1206 // going to move to next is also a
1207 // marginal point, so test for that
1208 // first.
1211 // Now see whether this marginal point
1212 // has been visited before.
1216 // Now see if the last edge of that other
1217 // marginal point actually ends up here.
1224 // That marginal point hasn't been visited
1225 // before, so we probably don't have to worry
1226 // about duplicate edges, except in the
1227 // case when that's the way into the end
1228 // of the facet, because that way always creates
1229 // an edge.
1234 }
1238 // The vertex hasn't been visited
1239 // before, but let's see if it's
1240 // marginal
1243 // If it is, we need to check
1244 // for the case that it's a
1245 // small branch, and that we're
1246 // heading right back to where
1247 // we came from
1254 }
1255 }
1257 // k now holds the number of edges of the new vertex
1258 // we are forming. Add memory for it if it doesn't exist
1259 // already.
1263 // Now create a new vertex with order k, or augment
1264 // the existing one
1267 // If we're augmenting a vertex but we don't
1268 // actually need any more edges, just skip this
1269 // routine to avoid memory confusion
1272 // Allocate memory and copy the edges
1273 // of the previous instance into it
1281 i++;
1282 }
1285 // Remove the previous instance with
1286 // fewer vertices from the memory
1287 // structure
1294 }
1300 // Allocate a new vertex of order k
1312 }
1315 // Unless the previous case was a double edge, connect
1316 // the first available edge of the new vertex to the
1317 // last one in the facet
1324 i++;
1325 }
1327 // Copy in the edges of the underlying vertex,
1328 // and do one less if this was a double edge
1339 i++;
1340 }
1346 // Update the double_edge flag, to pass it
1347 // to the next instance of this routine
1349 }
1350 }
1352 // Connect the final created vertex to the initial one
1358 }
1360 /** During the creation of a new facet in the plane routine, it is possible
1361 * that some order two vertices may arise. This routine removes them.
1362 * Suppose an order two vertex joins c and d. If there's a edge between
1363 * c and d already, then the order two vertex is just removed; otherwise,
1364 * the order two vertex is removed and c and d are joined together directly.
1365 * It is possible this process will create order two or order one vertices,
1366 * and the routine is continually run until all of them are removed.
1367 * \return False if the vertex removal was unsuccessful, indicative of the cell
1368 * reducing to zero volume and disappearing; true if the vertex removal
1369 * was successful. */
1376 // Pick a order 2 vertex and read in its edges
1380 #if VOROPP_VERBOSE >=1
1382 #endif
1384 }
1386 // Scan the edges of j to see if joins k
1389 }
1391 // If j doesn't already join k, join them together.
1392 // Otherwise delete the connection to the current
1393 // vertex from j and k.
1403 }
1405 // Compact the memory
1418 }
1420 // Collapse any order 1 vertices if they were created
1422 }
1424 }
1426 /** Order one vertices can potentially be created during the order two collapse
1427 * routine. This routine keeps removing them until there are none left.
1428 * \return False if the vertex removal was unsuccessful, indicative of the cell
1429 * having zero volume and disappearing; true if the vertex removal was
1430 * successful. */
1436 #if VOROPP_VERBOSE >=1
1438 #endif
1455 }
1456 }
1458 }
1460 /** This routine deletes the kth edge of vertex j and reorganizes the memory.
1461 * If the neighbor computation is enabled, we also have to supply an handedness
1462 * flag to decide whether to preserve the plane on the left or right of the
1463 * connection.
1464 * \return False if a zero order vertex was formed, indicative of the cell
1465 * disappearing; true if the vertex removal was successful. */
1470 #if VOROPP_VERBOSE >=1
1474 }
1475 #endif
1481 l++;
1482 }
1488 }
1495 l++;
1496 }
1507 }
1509 /** This routine is a fall-back, in case floating point errors caused the usual
1510 * search routine to fail. In the fall-back routine, we just test every edge to
1511 * find one straddling the plane. */
1515 /* qw=1;lw=0;
1516 for(qp=0;qp<p;qp++) {
1517 qw=m_test(qp,q);
1518 if(qw==1) {
1520 // The point is inside the cutting space. Now
1521 // see if we can find a neighbor which isn't.
1522 for(us=0;us<nu[qp];us++) {
1523 lp=ed[qp][us];
1524 if(lp<qp) {
1525 lw=m_test(lp,l);
1526 if(lw!=1) break;
1527 }
1528 }
1529 if(us<nu[qp]) {
1530 up=qp;
1531 if(lw==0) {
1532 complicated_setup=true;
1533 } else {
1534 complicated_setup=false;
1535 u=q;
1536 ls=ed[up][nu[up]+us];
1537 }
1538 break;
1539 }
1540 } else if(qw==-1) {
1542 // The point is outside the cutting space. See
1543 // if we can find a neighbor which isn't.
1544 for(ls=0;ls<nu[qp];ls++) {
1545 up=ed[qp][ls];
1546 if(up<qp) {
1547 uw=m_test(up,u);
1548 if(uw!=-1) break;
1549 }
1550 }
1551 if(ls<nu[qp]) {
1552 if(uw==0) {
1553 up=qp;
1554 complicated_setup=true;
1555 } else {
1556 complicated_setup=false;
1557 lp=qp;l=q;
1558 us=ed[lp][nu[lp]+ls];
1559 }
1560 break;
1561 }
1562 } else {
1564 // The point is in the plane, so we just
1565 // proceed with the complicated setup routine
1566 up=qp;
1567 complicated_setup=true;
1568 break;
1569 }
1570 }
1571 if(qp==p) return qw==-1?true:false;*/
1572 }
1574 /** Calculates the volume of the Voronoi cell, by decomposing the cell into
1575 * tetrahedra extending outward from the zeroth vertex, whose volumes are
1576 * evaluated using a scalar triple product.
1577 * \return A floating point number holding the calculated volume. */
1604 }
1605 }
1606 }
1607 }
1610 }
1612 /** Calculates the areas of each face of the Voronoi cell and prints the
1613 * results to an output stream.
1614 * \param[out] v the vector to store the results in. */
1641 }
1643 }
1644 }
1646 }
1649 /** Calculates the total surface area of the Voronoi cell.
1650 * \return The computed area. */
1675 }
1676 }
1677 }
1680 }
1683 /** Calculates the centroid of the Voronoi cell, by decomposing the cell into
1684 * tetrahedra extending outward from the zeroth vertex.
1685 * \param[out] (cx,cy,cz) references to floating point numbers in which to
1686 * pass back the centroid vector. */
1716 }
1717 }
1718 }
1719 }
1727 }
1729 /** Computes the maximum radius squared of a vertex from the center of the
1730 * cell. It can be used to determine when enough particles have been testing an
1731 * all planes that could cut the cell have been considered.
1732 * \return The maximum radius squared of a vertex.*/
1741 }
1743 }
1745 /** Calculates the total edge distance of the Voronoi cell.
1746 * \return A floating point number holding the calculated distance. */
1757 }
1758 }
1760 }
1762 /** Outputs the edges of the Voronoi cell in POV-Ray format to an open file
1763 * stream, displacing the cell by given vector.
1764 * \param[in] (x,y,z) a displacement vector to be added to the cell's position.
1765 * \param[in] fp a file handle to write to. */
1778 }
1779 }
1780 }
1781 }
1783 /** Outputs the edges of the Voronoi cell in gnuplot format to an output stream.
1784 * \param[in] (x,y,z) a displacement vector to be added to the cell's position.
1785 * \param[in] fp a file handle to write to. */
1800 }
1801 }
1803 }
1809 }
1811 }
1813 /** Outputs the Voronoi cell in the POV mesh2 format, described in section
1814 * 1.3.2.2 of the POV-Ray documentation. The mesh2 output consists of a list of
1815 * vertex vectors, followed by a list of triangular faces. The routine also
1816 * makes use of the optional inside_vector specification, which makes the mesh
1817 * object solid, so that the POV-Ray Constructive Solid Geometry (CSG) can be
1818 * applied.
1819 * \param[in] (x,y,z) a displacement vector to be added to the cell's position.
1820 * \param[in] fp a file handle to write to. */
1838 }
1839 }
1840 }
1843 }
1845 /** Several routines in the class that gather cell-based statistics internally
1846 * track their progress by flipping edges to negative so that they know what
1847 * parts of the cell have already been tested. This function resets them back
1848 * to positive. When it is called, it assumes that every edge in the routine
1849 * should have already been flipped to negative, and it bails out with an
1850 * internal error if it encounters a positive edge. */
1854 if(ed[i][j]>=0) voro_fatal_error("Edge reset routine found a previously untested edge",VOROPP_INTERNAL_ERROR);
1856 }
1857 }
1859 /** Checks to see if a given vertex is inside, outside or within the test
1860 * plane. If the point is far away from the test plane, the routine immediately
1861 * returns whether it is inside or outside. If the routine is close the the
1862 * plane and within the specified tolerance, then the special check_marginal()
1863 * routine is called.
1864 * \param[in] n the vertex to test.
1865 * \param[out] ans the result of the scalar product used in evaluating the
1866 * location of the point.
1867 * \return -1 if the point is inside the plane, 1 if the point is outside the
1868 * plane, or 0 if the point is within the plane. */
1874 }
1885 }
1888 /** Checks to see if a given vertex is inside, outside or within the test
1889 * plane. If the point is far away from the test plane, the routine immediately
1890 * returns whether it is inside or outside. If the routine is close the the
1891 * plane and within the specified tolerance, then the special check_marginal()
1892 * routine is called.
1893 * \param[in] n the vertex to test.
1894 * \param[out] ans the result of the scalar product used in evaluating the
1895 * location of the point.
1896 * \return -1 if the point is inside the plane, 1 if the point is outside the
1897 * plane, or 0 if the point is within the plane. */
1908 }
1910 }
1912 /** This routine calculates the unit normal vectors for every face.
1913 * \param[out] v the vector to store the results in. */
1920 }
1922 }
1924 /** This inline routine is called by normals(). It attempts to construct a
1925 * single normal vector that is associated with a particular face. It first
1926 * traces around the face, trying to find two vectors along the face edges
1927 * whose vector product is above the numerical tolerance. It then constructs
1928 * the normal vector using this product. If the face is too small, and none of
1929 * the vector products are large enough, the routine may return (0,0,0) as the
1930 * normal vector.
1931 * \param[in] v the vector to store the results in.
1932 * \param[in] i the initial vertex of the face to test.
1933 * \param[in] j the index of an edge of the vertex.
1934 * \param[in] k the neighboring vertex of i, set to ed[i][j]. */
1945 // Test to see if the length of this edge is above the tolerance
1954 // Construct the vector product of this edge with
1955 // the previous one
1961 // Test to see if this vector product of the
1962 // two edges is above the tolerance
1965 // Construct the normal vector and print it
1971 // Mark all of the remaining edges of this
1972 // face and exit
1976 }
1978 }
1979 }
1984 }
1991 }
1994 /** Returns the number of faces of a computed Voronoi cell.
1995 * \return The number of faces. */
2001 s++;
2011 }
2012 }
2015 }
2017 /** Returns a vector of the vertex orders.
2018 * \param[out] v the vector to store the results in. */
2022 }
2024 /** Outputs the vertex orders.
2025 * \param[out] fp the file handle to write to. */
2030 }
2031 }
2033 /** Returns a vector of the vertex vectors using the local coordinate system.
2034 * \param[out] v the vector to store the results in. */
2042 }
2043 }
2045 /** Outputs the vertex vectors using the local coordinate system.
2046 * \param[out] fp the file handle to write to. */
2050 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);
2051 }
2052 }
2055 /** Returns a vector of the vertex vectors in the global coordinate system.
2056 * \param[out] v the vector to store the results in.
2057 * \param[in] (x,y,z) the position vector of the particle in the global
2058 * coordinate system. */
2066 }
2067 }
2069 /** Outputs the vertex vectors using the global coordinate system.
2070 * \param[out] fp the file handle to write to.
2071 * \param[in] (x,y,z) the position vector of the particle in the global
2072 * coordinate system. */
2076 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);
2077 }
2078 }
2080 /** This routine returns the perimeters of each face.
2081 * \param[out] v the vector to store the results in. */
2106 }
2107 }
2109 }
2111 /** For each face, this routine outputs a bracketed sequence of numbers
2112 * containing a list of all the vertices that make up that face.
2113 * \param[out] v the vector to store the results in. */
2134 }
2135 }
2137 }
2139 /** Outputs a list of the number of edges in each face.
2140 * \param[out] v the vector to store the results in. */
2151 q++;
2158 }
2159 }
2161 }
2163 /** Computes the number of edges that each face has and outputs a frequency
2164 * table of the results.
2165 * \param[out] v the vector to store the results in. */
2176 q++;
2184 }
2185 }
2187 }
2189 /** This routine tests to see whether the cell intersects a plane by starting
2190 * from the guess point up. If up intersects, then it immediately returns true.
2191 * Otherwise, it calls the plane_intersects_track() routine.
2192 * \param[in] (x,y,z) the normal vector to the plane.
2193 * \param[in] rsq the distance along this vector of the plane.
2194 * \return False if the plane does not intersect the plane, true if it does. */
2199 }
2201 /** This routine tests to see if a cell intersects a plane. It first tests a
2202 * random sample of approximately sqrt(p)/4 points. If any of those are
2203 * intersect, then it immediately returns true. Otherwise, it takes the closest
2204 * point and passes that to plane_intersect_track() routine.
2205 * \param[in] (x,y,z) the normal vector to the plane.
2206 * \param[in] rsq the distance along this vector of the plane.
2207 * \return False if the plane does not intersect the plane, true if it does. */
2219 }
2221 }
2223 }
2225 }
2227 /* This routine tests to see if a cell intersects a plane, by tracing over the
2228 * cell from vertex to vertex, starting at up. It is meant to be called either
2229 * by plane_intersects() or plane_intersects_track(), when those routines
2230 * cannot immediately resolve the case.
2231 * \param[in] (x,y,z) the normal vector to the plane.
2232 * \param[in] rsq the distance along this vector of the plane.
2233 * \param[in] g the distance of up from the plane.
2234 * \return False if the plane does not intersect the plane, true if it does. */
2235 inline bool voronoicell_base::plane_intersects_track(double x,double y,double z,double rsq,double g) {
2239 /*
2240 int ls,us,lp;
2241 double l,u;
2242 unsigned int uw;
2244 // Initialize the safe testing routine
2245 px=x;py=y;pz=z;prsq=rsq;
2246 maskc+=4;
2247 if(maskc<4) reset_mask();
2249 return search_upward(uw,lp,ls,us,l,u);
2250 }*/
2251 /*
2252 int count=0,ls,us,tp;
2253 double t;
2254 // The test point is outside of the cutting space
2255 for(us=0;us<nu[up];us++) {
2256 tp=ed[up][us];
2257 t=x*pts[tp<<2]+y*pts[(tp<<2)+1]+z*pts[(tp<<2)+2];
2258 if(t>g) {
2259 ls=ed[up][nu[up]+us];
2260 up=tp;
2261 while (t<rsq) {
2262 if(++count>=p) {
2263 #if VOROPP_VERBOSE >=1
2264 fputs("Bailed out of convex calculation",stderr);
2265 #endif
2266 for(tp=0;tp<p;tp++) if(x*pts[tp<<2]+y*pts[(tp<<2)+1]+z*pts[(tp<<2)+2]>rsq) return true;
2267 return false;
2268 }
2270 // Test all the neighbors of the current point
2271 // and find the one which is closest to the
2272 // plane
2273 for(us=0;us<ls;us++) {
2274 tp=ed[up][us];double *pp=pts+(tp<<2);
2275 g=x*(*pp)+y*pp[1]+z*pp[2];
2276 if(g>t) break;
2277 }
2278 if(us==ls) {
2279 us++;
2280 while(us<nu[up]) {
2281 tp=ed[up][us];double *pp=pts+(tp<<2);
2282 g=x*(*pp)+y*pp[1]+z*pp[2];
2283 if(g>t) break;
2284 us++;
2285 }
2286 if(us==nu[up]) return false;
2287 }
2288 ls=ed[up][nu[up]+us];up=tp;t=g;
2289 }
2290 return true;
2291 }
2292 }
2293 return false;*/
2294 }
2296 /** Counts the number of edges of the Voronoi cell.
2297 * \return the number of edges. */
2302 }
2304 /** Outputs a custom string of information about the Voronoi cell. The string
2305 * of information follows a similar style as the C printf command, and detailed
2306 * information about its format is available at
2307 * http://math.lbl.gov/voro++/doc/custom.html.
2308 * \param[in] format the custom string to print.
2309 * \param[in] i the ID of the particle associated with this Voronoi cell.
2310 * \param[in] (x,y,z) the position of the particle associated with this Voronoi
2311 * cell.
2312 * \param[in] r a radius associated with the particle.
2313 * \param[in] fp the file handle to write to. */
2314 void voronoicell_base::output_custom(const char *format,int i,double x,double y,double z,double r,FILE *fp) {
2320 fmp++;
2323 // Particle-related output
2331 // Vertex-related output
2338 // Edge-related output
2343 // Face-related output
2363 // Volume-related output
2376 // End-of-string reached
2379 // The percent sign is not part of a
2380 // control sequence
2382 }
2384 fmp++;
2385 }
2387 }
2389 /** This initializes the class to be a rectangular box. It calls the base class
2390 * initialization routine to set up the edge and vertex information, and then
2391 * sets up the neighbor information, with initial faces being assigned ID
2392 * numbers from -1 to -6.
2393 * \param[in] (xmin,xmax) the minimum and maximum x coordinates.
2394 * \param[in] (ymin,ymax) the minimum and maximum y coordinates.
2395 * \param[in] (zmin,zmax) the minimum and maximum z coordinates. */
2396 void voronoicell_neighbor::init(double xmin,double xmax,double ymin,double ymax,double zmin,double zmax) {
2409 }
2411 /** This initializes the class to be an octahedron. It calls the base class
2412 * initialization routine to set up the edge and vertex information, and then
2413 * sets up the neighbor information, with the initial faces being assigned ID
2414 * numbers from -1 to -8.
2415 * \param[in] l The distance from the octahedron center to a vertex. Six
2416 * vertices are initialized at (-l,0,0), (l,0,0), (0,-l,0),
2417 * (0,l,0), (0,0,-l), and (0,0,l). */
2428 }
2430 /** This initializes the class to be a tetrahedron. It calls the base class
2431 * initialization routine to set up the edge and vertex information, and then
2432 * sets up the neighbor information, with the initial faces being assigned ID
2433 * numbers from -1 to -4.
2434 * \param (x0,y0,z0) a position vector for the first vertex.
2435 * \param (x1,y1,z1) a position vector for the second vertex.
2436 * \param (x2,y2,z2) a position vector for the third vertex.
2437 * \param (x3,y3,z3) a position vector for the fourth vertex. */
2438 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) {
2446 }
2448 /** This routine checks to make sure the neighbor information of each face is
2449 * consistent. */
2461 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);
2465 }
2466 }
2468 }
2470 /** The class constructor allocates memory for storing neighbor information. */
2478 }
2480 /** The class destructor frees the dynamically allocated memory for storing
2481 * neighbor information. */
2486 }
2488 /** Computes a vector list of neighbors. */
2504 }
2505 }
2507 }
2509 /** Prints the vertices, their edges, the relation table, and also notifies if
2510 * any memory errors are visible. */
2524 }
2525 }
2527 /** This prints out the neighbor information for vertex i. */
2535 }
2537 // Explicit instantiation
2541 template void voronoicell_base::check_memory_for_copy(voronoicell_neighbor&,voronoicell_base*);
2543 }