| blob | b831ff31a71bccf97882ebde5f66a37f15ad2078 |
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 <cstring>
19 /** Constructs a Voronoi cell and sets up the initial memory. */
33 }
39 }
40 }
42 /** The voronoicell destructor deallocates all the dynamic memory. */
50 }
52 /** Ensures that enough memory is allocated prior to carrying out a copy.
53 * \param[in] vc a reference to the specialized version of the calling class.
54 * \param[in] vb a pointered to the class to be copied. */
60 }
62 /** Copies the vertex and edge information from another class. The routine
63 * assumes that enough memory is available for the copy.
64 * \param[in] vb a pointer to the class to copy. */
72 }
75 }
77 /** Copies the information from another voronoicell class into this
78 * class, extending memory allocation if necessary.
79 * \param[in] c the class to copy. */
87 }
88 }
90 /** Copies the information from another voronoicell_neighbor class into this
91 * class, extending memory allocation if necessary.
92 * \param[in] c the class to copy. */
100 }
101 }
103 /** Translates the vertices of the Voronoi cell by a given vector.
104 * \param[in] (x,y,z) the coordinates of the vector. */
110 }
111 }
113 /** Increases the memory storage for a particular vertex order, by increasing
114 * the size of the of the corresponding mep array. If the arrays already exist,
115 * their size is doubled; if they don't exist, then new ones of size
116 * init_n_vertices are allocated. The routine also ensures that the pointers in
117 * the ed array are updated, by making use of the back pointers. For the cases
118 * where the back pointer has been temporarily overwritten in the marginal
119 * vertex code, the auxiliary delete stack is scanned to find out how to update
120 * the ed value. If the template has been instantiated with the neighbor
121 * tracking turned on, then the routine also reallocates the corresponding mne
122 * array.
123 * \param[in] i the order of the vertex memory to be increased. */
131 #if VOROPP_VERBOSE >=2
133 #endif
137 if(mem[i]>max_n_vertices) voro_fatal_error("Point memory allocation exceeded absolute maximum",VOROPP_MEMORY_ERROR);
138 #if VOROPP_VERBOSE >=2
140 #endif
156 }
157 }
159 #if VOROPP_VERBOSE >=3
161 #endif
162 }
165 }
169 }
170 }
172 /** Doubles the maximum number of vertices allowed, by reallocating the ed, nu,
173 * and pts arrays. If the allocation exceeds the absolute maximum set in
174 * max_vertices, then the routine exits with a fatal error. If the template has
175 * been instantiated with the neighbor tracking turned on, then the routine
176 * also reallocates the ne array. */
180 if(i>max_vertices) voro_fatal_error("Vertex memory allocation exceeded absolute maximum",VOROPP_MEMORY_ERROR);
181 #if VOROPP_VERBOSE >=2
183 #endif
196 }
198 /** Doubles the maximum allowed vertex order, by reallocating mem, mep, and mec
199 * arrays. If the allocation exceeds the absolute maximum set in
200 * max_vertex_order, then the routine causes a fatal error. If the template has
201 * been instantiated with the neighbor tracking turned on, then the routine
202 * also reallocates the mne array. */
206 if(i>max_vertex_order) voro_fatal_error("Vertex order memory allocation exceeded absolute maximum",VOROPP_MEMORY_ERROR);
207 #if VOROPP_VERBOSE >=2
209 #endif
221 }
223 /** Doubles the size allocation of the main delete stack. If the allocation
224 * exceeds the absolute maximum set in max_delete_size, then routine causes a
225 * fatal error. */
228 if(current_delete_size>max_delete_size) voro_fatal_error("Delete stack 1 memory allocation exceeded absolute maximum",VOROPP_MEMORY_ERROR);
229 #if VOROPP_VERBOSE >=2
231 #endif
236 }
238 /** Doubles the size allocation of the auxiliary delete stack. If the
239 * allocation exceeds the absolute maximum set in max_delete2_size, then the
240 * routine causes a fatal error. */
243 if(current_delete2_size>max_delete2_size) voro_fatal_error("Delete stack 2 memory allocation exceeded absolute maximum",VOROPP_MEMORY_ERROR);
244 #if VOROPP_VERBOSE >=2
246 #endif
251 }
253 /** Initializes a Voronoi cell as a rectangular box with the given dimensions.
254 * \param[in] (xmin,xmax) the minimum and maximum x coordinates.
255 * \param[in] (ymin,ymax) the minimum and maximum y coordinates.
256 * \param[in] (zmin,zmax) the minimum and maximum z coordinates. */
257 void voronoicell_base::init_base(double xmin,double xmax,double ymin,double ymax,double zmin,double zmax) {
280 }
282 /** Initializes a Voronoi cell as a regular octahedron.
283 * \param[in] l The distance from the octahedron center to a vertex. Six
284 * vertices are initialized at (-l,0,0), (l,0,0), (0,-l,0),
285 * (0,l,0), (0,0,-l), and (0,0,l). */
304 }
306 /** Initializes a Voronoi cell as a tetrahedron. It assumes that the normal to
307 * the face for the first three vertices points inside.
308 * \param (x0,y0,z0) a position vector for the first vertex.
309 * \param (x1,y1,z1) a position vector for the second vertex.
310 * \param (x2,y2,z2) a position vector for the third vertex.
311 * \param (x3,y3,z3) a position vector for the fourth vertex. */
312 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) {
326 }
328 /** Checks that the relational table of the Voronoi cell is accurate, and
329 * prints out any errors. This algorithm is O(p), so running it every time the
330 * plane routine is called will result in a significant slowdown. */
335 }
337 /** This routine checks for any two vertices that are connected by more than
338 * one edge. The plane algorithm is designed so that this should not happen, so
339 * any occurrences are most likely errors. Note that the routine is O(p), so
340 * running it every time the plane routine is called will result in a
341 * significant slowdown. */
346 }
348 /** Constructs the relational table if the edges have been specified. */
355 l++;
357 }
359 }
360 }
362 /** Starting from a point within the current cutting plane, this routine attempts
363 * to find an edge to a point outside the cutting plane. This prevents the plane
364 * routine from .
365 * \param[in] vc a reference to the specialized version of the calling class.
366 * \param[in,out] up */
379 }
380 }
382 }
384 /** Adds a point to the auxiliary delete stack if it is not already there.
385 * \param[in] vc a reference to the specialized version of the calling class.
386 * \param[in] lp the index of the point to add.
387 * \param[in,out] stackp2 a pointer to the end of the stack entries. */
393 }
395 /** Cuts the Voronoi cell by a particle whose center is at a separation of
396 * (x,y,z) from the cell center. The value of rsq should be initially set to
397 * \f$x^2+y^2+z^2\f$.
398 * \param[in] vc a reference to the specialized version of the calling class.
399 * \param[in] (x,y,z) the normal vector to the plane.
400 * \param[in] rsq the distance along this vector of the plane.
401 * \param[in] p_id the plane ID (for neighbor tracking only).
402 * \return False if the plane cut deleted the cell entirely, true otherwise. */
410 // Initialize the safe testing routine
413 // Test approximately sqrt(n)/4 points for their proximity to the plane
414 // and keep the one which is closest
417 // Starting from an initial guess, we now move from vertex to vertex,
418 // to try and find an edge which intersects the cutting plane,
419 // or a vertex which is on the plane
423 // The test point is inside the cutting plane.
429 us++;
434 }
444 }
446 us++;
451 us++;
452 }
455 }
456 }
458 }
460 // If the last point in the iteration is within the
461 // plane, we need to do the complicated setup
462 // routine. Otherwise, we use the regular iteration.
473 us++;
485 }
487 us++;
492 us++;
493 }
495 }
496 }
504 }
507 // Our original test point was on the plane, so we
508 // automatically head for the complicated setup
509 // routine
511 }
512 }
514 // This routine is a fall-back, in case floating point errors
515 // cause the usual search routine to fail. In the fall-back
516 // routine, we just test every edge to find one straddling
517 // the plane.
518 #if VOROPP_VERBOSE >=1
520 #endif
526 // The point is inside the cutting space. Now
527 // see if we can find a neighbor which isn't.
533 }
534 }
543 }
545 }
548 // The point is outside the cutting space. See
549 // if we can find a neighbor which isn't.
555 }
556 }
565 }
567 }
570 // The point is in the plane, so we just
571 // proceed with the complicated setup routine
575 }
576 }
578 }
580 // We're about to add the first point of the new facet. In either
581 // routine, we have to add a point, so first check there's space for
582 // it.
587 // We want to be strict about reaching the conclusion that the
588 // cell is entirely within the cutting plane. It's not enough
589 // to find a vertex that has edges which are all inside or on
590 // the plane. If the vertex has neighbors that are also on the
591 // plane, we should check those too.
594 // The search algorithm found a point which is on the cutting
595 // plane. We leave that point in place, and create a new one at
596 // the same location.
601 // Search for a collection of edges of the test vertex which
602 // are outside of the cutting space. Begin by testing the
603 // zeroth edge.
609 // The first edge is either inside the cutting space,
610 // or lies within the cutting plane. Test the edges
611 // sequentially until we find one that is outside.
614 i++;
616 // If we reached the last edge with no luck
617 // then all of the vertices are inside
618 // or on the plane, so the cell is completely
619 // deleted
626 // We found an edge outside the cutting space. Keep
627 // moving through these edges until we find one that's
628 // inside or on the plane.
633 j++;
634 }
636 // Compute the number of edges for the new vertex. In
637 // general it will be the number of outside edges
638 // found, plus two. But we need to recognize the
639 // special case when all but one edge is outside, and
640 // the remaining one is on the plane. For that case we
641 // have to reduce the edge count by one to prevent
642 // doubling up.
649 // Add memory for the new vertex if needed, and
650 // initialize
657 // Copy the edges of the original vertex into the new
658 // one. Delete the edges of the original vertex, and
659 // update the relational table.
671 }
675 // In this case, the zeroth edge is outside the cutting
676 // plane. Begin by searching backwards from the last
677 // edge until we find an edge which isn't outside.
682 i--;
684 // If i reaches zero, then we have a point in
685 // the plane all of whose edges are outside
686 // the cutting space, so we just exit
690 }
692 // Now search forwards from zero
697 j++;
700 }
702 // Compute the number of edges for the new vertex. In
703 // general it will be the number of outside edges
704 // found, plus two. But we need to recognize the
705 // special case when all but one edge is outside, and
706 // the remaining one is on the plane. For that case we
707 // have to reduce the edge count by one to prevent
708 // doubling up.
714 }
716 // Add memory to store the vertex if it doesn't exist
717 // already
722 // Copy the edges of the original vertex into the new
723 // one. Delete the edges of the original vertex, and
724 // update the relational table.
739 }
751 }
753 }
759 // Add this point to the auxiliary delete stack
763 // Look at the edges on either side of the group that was
764 // detected. We're going to commence facet computation by
765 // moving along one of them. We are going to end up coming back
766 // along the other one.
776 // The search algorithm found an intersected edge between the
777 // points lp and up. Create a new vertex between them which
778 // lies on the cutting plane. Since u and l differ by at least
779 // the tolerance, this division should never screw up.
787 // This point will always have three edges. Connect one of them
788 // to lp.
804 // Set the direction to move in
807 }
809 // When the code reaches here, we have initialized the first point, and
810 // we have a direction for moving it to construct the rest of the facet
814 // We're currently tracing round an intersected facet. Keep
815 // moving around it until we find a point or edge which
816 // intersects the plane.
822 // The point is still in the cutting space. Just add it
823 // to the delete stack and keep moving.
832 // The point is outside of the cutting space, so we've
833 // found an intersected edge. Introduce a regular point
834 // at the point of intersection. Connect it to the
835 // point we just tested. Also connect it to the previous
836 // new point in the facet we're constructing.
865 // We've found a point which is on the cutting plane.
866 // We're going to introduce a new point right here, but
867 // first we need to figure out the number of edges it
868 // has.
871 // If the previous vertex detected a double edge, our
872 // new vertex will have one less edge.
878 // Start testing the edges of the current point until
879 // we find one which isn't outside the cutting space
881 k++;
887 // Now we need to find out whether this marginal vertex
888 // we are on has been visited before, because if that's
889 // the case, we need to add vertices to the existing
890 // new vertex, rather than creating a fresh one. We also
891 // need to figure out whether we're in a case where we
892 // might be creating a duplicate edge.
896 // If we're heading into the final part of the
897 // new facet, then we never worry about the
898 // duplicate edge calculation.
904 // This vertex was visited before, so
905 // count those vertices to the ones we
906 // already have.
909 // The only time when we might make a
910 // duplicate edge is if the point we're
911 // going to move to next is also a
912 // marginal point, so test for that
913 // first.
916 // Now see whether this marginal point
917 // has been visited before.
921 // Now see if the last edge of that other
922 // marginal point actually ends up here.
929 // That marginal point hasn't been visited
930 // before, so we probably don't have to worry
931 // about duplicate edges, except in the
932 // case when that's the way into the end
933 // of the facet, because that way always creates
934 // an edge.
939 }
943 // The vertex hasn't been visited
944 // before, but let's see if it's
945 // marginal
948 // If it is, we need to check
949 // for the case that it's a
950 // small branch, and that we're
951 // heading right back to where
952 // we came from
959 }
960 }
962 // k now holds the number of edges of the new vertex
963 // we are forming. Add memory for it if it doesn't exist
964 // already.
968 // Now create a new vertex with order k, or augment
969 // the existing one
972 // If we're augmenting a vertex but we don't
973 // actually need any more edges, just skip this
974 // routine to avoid memory confusion
976 // Allocate memory and copy the edges
977 // of the previous instance into it
985 i++;
986 }
989 // Remove the previous instance with
990 // fewer vertices from the memory
991 // structure
998 }
1004 // Allocate a new vertex of order k
1016 }
1019 // Unless the previous case was a double edge, connect
1020 // the first available edge of the new vertex to the
1021 // last one in the facet
1028 i++;
1029 }
1031 // Copy in the edges of the underlying vertex,
1032 // and do one less if this was a double edge
1043 i++;
1044 }
1050 // Update the double_edge flag, to pass it
1051 // to the next instance of this routine
1053 }
1054 }
1056 // Connect the final created vertex to the initial one
1062 // Delete points: first, remove any duplicates
1068 dsp++;
1070 }
1072 // Add the points in the auxiliary delete stack,
1073 // and reset their back pointers
1081 }
1082 }
1084 // Scan connections and add in extras
1094 }
1097 }
1098 }
1099 }
1102 // Delete them from the array structure
1113 }
1117 // Vertex management
1122 // Memory management
1131 // Edge management
1137 }
1139 // Check for any vertices of zero order
1142 // Collapse any order 2 vertices and exit
1144 }
1146 /** During the creation of a new facet in the plane routine, it is possible
1147 * that some order two vertices may arise. This routine removes them.
1148 * Suppose an order two vertex joins c and d. If there's a edge between
1149 * c and d already, then the order two vertex is just removed; otherwise,
1150 * the order two vertex is removed and c and d are joined together directly.
1151 * It is possible this process will create order two or order one vertices,
1152 * and the routine is continually run until all of them are removed.
1153 * \return False if the vertex removal was unsuccessful, indicative of the cell
1154 * reducing to zero volume and disappearing; true if the vertex removal
1155 * was successful. */
1162 // Pick a order 2 vertex and read in its edges
1166 #if VOROPP_VERBOSE >=1
1168 #endif
1170 }
1172 // Scan the edges of j to see if joins k
1175 }
1177 // If j doesn't already join k, join them together.
1178 // Otherwise delete the connection to the current
1179 // vertex from j and k.
1189 }
1191 // Compact the memory
1204 }
1206 // Collapse any order 1 vertices if they were created
1208 }
1210 }
1212 /** Order one vertices can potentially be created during the order two collapse
1213 * routine. This routine keeps removing them until there are none left.
1214 * \return False if the vertex removal was unsuccessful, indicative of the cell
1215 * having zero volume and disappearing; true if the vertex removal was
1216 * successful. */
1222 #if VOROPP_VERBOSE >=1
1224 #endif
1241 }
1242 }
1244 }
1246 /** This routine deletes the kth edge of vertex j and reorganizes the memory.
1247 * If the neighbor computation is enabled, we also have to supply an handedness
1248 * flag to decide whether to preserve the plane on the left or right of the
1249 * connection.
1250 * \return False if a zero order vertex was formed, indicative of the cell
1251 * disappearing; true if the vertex removal was successful. */
1256 #if VOROPP_VERBOSE >=1
1260 }
1261 #endif
1267 l++;
1268 }
1274 }
1281 l++;
1282 }
1293 }
1295 /** Calculates the volume of the Voronoi cell, by decomposing the cell into
1296 * tetrahedra extending outward from the zeroth vertex, whose volumes are
1297 * evaluated using a scalar triple product.
1298 * \return A floating point number holding the calculated volume. */
1325 }
1326 }
1327 }
1328 }
1331 }
1333 /** Calculates the areas of each face of the Voronoi cell and prints the
1334 * results to an output stream.
1335 * \param[out] v the vector to store the results in. */
1362 }
1364 }
1365 }
1367 }
1370 /** Calculates the total surface area of the Voronoi cell.
1371 * \return The computed area. */
1396 }
1397 }
1398 }
1401 }
1404 /** Calculates the centroid of the Voronoi cell, by decomposing the cell into
1405 * tetrahedra extending outward from the zeroth vertex.
1406 * \param[out] (cx,cy,cz) references to floating point numbers in which to
1407 * pass back the centroid vector. */
1437 }
1438 }
1439 }
1440 }
1448 }
1450 /** Computes the maximum radius squared of a vertex from the center of the
1451 * cell. It can be used to determine when enough particles have been testing an
1452 * all planes that could cut the cell have been considered.
1453 * \return The maximum radius squared of a vertex.*/
1462 }
1464 }
1466 /** Calculates the total edge distance of the Voronoi cell.
1467 * \return A floating point number holding the calculated distance. */
1478 }
1479 }
1481 }
1483 /** Outputs the edges of the Voronoi cell in POV-Ray format to an open file
1484 * stream, displacing the cell by given vector.
1485 * \param[in] (x,y,z) a displacement vector to be added to the cell's position.
1486 * \param[in] fp a file handle to write to. */
1499 }
1500 }
1501 }
1502 }
1504 /** Outputs the edges of the Voronoi cell in gnuplot format to an output stream.
1505 * \param[in] (x,y,z) a displacement vector to be added to the cell's position.
1506 * \param[in] fp a file handle to write to. */
1521 }
1522 }
1524 }
1530 }
1532 }
1534 /** Outputs the Voronoi cell in the POV mesh2 format, described in section
1535 * 1.3.2.2 of the POV-Ray documentation. The mesh2 output consists of a list of
1536 * vertex vectors, followed by a list of triangular faces. The routine also
1537 * makes use of the optional inside_vector specification, which makes the mesh
1538 * object solid, so the the POV-Ray Constructive Solid Geometry (CSG) can be
1539 * applied.
1540 * \param[in] (x,y,z) a displacement vector to be added to the cell's position.
1541 * \param[in] fp a file handle to write to. */
1559 }
1560 }
1561 }
1564 }
1566 /** Several routines in the class that gather cell-based statistics internally
1567 * track their progress by flipping edges to negative so that they know what
1568 * parts of the cell have already been tested. This function resets them back
1569 * to positive. When it is called, it assumes that every edge in the routine
1570 * should have already been flipped to negative, and it bails out with an
1571 * internal error if it encounters a positive edge. */
1575 if(ed[i][j]>=0) voro_fatal_error("Edge reset routine found a previously untested edge",VOROPP_INTERNAL_ERROR);
1577 }
1578 }
1580 /** Checks to see if a given vertex is inside, outside or within the test
1581 * plane. If the point is far away from the test plane, the routine immediately
1582 * returns whether it is inside or outside. If the routine is close the the
1583 * plane and within the specified tolerance, then the special check_marginal()
1584 * routine is called.
1585 * \param[in] n the vertex to test.
1586 * \param[out] ans the result of the scalar product used in evaluating the
1587 * location of the point.
1588 * \return -1 if the point is inside the plane, 1 if the point is outside the
1589 * plane, or 0 if the point is within the plane. */
1599 }
1601 }
1603 /** Checks to see if a given vertex is inside, outside or within the test
1604 * plane, for the case when the point has been detected to be very close to the
1605 * plane. The routine ensures that the returned results are always consistent
1606 * with previous tests, by keeping a table of any marginal results. The routine
1607 * first sees if the vertex is in the table, and if it finds a previously
1608 * computed result it uses that. Otherwise, it computes a result for this
1609 * vertex and adds it the table.
1610 * \param[in] n the vertex to test.
1611 * \param[in] ans the result of the scalar product used in evaluating
1612 * the location of the point.
1613 * \return -1 if the point is inside the plane, 1 if the point is outside the
1614 * plane, or 0 if the point is within the plane. */
1621 voro_fatal_error("Marginal case buffer allocation exceeded absolute maximum",VOROPP_MEMORY_ERROR);
1622 #if VOROPP_VERBOSE >=2
1624 #endif
1629 }
1633 }
1635 /** For each face of the Voronoi cell, this routine prints the out the normal
1636 * vector of the face, and scales it to the distance from the cell center to
1637 * that plane.
1638 * \param[out] v the vector to store the results in. */
1645 }
1647 }
1649 /** This inline routine is called by normals(). It attempts to construct a
1650 * single normal vector that is associated with a particular face. It first
1651 * traces around the face, trying to find two vectors along the face edges
1652 * whose vector product is above the numerical tolerance. It then constructs
1653 * the normal vector using this product. If the face is too small, and none of
1654 * the vector products are large enough, the routine may return (0,0,0) as the
1655 * normal vector.
1656 * \param[in] v the vector to store the results in.
1657 * \param[in] i the initial vertex of the face to test.
1658 * \param[in] j the index of an edge of the vertex.
1659 * \param[in] k the neighboring vertex of i, set to ed[i][j]. */
1670 // Test to see if the length of this edge is above the tolerance
1679 // Construct the vector product of this edge with
1680 // the previous one
1686 // Test to see if this vector product of the
1687 // two edges is above the tolerance
1690 // Construct the normal vector and print it
1696 // Mark all of the remaining edges of this
1697 // face and exit
1701 }
1703 }
1704 }
1709 }
1716 }
1719 /** Returns the number of faces of a computed Voronoi cell.
1720 * \return The number of faces. */
1726 s++;
1736 }
1737 }
1740 }
1742 /** Returns a vector of the vertex orders.
1743 * \param[out] v the vector to store the results in. */
1747 }
1749 /** Outputs the vertex orders.
1750 * \param[out] fp the file handle to write to. */
1755 }
1756 }
1758 /** Returns a vector of the vertex vectors using the local coordinate system.
1759 * \param[out] v the vector to store the results in. */
1767 }
1768 }
1770 /** Outputs the vertex vectors using the local coordinate system.
1771 * \param[out] fp the file handle to write to. */
1775 for(double *ptsp=pts+3;ptsp<pts+3*p;ptsp+=3) fprintf(fp," (%g,%g,%g)",*ptsp*0.5,ptsp[1]*0.5,ptsp[2]*0.5);
1776 }
1777 }
1780 /** Returns a vector of the vertex vectors in the global coordinate system.
1781 * \param[out] v the vector to store the results in.
1782 * \param[in] (x,y,z) the position vector of the particle in the global
1783 * coordinate system. */
1791 }
1792 }
1794 /** Outputs the vertex vectors using the global coordinate system.
1795 * \param[out] fp the file handle to write to.
1796 * \param[in] (x,y,z) the position vector of the particle in the global
1797 * coordinate system. */
1801 for(double *ptsp=pts+3;ptsp<pts+3*p;ptsp+=3) fprintf(fp," (%g,%g,%g)",x+*ptsp*0.5,y+ptsp[1]*0.5,z+ptsp[2]*0.5);
1802 }
1803 }
1805 /** This routine returns the perimeters of each face.
1806 * \param[out] v the vector to store the results in. */
1831 }
1832 }
1834 }
1836 /** For each face, this routine outputs a bracketed sequence of numbers
1837 * containing a list of all the vertices that make up that face.
1838 * \param[out] v the vector to store the results in. */
1859 }
1860 }
1862 }
1864 /** Outputs a list of the number of edges in each face.
1865 * \param[out] v the vector to store the results in. */
1876 q++;
1883 }
1884 }
1886 }
1888 /** Computes the number of edges that each face has and outputs a frequency
1889 * table of the results.
1890 * \param[out] v the vector to store the results in. */
1901 q++;
1909 }
1910 }
1912 }
1914 /** This routine tests to see whether the cell intersects a plane by starting
1915 * from the guess point up. If up intersects, then it immediately returns true.
1916 * Otherwise, it calls the plane_intersects_track() routine.
1917 * \param[in] (x,y,z) the normal vector to the plane.
1918 * \param[in] rsq the distance along this vector of the plane.
1919 * \return False if the plane does not intersect the plane, true if it does. */
1924 }
1926 /** This routine tests to see if a cell intersects a plane. It first tests a
1927 * random sample of approximately sqrt(p)/4 points. If any of those are
1928 * intersect, then it immediately returns true. Otherwise, it takes the closest
1929 * point and passes that to plane_intersect_track() routine.
1930 * \param[in] (x,y,z) the normal vector to the plane.
1931 * \param[in] rsq the distance along this vector of the plane.
1932 * \return False if the plane does not intersect the plane, true if it does. */
1944 }
1946 }
1948 }
1950 }
1952 /* This routine tests to see if a cell intersects a plane, by tracing over the cell from
1953 * vertex to vertex, starting at up. It is meant to be called either by plane_intersects()
1954 * or plane_intersects_track(), when those routines cannot immediately resolve the case.
1955 * \param[in] (x,y,z) the normal vector to the plane.
1956 * \param[in] rsq the distance along this vector of the plane.
1957 * \param[in] g the distance of up from the plane.
1958 * \return False if the plane does not intersect the plane, true if it does. */
1959 inline bool voronoicell_base::plane_intersects_track(double x,double y,double z,double rsq,double g) {
1963 // The test point is outside of the cutting space
1972 #if VOROPP_VERBOSE >=1
1974 #endif
1977 }
1979 // Test all the neighbors of the current point
1980 // and find the one which is closest to the
1981 // plane
1986 }
1988 us++;
1993 us++;
1994 }
1996 }
1998 }
2000 }
2001 }
2003 }
2005 /** Counts the number of edges of the Voronoi cell.
2006 * \return the number of edges. */
2011 }
2013 /** Outputs a custom string of information about the Voronoi cell. The string
2014 * of information follows a similar style as the C printf command, and detailed
2015 * information about its format is available at
2016 * http://math.lbl.gov/voro++/doc/custom.html.
2017 * \param[in] format the custom string to print.
2018 * \param[in] i the ID of the particle associated with this Voronoi cell.
2019 * \param[in] (x,y,z) the position of the particle associated with this Voronoi
2020 * cell.
2021 * \param[in] r a radius associated with the particle.
2022 * \param[in] fp the file handle to write to. */
2023 void voronoicell_base::output_custom(const char *format,int i,double x,double y,double z,double r,FILE *fp) {
2029 fmp++;
2032 // Particle-related output
2040 // Vertex-related output
2047 // Edge-related output
2052 // Face-related output
2072 // Volume-related output
2085 // End-of-string reached
2088 // The percent sign is not part of a
2089 // control sequence
2091 }
2093 fmp++;
2094 }
2096 }
2098 /** This initializes the class to be a rectangular box. It calls the base class
2099 * initialization routine to set up the edge and vertex information, and then
2100 * sets up the neighbor information, with initial faces being assigned ID
2101 * numbers from -1 to -6.
2102 * \param[in] (xmin,xmax) the minimum and maximum x coordinates.
2103 * \param[in] (ymin,ymax) the minimum and maximum y coordinates.
2104 * \param[in] (zmin,zmax) the minimum and maximum z coordinates. */
2105 void voronoicell_neighbor::init(double xmin,double xmax,double ymin,double ymax,double zmin,double zmax) {
2118 }
2120 /** This initializes the class to be an octahedron. It calls the base class
2121 * initialization routine to set up the edge and vertex information, and then
2122 * sets up the neighbor information, with the initial faces being assigned ID
2123 * numbers from -1 to -8.
2124 * \param[in] l The distance from the octahedron center to a vertex. Six
2125 * vertices are initialized at (-l,0,0), (l,0,0), (0,-l,0),
2126 * (0,l,0), (0,0,-l), and (0,0,l). */
2137 }
2139 /** This initializes the class to be a tetrahedron. It calls the base class
2140 * initialization routine to set up the edge and vertex information, and then
2141 * sets up the neighbor information, with the initial faces being assigned ID
2142 * numbers from -1 to -4.
2143 * \param (x0,y0,z0) a position vector for the first vertex.
2144 * \param (x1,y1,z1) a position vector for the second vertex.
2145 * \param (x2,y2,z2) a position vector for the third vertex.
2146 * \param (x3,y3,z3) a position vector for the fourth vertex. */
2147 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) {
2155 }
2157 /** This routine checks to make sure the neighbor information of each face is
2158 * consistent. */
2170 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);
2174 }
2175 }
2177 }
2179 /** The class constructor allocates memory for storing neighbor information. */
2187 }
2189 /** The class destructor frees the dynamically allocated memory for storing
2190 * neighbor information. */
2195 }
2197 /** Computes a vector list of neighbors. */
2213 }
2214 }
2216 }
2218 /** Prints the vertices, their edges, the relation table, and also notifies if
2219 * any memory errors are visible. */
2233 }
2234 }
2236 /** This prints out the neighbor information for vertex i. */
2244 }
2246 // Explicit instantiation
2250 template void voronoicell_base::check_memory_for_copy(voronoicell_neighbor&,voronoicell_base*);
2252 }