| blob | 484cbdd53454c976275ff8224b4855ac9abef241 |
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 }
288 /** Initializes a Voronoi cell as a rectangular box with the given dimensions.
289 * \param[in] (xmin,xmax) the minimum and maximum x coordinates.
290 * \param[in] (ymin,ymax) the minimum and maximum y coordinates.
291 * \param[in] (zmin,zmax) the minimum and maximum z coordinates. */
292 void voronoicell_base::init_base(double xmin,double xmax,double ymin,double ymax,double zmin,double zmax) {
315 }
317 /** Initializes an L-shaped Voronoi cell of a fixed size for testing the
318 * 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 contributions to the Minkowski functionals for this Voronoi cell.
1613 * \param[in] r the radius to consider.
1614 * \param[out] ar the area functional.
1615 * \param[out] vo the volume functional. */
1630 }
1631 }
1632 }
1636 }
1638 inline void voronoicell_base::minkowski_contrib(int i,int k,int m,double r,double &ar,double &vo) {
1649 // Compute second orthonormal vector
1656 }
1660 // Compute third orthonormal vector
1674 }
1676 void voronoicell_base::minkowski_edge(double x0,double r1,double s1,double r2,double s2,double r,double &ar,double &vo) {
1684 }
1686 void voronoicell_base::minkowski_formula(double x0,double y0,double z0,double r,double &ar,double &vo) {
1692 double xs=x0*x0,ys=y0*y0,zs=z0*z0,res=xs+ys,rvs=res+zs,theta=atan(z0/y0),rs=r*r,rc=rs*r,temp,voc,arc;
1709 arc=x0*r*temp-0.5*temp2*y0*r/((rs-xs)*temp5)+x0*y0/6.*r/temp5+rs*0.5*temp3+rs*rs/3.*2*xs*ys/(res*(rs-xs)*sqrt((rs-xs)*(rs-xs)-(x2s-y2s-xs)*(x2s-y2s-xs)))-rs*0.5*temp4;
1713 }
1716 }
1718 /** Calculates the areas of each face of the Voronoi cell and prints the
1719 * results to an output stream.
1720 * \param[out] v the vector to store the results in. */
1747 }
1749 }
1750 }
1752 }
1754 /** Calculates the total surface area of the Voronoi cell.
1755 * \return The computed area. */
1780 }
1781 }
1782 }
1785 }
1787 /** Calculates the centroid of the Voronoi cell, by decomposing the cell into
1788 * tetrahedra extending outward from the zeroth vertex.
1789 * \param[out] (cx,cy,cz) references to floating point numbers in which to
1790 * pass back the centroid vector. */
1820 }
1821 }
1822 }
1823 }
1831 }
1833 /** Computes the maximum radius squared of a vertex from the center of the
1834 * cell. It can be used to determine when enough particles have been testing an
1835 * all planes that could cut the cell have been considered.
1836 * \return The maximum radius squared of a vertex.*/
1845 }
1847 }
1849 /** Calculates the total edge distance of the Voronoi cell.
1850 * \return A floating point number holding the calculated distance. */
1861 }
1862 }
1864 }
1866 /** Outputs the edges of the Voronoi cell in POV-Ray format to an open file
1867 * stream, displacing the cell by given vector.
1868 * \param[in] (x,y,z) a displacement vector to be added to the cell's position.
1869 * \param[in] fp a file handle to write to. */
1882 }
1883 }
1884 }
1885 }
1887 /** Outputs the edges of the Voronoi cell in gnuplot format to an output stream.
1888 * \param[in] (x,y,z) a displacement vector to be added to the cell's position.
1889 * \param[in] fp a file handle to write to. */
1904 }
1905 }
1907 }
1913 }
1915 }
1917 /** Outputs the Voronoi cell in the POV mesh2 format, described in section
1918 * 1.3.2.2 of the POV-Ray documentation. The mesh2 output consists of a list of
1919 * vertex vectors, followed by a list of triangular faces. The routine also
1920 * makes use of the optional inside_vector specification, which makes the mesh
1921 * object solid, so that the POV-Ray Constructive Solid Geometry (CSG) can be
1922 * applied.
1923 * \param[in] (x,y,z) a displacement vector to be added to the cell's position.
1924 * \param[in] fp a file handle to write to. */
1942 }
1943 }
1944 }
1947 }
1949 /** Several routines in the class that gather cell-based statistics internally
1950 * track their progress by flipping edges to negative so that they know what
1951 * parts of the cell have already been tested. This function resets them back
1952 * to positive. When it is called, it assumes that every edge in the routine
1953 * should have already been flipped to negative, and it bails out with an
1954 * internal error if it encounters a positive edge. */
1958 if(ed[i][j]>=0) voro_fatal_error("Edge reset routine found a previously untested edge",VOROPP_INTERNAL_ERROR);
1960 }
1961 }
1963 /** Checks to see if a given vertex is inside, outside or within the test
1964 * plane. If the point is far away from the test plane, the routine immediately
1965 * returns whether it is inside or outside. If the routine is close the the
1966 * plane and within the specified tolerance, then the special check_marginal()
1967 * routine is called.
1968 * \param[in] n the vertex to test.
1969 * \param[out] ans the result of the scalar product used in evaluating the
1970 * location of the point.
1971 * \return -1 if the point is inside the plane, 1 if the point is outside the
1972 * plane, or 0 if the point is within the plane. */
1978 }
1989 }
1991 /** Checks to see if a given vertex is inside, outside or within the test
1992 * plane. If the point is far away from the test plane, the routine immediately
1993 * returns whether it is inside or outside. If the routine is close the the
1994 * plane and within the specified tolerance, then the special check_marginal()
1995 * routine is called.
1996 * \param[in] n the vertex to test.
1997 * \param[out] ans the result of the scalar product used in evaluating the
1998 * location of the point.
1999 * \return -1 if the point is inside the plane, 1 if the point is outside the
2000 * plane, or 0 if the point is within the plane. */
2011 }
2013 }
2015 /** This routine calculates the unit normal vectors for every face.
2016 * \param[out] v the vector to store the results in. */
2023 }
2025 }
2027 /** This inline routine is called by normals(). It attempts to construct a
2028 * single normal vector that is associated with a particular face. It first
2029 * traces around the face, trying to find two vectors along the face edges
2030 * whose vector product is above the numerical tolerance. It then constructs
2031 * the normal vector using this product. If the face is too small, and none of
2032 * the vector products are large enough, the routine may return (0,0,0) as the
2033 * normal vector.
2034 * \param[in] v the vector to store the results in.
2035 * \param[in] i the initial vertex of the face to test.
2036 * \param[in] j the index of an edge of the vertex.
2037 * \param[in] k the neighboring vertex of i, set to ed[i][j]. */
2048 // Test to see if the length of this edge is above the tolerance
2057 // Construct the vector product of this edge with
2058 // the previous one
2064 // Test to see if this vector product of the
2065 // two edges is above the tolerance
2068 // Construct the normal vector and print it
2074 // Mark all of the remaining edges of this
2075 // face and exit
2079 }
2081 }
2082 }
2087 }
2094 }
2096 /** Returns the number of faces of a computed Voronoi cell.
2097 * \return The number of faces. */
2103 s++;
2113 }
2114 }
2117 }
2119 /** Returns a vector of the vertex orders.
2120 * \param[out] v the vector to store the results in. */
2124 }
2126 /** Outputs the vertex orders.
2127 * \param[out] fp the file handle to write to. */
2132 }
2133 }
2135 /** Returns a vector of the vertex vectors using the local coordinate system.
2136 * \param[out] v the vector to store the results in. */
2144 }
2145 }
2147 /** Outputs the vertex vectors using the local coordinate system.
2148 * \param[out] fp the file handle to write to. */
2152 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);
2153 }
2154 }
2156 /** Returns a vector of the vertex vectors in the global coordinate system.
2157 * \param[out] v the vector to store the results in.
2158 * \param[in] (x,y,z) the position vector of the particle in the global
2159 * coordinate system. */
2167 }
2168 }
2170 /** Outputs the vertex vectors using the global coordinate system.
2171 * \param[out] fp the file handle to write to.
2172 * \param[in] (x,y,z) the position vector of the particle in the global
2173 * coordinate system. */
2177 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);
2178 }
2179 }
2181 /** This routine returns the perimeters of each face.
2182 * \param[out] v the vector to store the results in. */
2207 }
2208 }
2210 }
2212 /** For each face, this routine outputs a bracketed sequence of numbers
2213 * containing a list of all the vertices that make up that face.
2214 * \param[out] v the vector to store the results in. */
2235 }
2236 }
2238 }
2240 /** Outputs a list of the number of edges in each face.
2241 * \param[out] v the vector to store the results in. */
2252 q++;
2259 }
2260 }
2262 }
2264 /** Computes the number of edges that each face has and outputs a frequency
2265 * table of the results.
2266 * \param[out] v the vector to store the results in. */
2277 q++;
2285 }
2286 }
2288 }
2290 /** This routine tests to see whether the cell intersects a plane by starting
2291 * from the guess point up. If up intersects, then it immediately returns true.
2292 * Otherwise, it calls the plane_intersects_track() routine.
2293 * \param[in] (x,y,z) the normal vector to the plane.
2294 * \param[in] rsq the distance along this vector of the plane.
2295 * \return False if the plane does not intersect the plane, true if it does. */
2300 }
2302 /** This routine tests to see if a cell intersects a plane. It first tests a
2303 * random sample of approximately sqrt(p)/4 points. If any of those are
2304 * intersect, then it immediately returns true. Otherwise, it takes the closest
2305 * point and passes that to plane_intersect_track() routine.
2306 * \param[in] (x,y,z) the normal vector to the plane.
2307 * \param[in] rsq the distance along this vector of the plane.
2308 * \return False if the plane does not intersect the plane, true if it does. */
2320 }
2322 }
2324 }
2326 }
2328 /* This routine tests to see if a cell intersects a plane, by tracing over the
2329 * cell from vertex to vertex, starting at up. It is meant to be called either
2330 * by plane_intersects() or plane_intersects_track(), when those routines
2331 * cannot immediately resolve the case.
2332 * \param[in] (x,y,z) the normal vector to the plane.
2333 * \param[in] rsq the distance along this vector of the plane.
2334 * \param[in] g the distance of up from the plane.
2335 * \return False if the plane does not intersect the plane, true if it does. */
2336 inline bool voronoicell_base::plane_intersects_track(double x,double y,double z,double rsq,double g) {
2340 /*
2341 int ls,us,lp;
2342 double l,u;
2343 unsigned int uw;
2345 // Initialize the safe testing routine
2346 px=x;py=y;pz=z;prsq=rsq;
2347 maskc+=4;
2348 if(maskc<4) reset_mask();
2350 return search_upward(uw,lp,ls,us,l,u);
2351 }*/
2352 /*
2353 int count=0,ls,us,tp;
2354 double t;
2355 // The test point is outside of the cutting space
2356 for(us=0;us<nu[up];us++) {
2357 tp=ed[up][us];
2358 t=x*pts[tp<<2]+y*pts[(tp<<2)+1]+z*pts[(tp<<2)+2];
2359 if(t>g) {
2360 ls=ed[up][nu[up]+us];
2361 up=tp;
2362 while (t<rsq) {
2363 if(++count>=p) {
2364 #if VOROPP_VERBOSE >=1
2365 fputs("Bailed out of convex calculation",stderr);
2366 #endif
2367 for(tp=0;tp<p;tp++) if(x*pts[tp<<2]+y*pts[(tp<<2)+1]+z*pts[(tp<<2)+2]>rsq) return true;
2368 return false;
2369 }
2371 // Test all the neighbors of the current point
2372 // and find the one which is closest to the
2373 // plane
2374 for(us=0;us<ls;us++) {
2375 tp=ed[up][us];double *pp=pts+(tp<<2);
2376 g=x*(*pp)+y*pp[1]+z*pp[2];
2377 if(g>t) break;
2378 }
2379 if(us==ls) {
2380 us++;
2381 while(us<nu[up]) {
2382 tp=ed[up][us];double *pp=pts+(tp<<2);
2383 g=x*(*pp)+y*pp[1]+z*pp[2];
2384 if(g>t) break;
2385 us++;
2386 }
2387 if(us==nu[up]) return false;
2388 }
2389 ls=ed[up][nu[up]+us];up=tp;t=g;
2390 }
2391 return true;
2392 }
2393 }
2394 return false;*/
2395 }
2397 /** Counts the number of edges of the Voronoi cell.
2398 * \return the number of edges. */
2403 }
2405 /** Outputs a custom string of information about the Voronoi cell. The string
2406 * of information follows a similar style as the C printf command, and detailed
2407 * information about its format is available at
2408 * http://math.lbl.gov/voro++/doc/custom.html.
2409 * \param[in] format the custom string to print.
2410 * \param[in] i the ID of the particle associated with this Voronoi cell.
2411 * \param[in] (x,y,z) the position of the particle associated with this Voronoi
2412 * cell.
2413 * \param[in] r a radius associated with the particle.
2414 * \param[in] fp the file handle to write to. */
2415 void voronoicell_base::output_custom(const char *format,int i,double x,double y,double z,double r,FILE *fp) {
2421 fmp++;
2424 // Particle-related output
2432 // Vertex-related output
2439 // Edge-related output
2444 // Face-related output
2464 // Volume-related output
2477 // End-of-string reached
2480 // The percent sign is not part of a
2481 // control sequence
2483 }
2485 fmp++;
2486 }
2488 }
2490 /** This initializes the class to be a rectangular box. It calls the base class
2491 * initialization routine to set up the edge and vertex information, and then
2492 * sets up the neighbor information, with initial faces being assigned ID
2493 * numbers from -1 to -6.
2494 * \param[in] (xmin,xmax) the minimum and maximum x coordinates.
2495 * \param[in] (ymin,ymax) the minimum and maximum y coordinates.
2496 * \param[in] (zmin,zmax) the minimum and maximum z coordinates. */
2497 void voronoicell_neighbor::init(double xmin,double xmax,double ymin,double ymax,double zmin,double zmax) {
2510 }
2512 /** This initializes the class to be an octahedron. It calls the base class
2513 * initialization routine to set up the edge and vertex information, and then
2514 * sets up the neighbor information, with the initial faces being assigned ID
2515 * numbers from -1 to -8.
2516 * \param[in] l The distance from the octahedron center to a vertex. Six
2517 * vertices are initialized at (-l,0,0), (l,0,0), (0,-l,0),
2518 * (0,l,0), (0,0,-l), and (0,0,l). */
2529 }
2531 /** This initializes the class to be a tetrahedron. It calls the base class
2532 * initialization routine to set up the edge and vertex information, and then
2533 * sets up the neighbor information, with the initial faces being assigned ID
2534 * numbers from -1 to -4.
2535 * \param (x0,y0,z0) a position vector for the first vertex.
2536 * \param (x1,y1,z1) a position vector for the second vertex.
2537 * \param (x2,y2,z2) a position vector for the third vertex.
2538 * \param (x3,y3,z3) a position vector for the fourth vertex. */
2539 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) {
2547 }
2549 /** This routine checks to make sure the neighbor information of each face is
2550 * consistent. */
2562 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);
2566 }
2567 }
2569 }
2571 /** The class constructor allocates memory for storing neighbor information. */
2579 }
2581 /** The class destructor frees the dynamically allocated memory for storing
2582 * neighbor information. */
2587 }
2589 /** Computes a vector list of neighbors. */
2605 }
2606 }
2608 }
2610 /** Prints the vertices, their edges, the relation table, and also notifies if
2611 * any memory errors are visible. */
2625 }
2626 }
2628 /** This prints out the neighbor information for vertex i. */
2636 }
2638 // Explicit instantiation
2642 template void voronoicell_base::check_memory_for_copy(voronoicell_neighbor&,voronoicell_base*);
2644 }