| blob | 9aff884fcc950f98f1d162ba1be7cbf4cdae2839 |
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. */
34 }
40 }
41 }
43 /** The voronoicell destructor deallocates all the dynamic memory. */
51 }
53 /** Ensures that enough memory is allocated prior to carrying out a copy.
54 * \param[in] vc a reference to the specialized version of the calling class.
55 * \param[in] vb a pointered to the class to be copied. */
61 }
63 /** Copies the vertex and edge information from another class. The routine
64 * assumes that enough memory is available for the copy.
65 * \param[in] vb a pointer to the class to copy. */
73 }
77 }
79 /** Copies the information from another voronoicell class into this
80 * class, extending memory allocation if necessary.
81 * \param[in] c the class to copy. */
89 }
90 }
92 /** Copies the information from another voronoicell_neighbor class into this
93 * class, extending memory allocation if necessary.
94 * \param[in] c the class to copy. */
102 }
103 }
105 /** Translates the vertices of the Voronoi cell by a given vector.
106 * \param[in] (x,y,z) the coordinates of the vector. */
113 }
114 }
116 /** Increases the memory storage for a particular vertex order, by increasing
117 * the size of the of the corresponding mep array. If the arrays already exist,
118 * their size is doubled; if they don't exist, then new ones of size
119 * init_n_vertices are allocated. The routine also ensures that the pointers in
120 * the ed array are updated, by making use of the back pointers. For the cases
121 * where the back pointer has been temporarily overwritten in the marginal
122 * vertex code, the auxiliary delete stack is scanned to find out how to update
123 * the ed value. If the template has been instantiated with the neighbor
124 * tracking turned on, then the routine also reallocates the corresponding mne
125 * array.
126 * \param[in] i the order of the vertex memory to be increased. */
134 #if VOROPP_VERBOSE >=2
136 #endif
140 if(mem[i]>max_n_vertices) voro_fatal_error("Point memory allocation exceeded absolute maximum",VOROPP_MEMORY_ERROR);
141 #if VOROPP_VERBOSE >=2
143 #endif
159 }
160 }
162 #if VOROPP_VERBOSE >=3
164 #endif
165 }
168 }
172 }
173 }
175 /** Doubles the maximum number of vertices allowed, by reallocating the ed, nu,
176 * and pts arrays. If the allocation exceeds the absolute maximum set in
177 * max_vertices, then the routine exits with a fatal error. If the template has
178 * been instantiated with the neighbor tracking turned on, then the routine
179 * also reallocates the ne array. */
183 if(i>max_vertices) voro_fatal_error("Vertex memory allocation exceeded absolute maximum",VOROPP_MEMORY_ERROR);
184 #if VOROPP_VERBOSE >=2
186 #endif
201 }
203 /** Doubles the maximum allowed vertex order, by reallocating mem, mep, and mec
204 * arrays. If the allocation exceeds the absolute maximum set in
205 * max_vertex_order, then the routine causes a fatal error. If the template has
206 * been instantiated with the neighbor tracking turned on, then the routine
207 * also reallocates the mne array. */
211 if(i>max_vertex_order) voro_fatal_error("Vertex order memory allocation exceeded absolute maximum",VOROPP_MEMORY_ERROR);
212 #if VOROPP_VERBOSE >=2
214 #endif
226 }
228 /** Doubles the size allocation of the main delete stack. If the allocation
229 * exceeds the absolute maximum set in max_delete_size, then routine causes a
230 * fatal error. */
233 if(current_delete_size>max_delete_size) voro_fatal_error("Delete stack 1 memory allocation exceeded absolute maximum",VOROPP_MEMORY_ERROR);
234 #if VOROPP_VERBOSE >=2
236 #endif
241 }
243 /** Doubles the size allocation of the auxiliary delete stack. If the
244 * allocation exceeds the absolute maximum set in max_delete2_size, then the
245 * routine causes a fatal error. */
248 if(current_delete2_size>max_delete2_size) voro_fatal_error("Delete stack 2 memory allocation exceeded absolute maximum",VOROPP_MEMORY_ERROR);
249 #if VOROPP_VERBOSE >=2
251 #endif
256 }
258 /** Initializes a Voronoi cell as a rectangular box with the given dimensions.
259 * \param[in] (xmin,xmax) the minimum and maximum x coordinates.
260 * \param[in] (ymin,ymax) the minimum and maximum y coordinates.
261 * \param[in] (zmin,zmax) the minimum and maximum z coordinates. */
262 void voronoicell_base::init_base(double xmin,double xmax,double ymin,double ymax,double zmin,double zmax) {
286 }
288 /** Initializes a Voronoi cell as a regular octahedron.
289 * \param[in] l The distance from the octahedron center to a vertex. Six
290 * vertices are initialized at (-l,0,0), (l,0,0), (0,-l,0),
291 * (0,l,0), (0,0,-l), and (0,0,l). */
311 }
313 /** Initializes a Voronoi cell as a tetrahedron. It assumes that the normal to
314 * the face for the first three vertices points inside.
315 * \param (x0,y0,z0) a position vector for the first vertex.
316 * \param (x1,y1,z1) a position vector for the second vertex.
317 * \param (x2,y2,z2) a position vector for the third vertex.
318 * \param (x3,y3,z3) a position vector for the fourth vertex. */
319 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) {
333 }
335 /** Checks that the relational table of the Voronoi cell is accurate, and
336 * prints out any errors. This algorithm is O(p), so running it every time the
337 * plane routine is called will result in a significant slowdown. */
342 }
344 /** This routine checks for any two vertices that are connected by more than
345 * one edge. The plane algorithm is designed so that this should not happen, so
346 * any occurrences are most likely errors. Note that the routine is O(p), so
347 * running it every time the plane routine is called will result in a
348 * significant slowdown. */
353 }
355 /** Constructs the relational table if the edges have been specified. */
362 l++;
364 }
366 }
367 }
369 /** Cuts the Voronoi cell by a particle whose center is at a separation of
370 * (x,y,z) from the cell center. The value of rsq should be initially set to
371 * \f$x^2+y^2+z^2\f$.
372 * \param[in] vc a reference to the specialized version of the calling class.
373 * \param[in] (x,y,z) the normal vector to the plane.
374 * \param[in] rsq the distance along this vector of the plane.
375 * \param[in] p_id the plane ID (for neighbor tracking only).
376 * \return False if the plane cut deleted the cell entirely, true otherwise. */
384 // Initialize the safe testing routine
387 // Test approximately sqrt(n)/4 points for their proximity to the plane
388 // and keep the one which is closest
391 // Starting from an initial guess, we now move from vertex to vertex,
392 // to try and find an edge which intersects the cutting plane,
393 // or a vertex which is on the plane
397 // The test point is inside the cutting plane.
403 us++;
408 }
418 }
420 us++;
425 us++;
426 }
429 }
430 }
432 }
434 // If the last point in the iteration is within the
435 // plane, we need to do the complicated setup
436 // routine. Otherwise, we use the regular iteration.
447 us++;
459 }
461 us++;
466 us++;
467 }
469 }
470 }
478 }
481 // Our original test point was on the plane, so we
482 // automatically head for the complicated setup
483 // routine
485 }
486 }
488 // This routine is a fall-back, in case floating point errors
489 // cause the usual search routine to fail. In the fall-back
490 // routine, we just test every edge to find one straddling
491 // the plane.
492 #if VOROPP_VERBOSE >=1
494 #endif
500 // The point is inside the cutting space. Now
501 // see if we can find a neighbor which isn't.
507 }
508 }
517 }
519 }
522 // The point is outside the cutting space. See
523 // if we can find a neighbor which isn't.
529 }
530 }
539 }
541 }
544 // The point is in the plane, so we just
545 // proceed with the complicated setup routine
549 }
550 }
552 }
556 // We're about to add the first point of the new facet. In either
557 // routine, we have to add a point, so first check there's space for
558 // it.
563 // The search algorithm found a point which is on the cutting
564 // plane. We leave that point in place, and create a new one at
565 // the same location.
570 // Search for a collection of edges of the test vertex which
571 // are outside of the cutting space. Begin by testing the
572 // zeroth edge.
578 // The first edge is either inside the cutting space,
579 // or lies within the cutting plane. Test the edges
580 // sequentially until we find one that is outside.
583 i++;
585 // If we reached the last edge with no luck
586 // then all of the vertices are inside
587 // or on the plane, so the cell is completely
588 // deleted
595 // We found an edge outside the cutting space. Keep
596 // moving through these edges until we find one that's
597 // inside or on the plane.
602 j++;
603 }
605 // Compute the number of edges for the new vertex. In
606 // general it will be the number of outside edges
607 // found, plus two. But we need to recognize the
608 // special case when all but one edge is outside, and
609 // the remaining one is on the plane. For that case we
610 // have to reduce the edge count by one to prevent
611 // doubling up.
618 // Add memory for the new vertex if needed, and
619 // initialize
626 // Copy the edges of the original vertex into the new
627 // one. Delete the edges of the original vertex, and
628 // update the relational table.
640 }
644 // In this case, the zeroth edge is outside the cutting
645 // plane. Begin by searching backwards from the last
646 // edge until we find an edge which isn't outside.
651 i--;
653 // If i reaches zero, then we have a point in
654 // the plane all of whose edges are outside
655 // the cutting space, so we just exit
659 }
661 // Now search forwards from zero
666 j++;
669 }
671 // Compute the number of edges for the new vertex. In
672 // general it will be the number of outside edges
673 // found, plus two. But we need to recognize the
674 // special case when all but one edge is outside, and
675 // the remaining one is on the plane. For that case we
676 // have to reduce the edge count by one to prevent
677 // doubling up.
683 }
685 // Add memory to store the vertex if it doesn't exist
686 // already
691 // Copy the edges of the original vertex into the new
692 // one. Delete the edges of the original vertex, and
693 // update the relational table.
708 }
720 }
722 }
728 // Add this point to the auxiliary delete stack
732 // Look at the edges on either side of the group that was
733 // detected. We're going to commence facet computation by
734 // moving along one of them. We are going to end up coming back
735 // along the other one.
745 // The search algorithm found an intersected edge between the
746 // points lp and up. Create a new vertex between them which
747 // lies on the cutting plane. Since u and l differ by at least
748 // the tolerance, this division should never screw up.
756 // This point will always have three edges. Connect one of them
757 // to lp.
773 // Set the direction to move in
776 }
778 // When the code reaches here, we have initialized the first point, and
779 // we have a direction for moving it to construct the rest of the facet
783 // We're currently tracing round an intersected facet. Keep
784 // moving around it until we find a point or edge which
785 // intersects the plane.
791 // The point is still in the cutting space. Just add it
792 // to the delete stack and keep moving.
801 // The point is outside of the cutting space, so we've
802 // found an intersected edge. Introduce a regular point
803 // at the point of intersection. Connect it to the
804 // point we just tested. Also connect it to the previous
805 // new point in the facet we're constructing.
834 // We've found a point which is on the cutting plane.
835 // We're going to introduce a new point right here, but
836 // first we need to figure out the number of edges it
837 // has.
840 // If the previous vertex detected a double edge, our
841 // new vertex will have one less edge.
847 // Start testing the edges of the current point until
848 // we find one which isn't outside the cutting space
850 k++;
856 // Now we need to find out whether this marginal vertex
857 // we are on has been visited before, because if that's
858 // the case, we need to add vertices to the existing
859 // new vertex, rather than creating a fresh one. We also
860 // need to figure out whether we're in a case where we
861 // might be creating a duplicate edge.
865 // If we're heading into the final part of the
866 // new facet, then we never worry about the
867 // duplicate edge calculation.
873 // This vertex was visited before, so
874 // count those vertices to the ones we
875 // already have.
878 // The only time when we might make a
879 // duplicate edge is if the point we're
880 // going to move to next is also a
881 // marginal point, so test for that
882 // first.
885 // Now see whether this marginal point
886 // has been visited before.
890 // Now see if the last edge of that other
891 // marginal point actually ends up here.
898 // That marginal point hasn't been visited
899 // before, so we probably don't have to worry
900 // about duplicate edges, except in the
901 // case when that's the way into the end
902 // of the facet, because that way always creates
903 // an edge.
908 }
912 // The vertex hasn't been visited
913 // before, but let's see if it's
914 // marginal
917 // If it is, we need to check
918 // for the case that it's a
919 // small branch, and that we're
920 // heading right back to where
921 // we came from
928 }
929 }
931 // k now holds the number of edges of the new vertex
932 // we are forming. Add memory for it if it doesn't exist
933 // already.
937 // Now create a new vertex with order k, or augment
938 // the existing one
941 // If we're augmenting a vertex but we don't
942 // actually need any more edges, just skip this
943 // routine to avoid memory confusion
945 // Allocate memory and copy the edges
946 // of the previous instance into it
954 i++;
955 }
958 // Remove the previous instance with
959 // fewer vertices from the memory
960 // structure
967 }
973 // Allocate a new vertex of order k
985 }
988 // Unless the previous case was a double edge, connect
989 // the first available edge of the new vertex to the
990 // last one in the facet
997 i++;
998 }
1000 // Copy in the edges of the underlying vertex,
1001 // and do one less if this was a double edge
1012 i++;
1013 }
1019 // Update the double_edge flag, to pass it
1020 // to the next instance of this routine
1022 }
1023 }
1025 // Connect the final created vertex to the initial one
1031 // Delete points: first, remove any duplicates
1037 dsp++;
1039 }
1041 // Add the points in the auxiliary delete stack,
1042 // and reset their back pointers
1050 }
1051 }
1053 // Scan connections and add in extras
1063 }
1066 }
1067 }
1068 }
1071 // Delete them from the array structure
1082 }
1086 // Vertex management
1091 // Memory management
1100 // Edge management
1106 }
1108 // Check for any vertices of zero order
1111 // Collapse any order 2 vertices and exit
1113 }
1115 /** During the creation of a new facet in the plane routine, it is possible
1116 * that some order two vertices may arise. This routine removes them.
1117 * Suppose an order two vertex joins c and d. If there's a edge between
1118 * c and d already, then the order two vertex is just removed; otherwise,
1119 * the order two vertex is removed and c and d are joined together directly.
1120 * It is possible this process will create order two or order one vertices,
1121 * and the routine is continually run until all of them are removed.
1122 * \return False if the vertex removal was unsuccessful, indicative of the cell
1123 * reducing to zero volume and disappearing; true if the vertex removal
1124 * was successful. */
1131 // Pick a order 2 vertex and read in its edges
1135 #if VOROPP_VERBOSE >=1
1137 #endif
1139 }
1141 // Scan the edges of j to see if joins k
1144 }
1146 // If j doesn't already join k, join them together.
1147 // Otherwise delete the connection to the current
1148 // vertex from j and k.
1158 }
1160 // Compact the memory
1173 }
1175 // Collapse any order 1 vertices if they were created
1177 }
1179 }
1181 /** Order one vertices can potentially be created during the order two collapse
1182 * routine. This routine keeps removing them until there are none left.
1183 * \return False if the vertex removal was unsuccessful, indicative of the cell
1184 * having zero volume and disappearing; true if the vertex removal was
1185 * successful. */
1191 #if VOROPP_VERBOSE >=1
1193 #endif
1210 }
1211 }
1213 }
1215 /** This routine deletes the kth edge of vertex j and reorganizes the memory.
1216 * If the neighbor computation is enabled, we also have to supply an handedness
1217 * flag to decide whether to preserve the plane on the left or right of the
1218 * connection.
1219 * \return False if a zero order vertex was formed, indicative of the cell
1220 * disappearing; true if the vertex removal was successful. */
1225 #if VOROPP_VERBOSE >=1
1229 }
1230 #endif
1236 l++;
1237 }
1243 }
1250 l++;
1251 }
1262 }
1264 /** Calculates the volume of the Voronoi cell, by decomposing the cell into
1265 * tetrahedra extending outward from the zeroth vertex, whose volumes are
1266 * evaluated using a scalar triple product.
1267 * \return A floating point number holding the calculated volume. */
1295 }
1296 }
1297 }
1298 }
1301 }
1303 /** Calculates the areas of each face of the Voronoi cell and prints the
1304 * results to an output stream.
1305 * \param[out] v the vector to store the results in. */
1333 }
1335 }
1336 }
1338 }
1341 /** Calculates the total surface area of the Voronoi cell.
1342 * \return The computed area. */
1368 }
1369 }
1370 }
1373 }
1376 /** Calculates the centroid of the Voronoi cell, by decomposing the cell into
1377 * tetrahedra extending outward from the zeroth vertex.
1378 * \param[out] (cx,cy,cz) references to floating point numbers in which to
1379 * pass back the centroid vector. */
1410 }
1411 }
1412 }
1413 }
1421 }
1423 /** Computes the maximum radius squared of a vertex from the center of the
1424 * cell. It can be used to determine when enough particles have been testing an
1425 * all planes that could cut the cell have been considered.
1426 * \return The maximum radius squared of a vertex.*/
1436 }
1438 }
1440 /** Calculates the total edge distance of the Voronoi cell.
1441 * \return A floating point number holding the calculated distance. */
1453 }
1454 }
1456 }
1458 /** Outputs the edges of the Voronoi cell in POV-Ray format to an open file
1459 * stream, displacing the cell by given vector.
1460 * \param[in] (x,y,z) a displacement vector to be added to the cell's position.
1461 * \param[in] fp a file handle to write to. */
1475 }
1476 }
1477 }
1478 }
1480 /** Outputs the edges of the Voronoi cell in gnuplot format to an output stream.
1481 * \param[in] (x,y,z) a displacement vector to be added to the cell's position.
1482 * \param[in] fp a file handle to write to. */
1498 }
1499 }
1501 }
1507 }
1509 }
1511 /** Outputs the Voronoi cell in the POV mesh2 format, described in section
1512 * 1.3.2.2 of the POV-Ray documentation. The mesh2 output consists of a list of
1513 * vertex vectors, followed by a list of triangular faces. The routine also
1514 * makes use of the optional inside_vector specification, which makes the mesh
1515 * object solid, so the the POV-Ray Constructive Solid Geometry (CSG) can be
1516 * applied.
1517 * \param[in] (x,y,z) a displacement vector to be added to the cell's position.
1518 * \param[in] fp a file handle to write to. */
1537 }
1538 }
1539 }
1542 }
1544 /** Several routines in the class that gather cell-based statistics internally
1545 * track their progress by flipping edges to negative so that they know what
1546 * parts of the cell have already been tested. This function resets them back
1547 * to positive. When it is called, it assumes that every edge in the routine
1548 * should have already been flipped to negative, and it bails out with an
1549 * internal error if it encounters a positive edge. */
1553 if(ed[i][j]>=0) voro_fatal_error("Edge reset routine found a previously untested edge",VOROPP_INTERNAL_ERROR);
1555 }
1556 }
1558 /** Checks to see if a given vertex is inside, outside or within the test
1559 * plane. If the point is far away from the test plane, the routine immediately
1560 * returns whether it is inside or outside. If the routine is close the the
1561 * plane and within the specified tolerance, then the special check_marginal()
1562 * routine is called.
1563 * \param[in] n the vertex to test.
1564 * \param[out] ans the result of the scalar product used in evaluating the
1565 * location of the point.
1566 * \return -1 if the point is inside the plane, 1 if the point is outside the
1567 * plane, or 0 if the point is within the plane. */
1578 }
1580 /** For each face of the Voronoi cell, this routine prints the out the normal
1581 * vector of the face, and scales it to the distance from the cell center to
1582 * that plane.
1583 * \param[out] v the vector to store the results in. */
1591 }
1593 }
1595 /** This inline routine is called by normals(). It attempts to construct a
1596 * single normal vector that is associated with a particular face. It first
1597 * traces around the face, trying to find two vectors along the face edges
1598 * whose vector product is above the numerical tolerance. It then constructs
1599 * the normal vector using this product. If the face is too small, and none of
1600 * the vector products are large enough, the routine may return (0,0,0) as the
1601 * normal vector.
1602 * \param[in] v the vector to store the results in.
1603 * \param[in] i the initial vertex of the face to test.
1604 * \param[in] j the index of an edge of the vertex.
1605 * \param[in] k the neighboring vertex of i, set to ed[i][j]. */
1616 // Test to see if the length of this edge is above the tolerance
1625 // Construct the vector product of this edge with
1626 // the previous one
1632 // Test to see if this vector product of the
1633 // two edges is above the tolerance
1636 // Construct the normal vector and print it
1642 // Mark all of the remaining edges of this
1643 // face and exit
1647 }
1649 }
1650 }
1655 }
1662 }
1665 /** Returns the number of faces of a computed Voronoi cell.
1666 * \return The number of faces. */
1672 s++;
1682 }
1683 }
1686 }
1688 /** Returns a vector of the vertex orders.
1689 * \param[out] v the vector to store the results in. */
1693 }
1695 /** Outputs the vertex orders.
1696 * \param[out] fp the file handle to write to. */
1701 }
1702 }
1704 /** Returns a vector of the vertex vectors using the local coordinate system.
1705 * \param[out] v the vector to store the results in. */
1714 }
1715 }
1717 /** Outputs the vertex vectors using the local coordinate system.
1718 * \param[out] fp the file handle to write to. */
1723 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);
1724 }
1725 }
1728 /** Returns a vector of the vertex vectors in the global coordinate system.
1729 * \param[out] v the vector to store the results in.
1730 * \param[in] (x,y,z) the position vector of the particle in the global
1731 * coordinate system. */
1740 }
1741 }
1743 /** Outputs the vertex vectors using the global coordinate system.
1744 * \param[out] fp the file handle to write to.
1745 * \param[in] (x,y,z) the position vector of the particle in the global
1746 * coordinate system. */
1751 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);
1752 }
1753 }
1755 /** This routine returns the perimeters of each face.
1756 * \param[out] v the vector to store the results in. */
1782 }
1783 }
1785 }
1787 /** For each face, this routine outputs a bracketed sequence of numbers
1788 * containing a list of all the vertices that make up that face.
1789 * \param[out] v the vector to store the results in. */
1810 }
1811 }
1813 }
1815 /** Outputs a list of the number of edges in each face.
1816 * \param[out] v the vector to store the results in. */
1827 q++;
1834 }
1835 }
1837 }
1839 /** Computes the number of edges that each face has and outputs a frequency
1840 * table of the results.
1841 * \param[out] v the vector to store the results in. */
1852 q++;
1860 }
1861 }
1863 }
1865 /** This routine tests to see whether the cell intersects a plane by starting
1866 * from the guess point up. If up intersects, then it immediately returns true.
1867 * Otherwise, it calls the plane_intersects_track() routine.
1868 * \param[in] (x,y,z) the normal vector to the plane.
1869 * \param[in] rsq the distance along this vector of the plane.
1870 * \return False if the plane does not intersect the plane, true if it does. */
1877 }
1879 /** This routine tests to see if a cell intersects a plane. It first tests a
1880 * random sample of approximately sqrt(p)/4 points. If any of those are
1881 * intersect, then it immediately returns true. Otherwise, it takes the closest
1882 * point and passes that to plane_intersect_track() routine.
1883 * \param[in] (x,y,z) the normal vector to the plane.
1884 * \param[in] rsq the distance along this vector of the plane.
1885 * \return False if the plane does not intersect the plane, true if it does. */
1898 }
1900 }
1902 }
1904 }
1906 /* This routine tests to see if a cell intersects a plane, by tracing over the cell from
1907 * vertex to vertex, starting at up. It is meant to be called either by plane_intersects()
1908 * or plane_intersects_track(), when those routines cannot immediately resolve the case.
1909 * \param[in] (x,y,z) the normal vector to the plane.
1910 * \param[in] rsq the distance along this vector of the plane.
1911 * \param[in] g the distance of up from the plane.
1912 * \return False if the plane does not intersect the plane, true if it does. */
1913 inline bool voronoicell_base::plane_intersects_track(mpq_class x,mpq_class y,mpq_class z,mpq_class rsq,mpq_class g) {
1915 mpq_class t;
1917 // The test point is outside of the cutting space
1926 #if VOROPP_VERBOSE >=1
1928 #endif
1931 }
1933 // Test all the neighbors of the current point
1934 // and find the one which is closest to the
1935 // plane
1940 }
1942 us++;
1947 us++;
1948 }
1950 }
1952 }
1954 }
1955 }
1957 }
1959 /** Counts the number of edges of the Voronoi cell.
1960 * \return the number of edges. */
1965 }
1967 /** Outputs a custom string of information about the Voronoi cell. The string
1968 * of information follows a similar style as the C printf command, and detailed
1969 * information about its format is available at
1970 * http://math.lbl.gov/voro++/doc/custom.html.
1971 * \param[in] format the custom string to print.
1972 * \param[in] i the ID of the particle associated with this Voronoi cell.
1973 * \param[in] (x,y,z) the position of the particle associated with this Voronoi
1974 * cell.
1975 * \param[in] r a radius associated with the particle.
1976 * \param[in] fp the file handle to write to. */
1977 void voronoicell_base::output_custom(const char *format,int i,double x,double y,double z,double r,FILE *fp) {
1983 fmp++;
1986 // Particle-related output
1994 // Vertex-related output
2001 // Edge-related output
2006 // Face-related output
2026 // Volume-related output
2039 // End-of-string reached
2042 // The percent sign is not part of a
2043 // control sequence
2045 }
2047 fmp++;
2048 }
2050 }
2052 /** This initializes the class to be a rectangular box. It calls the base class
2053 * initialization routine to set up the edge and vertex information, and then
2054 * sets up the neighbor information, with initial faces being assigned ID
2055 * numbers from -1 to -6.
2056 * \param[in] (xmin,xmax) the minimum and maximum x coordinates.
2057 * \param[in] (ymin,ymax) the minimum and maximum y coordinates.
2058 * \param[in] (zmin,zmax) the minimum and maximum z coordinates. */
2059 void voronoicell_neighbor::init(double xmin,double xmax,double ymin,double ymax,double zmin,double zmax) {
2072 }
2074 /** This initializes the class to be an octahedron. It calls the base class
2075 * initialization routine to set up the edge and vertex information, and then
2076 * sets up the neighbor information, with the initial faces being assigned ID
2077 * numbers from -1 to -8.
2078 * \param[in] l The distance from the octahedron center to a vertex. Six
2079 * vertices are initialized at (-l,0,0), (l,0,0), (0,-l,0),
2080 * (0,l,0), (0,0,-l), and (0,0,l). */
2091 }
2093 /** This initializes the class to be a tetrahedron. It calls the base class
2094 * initialization routine to set up the edge and vertex information, and then
2095 * sets up the neighbor information, with the initial faces being assigned ID
2096 * numbers from -1 to -4.
2097 * \param (x0,y0,z0) a position vector for the first vertex.
2098 * \param (x1,y1,z1) a position vector for the second vertex.
2099 * \param (x2,y2,z2) a position vector for the third vertex.
2100 * \param (x3,y3,z3) a position vector for the fourth vertex. */
2101 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) {
2109 }
2111 /** This routine checks to make sure the neighbor information of each face is
2112 * consistent. */
2124 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);
2128 }
2129 }
2131 }
2133 /** The class constructor allocates memory for storing neighbor information. */
2141 }
2143 /** The class destructor frees the dynamically allocated memory for storing
2144 * neighbor information. */
2149 }
2151 /** Computes a vector list of neighbors. */
2167 }
2168 }
2170 }
2172 /** Prints the vertices, their edges, the relation table, and also notifies if
2173 * any memory errors are visible. */
2188 }
2189 }
2191 /** This prints out the neighbor information for vertex i. */
2199 }
2203 }
2205 // Explicit instantiation
2209 template void voronoicell_base::check_memory_for_copy(voronoicell_neighbor&,voronoicell_base*);
2211 }