255#define FILENAMESIZE 2048
260#define INPUTLINESIZE 1024
266#define TRIPERBLOCK 4092
267#define SUBSEGPERBLOCK 508
268#define VERTEXPERBLOCK 4092
269#define VIRUSPERBLOCK 1020
271#define BADSUBSEGPERBLOCK 252
273#define BADTRIPERBLOCK 4092
275#define FLIPSTACKERPERBLOCK 252
277#define SPLAYNODEPERBLOCK 508
284#define SEGMENTVERTEX 1
286#define DEADVERTEX -32768
287#define UNDEADVERTEX -32767
295#define SAMPLEFACTOR 11
305#define PI 3.141592653589793238462643383279502884197169399375105820974944592308
309#define SQUAREROOTTWO 1.4142135623730950488016887242096980785696718753769480732
313#define ONETHIRD 0.333333333333333333333333333333333333333333333333333333333333
913#define decode(ptr, otri) \
914 (otri).orient = (int) ((uintptr_t) (ptr) & (uintptr_t) 3l); \
915 (otri).tri = (triangle *) \
916 ((uintptr_t) (ptr) ^ (uintptr_t) (otri).orient)
922#define encode(otri) \
923 (triangle) ((uintptr_t) (otri).tri | (uintptr_t) (otri).orient)
933#define sym(otri1, otri2) \
934 ptr = (otri1).tri[(otri1).orient]; \
937#define symself(otri) \
938 ptr = (otri).tri[(otri).orient]; \
943#define lnext(otri1, otri2) \
944 (otri2).tri = (otri1).tri; \
945 (otri2).orient = plus1mod3[(otri1).orient]
947#define lnextself(otri) \
948 (otri).orient = plus1mod3[(otri).orient]
952#define lprev(otri1, otri2) \
953 (otri2).tri = (otri1).tri; \
954 (otri2).orient = minus1mod3[(otri1).orient]
956#define lprevself(otri) \
957 (otri).orient = minus1mod3[(otri).orient]
963#define onext(otri1, otri2) \
964 lprev(otri1, otri2); \
967#define onextself(otri) \
975#define oprev(otri1, otri2) \
979#define oprevself(otri) \
987#define dnext(otri1, otri2) \
991#define dnextself(otri) \
999#define dprev(otri1, otri2) \
1000 lnext(otri1, otri2); \
1003#define dprevself(otri) \
1011#define rnext(otri1, otri2) \
1012 sym(otri1, otri2); \
1016#define rnextself(otri) \
1025#define rprev(otri1, otri2) \
1026 sym(otri1, otri2); \
1030#define rprevself(otri) \
1038#define org(otri, vertexptr) \
1039 vertexptr = (vertex) (otri).tri[plus1mod3[(otri).orient] + 3]
1041#define dest(otri, vertexptr) \
1042 vertexptr = (vertex) (otri).tri[minus1mod3[(otri).orient] + 3]
1044#define apex(otri, vertexptr) \
1045 vertexptr = (vertex) (otri).tri[(otri).orient + 3]
1047#define setorg(otri, vertexptr) \
1048 (otri).tri[plus1mod3[(otri).orient] + 3] = (triangle) vertexptr
1050#define setdest(otri, vertexptr) \
1051 (otri).tri[minus1mod3[(otri).orient] + 3] = (triangle) vertexptr
1053#define setapex(otri, vertexptr) \
1054 (otri).tri[(otri).orient + 3] = (triangle) vertexptr
1058#define bond(otri1, otri2) \
1059 (otri1).tri[(otri1).orient] = encode(otri2); \
1060 (otri2).tri[(otri2).orient] = encode(otri1)
1067#define dissolve(otri) \
1068 (otri).tri[(otri).orient] = (triangle) m->dummytri
1072#define otricopy(otri1, otri2) \
1073 (otri2).tri = (otri1).tri; \
1074 (otri2).orient = (otri1).orient
1078#define otriequal(otri1, otri2) \
1079 (((otri1).tri == (otri2).tri) && \
1080 ((otri1).orient == (otri2).orient))
1085#define infect(otri) \
1086 (otri).tri[6] = (triangle) \
1087 ((uintptr_t) (otri).tri[6] | (uintptr_t) 2l)
1089#define uninfect(otri) \
1090 (otri).tri[6] = (triangle) \
1091 ((uintptr_t) (otri).tri[6] & ~ (uintptr_t) 2l)
1095#define infected(otri) \
1096 (((uintptr_t) (otri).tri[6] & (uintptr_t) 2l) != 0l)
1100#define elemattribute(otri, attnum) \
1101 ((REAL *) (otri).tri)[m->elemattribindex + (attnum)]
1103#define setelemattribute(otri, attnum, value) \
1104 ((REAL *) (otri).tri)[m->elemattribindex + (attnum)] = value
1108#define areabound(otri) ((REAL *) (otri).tri)[m->areaboundindex]
1110#define setareabound(otri, value) \
1111 ((REAL *) (otri).tri)[m->areaboundindex] = value
1118#define deadtri(tria) ((tria)[1] == (triangle) NULL)
1120#define killtri(tria) \
1121 (tria)[1] = (triangle) NULL; \
1122 (tria)[3] = (triangle) NULL
1133#define sdecode(sptr, osub) \
1134 (osub).ssorient = (int) ((uintptr_t) (sptr) & (uintptr_t) 1l); \
1135 (osub).ss = (subseg *) \
1136 ((uintptr_t) (sptr) & ~ (uintptr_t) 3l)
1142#define sencode(osub) \
1143 (subseg) ((uintptr_t) (osub).ss | (uintptr_t) (osub).ssorient)
1147#define ssym(osub1, osub2) \
1148 (osub2).ss = (osub1).ss; \
1149 (osub2).ssorient = 1 - (osub1).ssorient
1151#define ssymself(osub) \
1152 (osub).ssorient = 1 - (osub).ssorient
1157#define spivot(osub1, osub2) \
1158 sptr = (osub1).ss[(osub1).ssorient]; \
1159 sdecode(sptr, osub2)
1161#define spivotself(osub) \
1162 sptr = (osub).ss[(osub).ssorient]; \
1168#define snext(osub1, osub2) \
1169 sptr = (osub1).ss[1 - (osub1).ssorient]; \
1170 sdecode(sptr, osub2)
1172#define snextself(osub) \
1173 sptr = (osub).ss[1 - (osub).ssorient]; \
1179#define sorg(osub, vertexptr) \
1180 vertexptr = (vertex) (osub).ss[2 + (osub).ssorient]
1182#define sdest(osub, vertexptr) \
1183 vertexptr = (vertex) (osub).ss[3 - (osub).ssorient]
1185#define setsorg(osub, vertexptr) \
1186 (osub).ss[2 + (osub).ssorient] = (subseg) vertexptr
1188#define setsdest(osub, vertexptr) \
1189 (osub).ss[3 - (osub).ssorient] = (subseg) vertexptr
1191#define segorg(osub, vertexptr) \
1192 vertexptr = (vertex) (osub).ss[4 + (osub).ssorient]
1194#define segdest(osub, vertexptr) \
1195 vertexptr = (vertex) (osub).ss[5 - (osub).ssorient]
1197#define setsegorg(osub, vertexptr) \
1198 (osub).ss[4 + (osub).ssorient] = (subseg) vertexptr
1200#define setsegdest(osub, vertexptr) \
1201 (osub).ss[5 - (osub).ssorient] = (subseg) vertexptr
1207#define mark(osub) (* (int *) ((osub).ss + 8))
1209#define setmark(osub, value) \
1210 * (int *) ((osub).ss + 8) = value
1214#define sbond(osub1, osub2) \
1215 (osub1).ss[(osub1).ssorient] = sencode(osub2); \
1216 (osub2).ss[(osub2).ssorient] = sencode(osub1)
1221#define sdissolve(osub) \
1222 (osub).ss[(osub).ssorient] = (subseg) m->dummysub
1226#define subsegcopy(osub1, osub2) \
1227 (osub2).ss = (osub1).ss; \
1228 (osub2).ssorient = (osub1).ssorient
1232#define subsegequal(osub1, osub2) \
1233 (((osub1).ss == (osub2).ss) && \
1234 ((osub1).ssorient == (osub2).ssorient))
1241#define deadsubseg(sub) ((sub)[1] == (subseg) NULL)
1243#define killsubseg(sub) \
1244 (sub)[1] = (subseg) NULL; \
1245 (sub)[2] = (subseg) NULL
1253#define tspivot(otri, osub) \
1254 sptr = (subseg) (otri).tri[6 + (otri).orient]; \
1260#define stpivot(osub, otri) \
1261 ptr = (triangle) (osub).ss[6 + (osub).ssorient]; \
1266#define tsbond(otri, osub) \
1267 (otri).tri[6 + (otri).orient] = (triangle) sencode(osub); \
1268 (osub).ss[6 + (osub).ssorient] = (subseg) encode(otri)
1272#define tsdissolve(otri) \
1273 (otri).tri[6 + (otri).orient] = (triangle) m->dummysub
1277#define stdissolve(osub) \
1278 (osub).ss[6 + (osub).ssorient] = (subseg) m->dummytri
1284#define vertexmark(vx) ((int *) (vx))[m->vertexmarkindex]
1286#define setvertexmark(vx, value) \
1287 ((int *) (vx))[m->vertexmarkindex] = value
1289#define vertextype(vx) ((int *) (vx))[m->vertexmarkindex + 1]
1291#define setvertextype(vx, value) \
1292 ((int *) (vx))[m->vertexmarkindex + 1] = value
1294#define vertex2tri(vx) ((triangle *) (vx))[m->vertex2triindex]
1296#define setvertex2tri(vx, value) \
1297 ((triangle *) (vx))[m->vertex2triindex] = value
1336#ifdef ANSI_DECLARATORS
1347 REAL dxoa, dxda, dxod;
1348 REAL dyoa, dyda, dyod;
1349 REAL oalen, dalen, odlen;
1354 dxoa = triorg[0] - triapex[0];
1355 dyoa = triorg[1] - triapex[1];
1356 dxda = tridest[0] - triapex[0];
1357 dyda = tridest[1] - triapex[1];
1358 dxod = triorg[0] - tridest[0];
1359 dyod = triorg[1] - tridest[1];
1361 oalen = dxoa * dxoa + dyoa * dyoa;
1362 dalen = dxda * dxda + dyda * dyda;
1363 odlen = dxod * dxod + dyod * dyod;
1365 maxlen = (dalen > oalen) ? dalen : oalen;
1366 maxlen = (odlen > maxlen) ? odlen : maxlen;
1368 if (maxlen > 0.05 * (triorg[0] * triorg[0] + triorg[1] * triorg[1]) + 0.02) {
1385#ifdef ANSI_DECLARATORS
1396#ifdef ANSI_DECLARATORS
1407 if (memptr == (
VOID *) NULL) {
1408 printf(
"Error: Out of memory.\n");
1414#ifdef ANSI_DECLARATORS
1445 printf(
"triangle [-pAcjevngBPNEIOXzo_lQVh] input_file\n");
1447 printf(
"triangle [-pAcjevngBPNEIOXzo_iFlCQVh] input_file\n");
1451 printf(
"triangle [-prq__a__uAcDjevngBPNEIOXzo_YS__lQVh] input_file\n");
1453 printf(
"triangle [-prq__a__uAcDjevngBPNEIOXzo_YS__iFlsCQVh] input_file\n");
1457 printf(
" -p Triangulates a Planar Straight Line Graph (.poly file).\n");
1459 printf(
" -r Refines a previously generated mesh.\n");
1461 " -q Quality mesh generation. A minimum angle may be specified.\n");
1462 printf(
" -a Applies a maximum triangle area constraint.\n");
1463 printf(
" -u Applies a user-defined triangle constraint.\n");
1466 " -A Applies attributes to identify triangles in certain regions.\n");
1467 printf(
" -c Encloses the convex hull with segments.\n");
1469 printf(
" -D Conforming Delaunay: all triangles are truly Delaunay.\n");
1475 printf(
" -j Jettison unused vertices from output .node file.\n");
1476 printf(
" -e Generates an edge list.\n");
1477 printf(
" -v Generates a Voronoi diagram.\n");
1478 printf(
" -n Generates a list of triangle neighbors.\n");
1479 printf(
" -g Generates an .off file for Geomview.\n");
1480 printf(
" -B Suppresses output of boundary information.\n");
1481 printf(
" -P Suppresses output of .poly file.\n");
1482 printf(
" -N Suppresses output of .node file.\n");
1483 printf(
" -E Suppresses output of .ele file.\n");
1484 printf(
" -I Suppresses mesh iteration numbers.\n");
1485 printf(
" -O Ignores holes in .poly file.\n");
1486 printf(
" -X Suppresses use of exact arithmetic.\n");
1487 printf(
" -z Numbers all items starting from zero (rather than one).\n");
1488 printf(
" -o2 Generates second-order subparametric elements.\n");
1490 printf(
" -Y Suppresses boundary segment splitting.\n");
1491 printf(
" -S Specifies maximum number of added Steiner points.\n");
1494 printf(
" -i Uses incremental method, rather than divide-and-conquer.\n");
1495 printf(
" -F Uses Fortune's sweepline algorithm, rather than d-and-c.\n");
1497 printf(
" -l Uses vertical cuts only, rather than alternating cuts.\n");
1501 " -s Force segments into mesh by splitting (instead of using CDT).\n");
1503 printf(
" -C Check consistency of final mesh.\n");
1505 printf(
" -Q Quiet: No terminal output except errors.\n");
1506 printf(
" -V Verbose: Detailed information on what I'm doing.\n");
1507 printf(
" -h Help: Detailed instructions for Triangle.\n");
1523 printf(
"Triangle\n");
1525"A Two-Dimensional Quality Mesh Generator and Delaunay Triangulator.\n");
1526 printf(
"Version 1.6\n\n");
1528"Copyright 1993, 1995, 1997, 1998, 2002, 2005 Jonathan Richard Shewchuk\n");
1529 printf(
"2360 Woolsey #H / Berkeley, California 94705-1927\n");
1530 printf(
"Bugs/comments to jrs@cs.berkeley.edu\n");
1532"Created as part of the Quake project (tools for earthquake simulation).\n");
1534"Supported in part by NSF Grant CMS-9318163 and an NSERC 1967 Scholarship.\n");
1535 printf(
"There is no warranty whatsoever. Use at your own risk.\n");
1537 printf(
"This executable is compiled for single precision arithmetic.\n\n\n");
1539 printf(
"This executable is compiled for double precision arithmetic.\n\n\n");
1542"Triangle generates exact Delaunay triangulations, constrained Delaunay\n");
1544"triangulations, conforming Delaunay triangulations, Voronoi diagrams, and\n");
1546"high-quality triangular meshes. The latter can be generated with no small\n"
1549"or large angles, and are thus suitable for finite element analysis. If no\n"
1552"command line switch is specified, your .node input file is read, and the\n");
1554"Delaunay triangulation is returned in .node and .ele output files. The\n");
1555 printf(
"command syntax is:\n\n");
1556 printf(
"triangle [-prq__a__uAcDjevngBPNEIOXzo_YS__iFlsCQVh] input_file\n\n");
1558"Underscores indicate that numbers may optionally follow certain switches.\n");
1560"Do not leave any space between a switch and its numeric parameter.\n");
1562"input_file must be a file with extension .node, or extension .poly if the\n");
1564"-p switch is used. If -r is used, you must supply .node and .ele files,\n");
1566"and possibly a .poly file and an .area file as well. The formats of these\n"
1568 printf(
"files are described below.\n\n");
1569 printf(
"Command Line Switches:\n\n");
1571" -p Reads a Planar Straight Line Graph (.poly file), which can specify\n"
1574" vertices, segments, holes, regional attributes, and regional area\n");
1576" constraints. Generates a constrained Delaunay triangulation (CDT)\n"
1579" fitting the input; or, if -s, -q, -a, or -u is used, a conforming\n");
1581" constrained Delaunay triangulation (CCDT). If you want a truly\n");
1583" Delaunay (not just constrained Delaunay) triangulation, use -D as\n");
1585" well. When -p is not used, Triangle reads a .node file by default.\n"
1588" -r Refines a previously generated mesh. The mesh is read from a .node\n"
1591" file and an .ele file. If -p is also used, a .poly file is read\n");
1593" and used to constrain segments in the mesh. If -a is also used\n");
1595" (with no number following), an .area file is read and used to\n");
1597" impose area constraints on the mesh. Further details on refinement\n"
1599 printf(
" appear below.\n");
1601" -q Quality mesh generation by Delaunay refinement (a hybrid of Paul\n");
1603" Chew's and Jim Ruppert's algorithms). Adds vertices to the mesh to\n"
1606" ensure that all angles are between 20 and 140 degrees. An\n");
1608" alternative bound on the minimum angle, replacing 20 degrees, may\n");
1610" be specified after the `q'. The specified angle may include a\n");
1612" decimal point, but not exponential notation. Note that a bound of\n"
1615" theta degrees on the smallest angle also implies a bound of\n");
1617" (180 - 2 theta) on the largest angle. If the minimum angle is 28.6\n"
1620" degrees or smaller, Triangle is mathematically guaranteed to\n");
1622" terminate (assuming infinite precision arithmetic--Triangle may\n");
1624" fail to terminate if you run out of precision). In practice,\n");
1626" Triangle often succeeds for minimum angles up to 34 degrees. For\n");
1628" some meshes, however, you might need to reduce the minimum angle to\n"
1631" avoid problems associated with insufficient floating-point\n");
1632 printf(
" precision.\n");
1634" -a Imposes a maximum triangle area. If a number follows the `a', no\n");
1636" triangle is generated whose area is larger than that number. If no\n"
1639" number is specified, an .area file (if -r is used) or .poly file\n");
1641" (if -r is not used) specifies a set of maximum area constraints.\n");
1643" An .area file contains a separate area constraint for each\n");
1645" triangle, and is useful for refining a finite element mesh based on\n"
1648" a posteriori error estimates. A .poly file can optionally contain\n"
1651" an area constraint for each segment-bounded region, thereby\n");
1653" controlling triangle densities in a first triangulation of a PSLG.\n"
1656" You can impose both a fixed area constraint and a varying area\n");
1658" constraint by invoking the -a switch twice, once with and once\n");
1660" without a number following. Each area specified may include a\n");
1661 printf(
" decimal point.\n");
1663" -u Imposes a user-defined constraint on triangle size. There are two\n"
1666" ways to use this feature. One is to edit the triunsuitable()\n");
1668" procedure in triangle.c to encode any constraint you like, then\n");
1670" recompile Triangle. The other is to compile triangle.c with the\n");
1672" EXTERNAL_TEST symbol set (compiler switch -DEXTERNAL_TEST), then\n");
1674" link Triangle with a separate object file that implements\n");
1676" triunsuitable(). In either case, the -u switch causes the user-\n");
1677 printf(
" defined test to be applied to every triangle.\n");
1679" -A Assigns an additional floating-point attribute to each triangle\n");
1681" that identifies what segment-bounded region each triangle belongs\n");
1683" to. Attributes are assigned to regions by the .poly file. If a\n");
1685" region is not explicitly marked by the .poly file, triangles in\n");
1687" that region are assigned an attribute of zero. The -A switch has\n");
1689" an effect only when the -p switch is used and the -r switch is not.\n"
1692" -c Creates segments on the convex hull of the triangulation. If you\n");
1694" are triangulating a vertex set, this switch causes a .poly file to\n"
1697" be written, containing all edges of the convex hull. If you are\n");
1699" triangulating a PSLG, this switch specifies that the whole convex\n");
1701" hull of the PSLG should be triangulated, regardless of what\n");
1703" segments the PSLG has. If you do not use this switch when\n");
1705" triangulating a PSLG, Triangle assumes that you have identified the\n"
1708" region to be triangulated by surrounding it with segments of the\n");
1710" input PSLG. Beware: if you are not careful, this switch can cause\n"
1713" the introduction of an extremely thin angle between a PSLG segment\n"
1716" and a convex hull segment, which can cause overrefinement (and\n");
1718" possibly failure if Triangle runs out of precision). If you are\n");
1720" refining a mesh, the -c switch works differently: it causes a\n");
1722" .poly file to be written containing the boundary edges of the mesh\n"
1724 printf(
" (useful if no .poly file was read).\n");
1726" -D Conforming Delaunay triangulation: use this switch if you want to\n"
1729" ensure that all the triangles in the mesh are Delaunay, and not\n");
1731" merely constrained Delaunay; or if you want to ensure that all the\n"
1734" Voronoi vertices lie within the triangulation. (Some finite volume\n"
1737" methods have this requirement.) This switch invokes Ruppert's\n");
1739" original algorithm, which splits every subsegment whose diametral\n");
1741" circle is encroached. It usually increases the number of vertices\n"
1743 printf(
" and triangles.\n");
1745" -j Jettisons vertices that are not part of the final triangulation\n");
1747" from the output .node file. By default, Triangle copies all\n");
1749" vertices in the input .node file to the output .node file, in the\n");
1751" same order, so their indices do not change. The -j switch prevents\n"
1754" duplicated input vertices, or vertices `eaten' by holes, from\n");
1756" appearing in the output .node file. Thus, if two input vertices\n");
1758" have exactly the same coordinates, only the first appears in the\n");
1760" output. If any vertices are jettisoned, the vertex numbering in\n");
1762" the output .node file differs from that of the input .node file.\n");
1764" -e Outputs (to an .edge file) a list of edges of the triangulation.\n");
1766" -v Outputs the Voronoi diagram associated with the triangulation.\n");
1768" Does not attempt to detect degeneracies, so some Voronoi vertices\n");
1770" may be duplicated. See the discussion of Voronoi diagrams below.\n");
1772" -n Outputs (to a .neigh file) a list of triangles neighboring each\n");
1773 printf(
" triangle.\n");
1775" -g Outputs the mesh to an Object File Format (.off) file, suitable for\n"
1777 printf(
" viewing with the Geometry Center's Geomview package.\n");
1779" -B No boundary markers in the output .node, .poly, and .edge output\n");
1781" files. See the detailed discussion of boundary markers below.\n");
1783" -P No output .poly file. Saves disk space, but you lose the ability\n");
1785" to maintain constraining segments on later refinements of the mesh.\n"
1787 printf(
" -N No output .node file.\n");
1788 printf(
" -E No output .ele file.\n");
1790" -I No iteration numbers. Suppresses the output of .node and .poly\n");
1792" files, so your input files won't be overwritten. (If your input is\n"
1795" a .poly file only, a .node file is written.) Cannot be used with\n");
1797" the -r switch, because that would overwrite your input .ele file.\n");
1799" Shouldn't be used with the -q, -a, -u, or -s switch if you are\n");
1801" using a .node file for input, because no .node file is written, so\n"
1803 printf(
" there is no record of any added Steiner points.\n");
1804 printf(
" -O No holes. Ignores the holes in the .poly file.\n");
1806" -X No exact arithmetic. Normally, Triangle uses exact floating-point\n"
1809" arithmetic for certain tests if it thinks the inexact tests are not\n"
1812" accurate enough. Exact arithmetic ensures the robustness of the\n");
1814" triangulation algorithms, despite floating-point roundoff error.\n");
1816" Disabling exact arithmetic with the -X switch causes a small\n");
1818" improvement in speed and creates the possibility that Triangle will\n"
1820 printf(
" fail to produce a valid mesh. Not recommended.\n");
1822" -z Numbers all items starting from zero (rather than one). Note that\n"
1825" this switch is normally overridden by the value used to number the\n"
1828" first vertex of the input .node or .poly file. However, this\n");
1830" switch is useful when calling Triangle from another program.\n");
1832" -o2 Generates second-order subparametric elements with six nodes each.\n"
1835" -Y No new vertices on the boundary. This switch is useful when the\n");
1837" mesh boundary must be preserved so that it conforms to some\n");
1839" adjacent mesh. Be forewarned that you will probably sacrifice much\n"
1842" of the quality of the mesh; Triangle will try, but the resulting\n");
1844" mesh may contain poorly shaped triangles. Works well if all the\n");
1846" boundary vertices are closely spaced. Specify this switch twice\n");
1848" (`-YY') to prevent all segment splitting, including internal\n");
1849 printf(
" boundaries.\n");
1851" -S Specifies the maximum number of Steiner points (vertices that are\n");
1853" not in the input, but are added to meet the constraints on minimum\n"
1856" angle and maximum area). The default is to allow an unlimited\n");
1858" number. If you specify this switch with no number after it,\n");
1860" the limit is set to zero. Triangle always adds vertices at segment\n"
1863" intersections, even if it needs to use more vertices than the limit\n"
1866" you set. When Triangle inserts segments by splitting (-s), it\n");
1868" always adds enough vertices to ensure that all the segments of the\n"
1870 printf(
" PLSG are recovered, ignoring the limit if necessary.\n");
1872" -i Uses an incremental rather than a divide-and-conquer algorithm to\n");
1874" construct a Delaunay triangulation. Try it if the divide-and-\n");
1875 printf(
" conquer algorithm fails.\n");
1877" -F Uses Steven Fortune's sweepline algorithm to construct a Delaunay\n");
1879" triangulation. Warning: does not use exact arithmetic for all\n");
1880 printf(
" calculations. An exact result is not guaranteed.\n");
1882" -l Uses only vertical cuts in the divide-and-conquer algorithm. By\n");
1884" default, Triangle alternates between vertical and horizontal cuts,\n"
1887" which usually improve the speed except with vertex sets that are\n");
1889" small or short and wide. This switch is primarily of theoretical\n");
1890 printf(
" interest.\n");
1892" -s Specifies that segments should be forced into the triangulation by\n"
1895" recursively splitting them at their midpoints, rather than by\n");
1897" generating a constrained Delaunay triangulation. Segment splitting\n"
1900" is true to Ruppert's original algorithm, but can create needlessly\n"
1903" small triangles. This switch is primarily of theoretical interest.\n"
1906" -C Check the consistency of the final mesh. Uses exact arithmetic for\n"
1909" checking, even if the -X switch is used. Useful if you suspect\n");
1910 printf(
" Triangle is buggy.\n");
1912" -Q Quiet: Suppresses all explanation of what Triangle is doing,\n");
1913 printf(
" unless an error occurs.\n");
1915" -V Verbose: Gives detailed information about what Triangle is doing.\n"
1918" Add more `V's for increasing amount of detail. `-V' is most\n");
1920" useful; itgives information on algorithmic progress and much more\n");
1922" detailed statistics. `-VV' gives vertex-by-vertex details, and\n");
1924" prints so much that Triangle runs much more slowly. `-VVVV' gives\n"
1926 printf(
" information only a debugger could love.\n");
1927 printf(
" -h Help: Displays these instructions.\n");
1929 printf(
"Definitions:\n");
1932" A Delaunay triangulation of a vertex set is a triangulation whose\n");
1934" vertices are the vertex set, that covers the convex hull of the vertex\n");
1936" set. A Delaunay triangulation has the property that no vertex lies\n");
1938" inside the circumscribing circle (circle that passes through all three\n");
1939 printf(
" vertices) of any triangle in the triangulation.\n\n");
1941" A Voronoi diagram of a vertex set is a subdivision of the plane into\n");
1943" polygonal cells (some of which may be unbounded, meaning infinitely\n");
1945" large), where each cell is the set of points in the plane that are closer\n"
1948" to some input vertex than to any other input vertex. The Voronoi diagram\n"
1950 printf(
" is a geometric dual of the Delaunay triangulation.\n\n");
1952" A Planar Straight Line Graph (PSLG) is a set of vertices and segments.\n");
1954" Segments are simply edges, whose endpoints are all vertices in the PSLG.\n"
1957" Segments may intersect each other only at their endpoints. The file\n");
1958 printf(
" format for PSLGs (.poly files) is described below.\n\n");
1960" A constrained Delaunay triangulation (CDT) of a PSLG is similar to a\n");
1962" Delaunay triangulation, but each PSLG segment is present as a single edge\n"
1965" of the CDT. (A constrained Delaunay triangulation is not truly a\n");
1967" Delaunay triangulation, because some of its triangles might not be\n");
1969" Delaunay.) By definition, a CDT does not have any vertices other than\n");
1971" those specified in the input PSLG. Depending on context, a CDT might\n");
1973" cover the convex hull of the PSLG, or it might cover only a segment-\n");
1974 printf(
" bounded region (e.g. a polygon).\n\n");
1976" A conforming Delaunay triangulation of a PSLG is a triangulation in which\n"
1979" each triangle is truly Delaunay, and each PSLG segment is represented by\n"
1982" a linear contiguous sequence of edges of the triangulation. New vertices\n"
1985" (not part of the PSLG) may appear, and each input segment may have been\n");
1987" subdivided into shorter edges (subsegments) by these additional vertices.\n"
1990" The new vertices are frequently necessary to maintain the Delaunay\n");
1991 printf(
" property while ensuring that every segment is represented.\n\n");
1993" A conforming constrained Delaunay triangulation (CCDT) of a PSLG is a\n");
1995" triangulation of a PSLG whose triangles are constrained Delaunay. New\n");
1996 printf(
" vertices may appear, and input segments may be subdivided into\n");
1998" subsegments, but not to guarantee that segments are respected; rather, to\n"
2001" improve the quality of the triangles. The high-quality meshes produced\n");
2003" by the -q switch are usually CCDTs, but can be made conforming Delaunay\n");
2004 printf(
" with the -D switch.\n\n");
2005 printf(
"File Formats:\n\n");
2007" All files may contain comments prefixed by the character '#'. Vertices,\n"
2010" triangles, edges, holes, and maximum area constraints must be numbered\n");
2012" consecutively, starting from either 1 or 0. Whichever you choose, all\n");
2014" input files must be consistent; if the vertices are numbered from 1, so\n");
2016" must be all other objects. Triangle automatically detects your choice\n");
2018" while reading the .node (or .poly) file. (When calling Triangle from\n");
2020" another program, use the -z switch if you wish to number objects from\n");
2021 printf(
" zero.) Examples of these file formats are given below.\n\n");
2022 printf(
" .node files:\n");
2024" First line: <# of vertices> <dimension (must be 2)> <# of attributes>\n"
2027" <# of boundary markers (0 or 1)>\n"
2030" Remaining lines: <vertex #> <x> <y> [attributes] [boundary marker]\n");
2033" The attributes, which are typically floating-point values of physical\n");
2035" quantities (such as mass or conductivity) associated with the nodes of\n"
2038" a finite element mesh, are copied unchanged to the output mesh. If -q,\n"
2041" -a, -u, -D, or -s is selected, each new Steiner point added to the mesh\n"
2043 printf(
" has attributes assigned to it by linear interpolation.\n\n");
2045" If the fourth entry of the first line is `1', the last column of the\n");
2047" remainder of the file is assumed to contain boundary markers. Boundary\n"
2050" markers are used to identify boundary vertices and vertices resting on\n"
2053" PSLG segments; a complete description appears in a section below. The\n"
2056" .node file produced by Triangle contains boundary markers in the last\n");
2057 printf(
" column unless they are suppressed by the -B switch.\n\n");
2058 printf(
" .ele files:\n");
2060" First line: <# of triangles> <nodes per triangle> <# of attributes>\n");
2062" Remaining lines: <triangle #> <node> <node> <node> ... [attributes]\n");
2065" Nodes are indices into the corresponding .node file. The first three\n");
2067" nodes are the corner vertices, and are listed in counterclockwise order\n"
2070" around each triangle. (The remaining nodes, if any, depend on the type\n"
2072 printf(
" of finite element used.)\n\n");
2074" The attributes are just like those of .node files. Because there is no\n"
2077" simple mapping from input to output triangles, Triangle attempts to\n");
2079" interpolate attributes, and may cause a lot of diffusion of attributes\n"
2082" among nearby triangles as the triangulation is refined. Attributes do\n"
2084 printf(
" not diffuse across segments, so attributes used to identify\n");
2085 printf(
" segment-bounded regions remain intact.\n\n");
2087" In .ele files produced by Triangle, each triangular element has three\n");
2089" nodes (vertices) unless the -o2 switch is used, in which case\n");
2091" subparametric quadratic elements with six nodes each are generated.\n");
2093" The first three nodes are the corners in counterclockwise order, and\n");
2095" the fourth, fifth, and sixth nodes lie on the midpoints of the edges\n");
2097" opposite the first, second, and third vertices, respectively.\n");
2099 printf(
" .poly files:\n");
2101" First line: <# of vertices> <dimension (must be 2)> <# of attributes>\n"
2104" <# of boundary markers (0 or 1)>\n"
2107" Following lines: <vertex #> <x> <y> [attributes] [boundary marker]\n");
2108 printf(
" One line: <# of segments> <# of boundary markers (0 or 1)>\n");
2110" Following lines: <segment #> <endpoint> <endpoint> [boundary marker]\n");
2111 printf(
" One line: <# of holes>\n");
2112 printf(
" Following lines: <hole #> <x> <y>\n");
2114" Optional line: <# of regional attributes and/or area constraints>\n");
2116" Optional following lines: <region #> <x> <y> <attribute> <max area>\n");
2119" A .poly file represents a PSLG, as well as some additional information.\n"
2122" The first section lists all the vertices, and is identical to the\n");
2124" format of .node files. <# of vertices> may be set to zero to indicate\n"
2127" that the vertices are listed in a separate .node file; .poly files\n");
2129" produced by Triangle always have this format. A vertex set represented\n"
2132" this way has the advantage that it may easily be triangulated with or\n");
2134" without segments (depending on whether the -p switch is invoked).\n");
2137" The second section lists the segments. Segments are edges whose\n");
2139" presence in the triangulation is enforced. (Depending on the choice of\n"
2142" switches, segment might be subdivided into smaller edges). Each\n");
2144" segment is specified by listing the indices of its two endpoints. This\n"
2147" means that you must include its endpoints in the vertex list. Each\n");
2148 printf(
" segment, like each point, may have a boundary marker.\n\n");
2150" If -q, -a, -u, and -s are not selected, Triangle produces a constrained\n"
2153" Delaunay triangulation (CDT), in which each segment appears as a single\n"
2156" edge in the triangulation. If -q, -a, -u, or -s is selected, Triangle\n"
2159" produces a conforming constrained Delaunay triangulation (CCDT), in\n");
2161" which segments may be subdivided into smaller edges. If -D is\n");
2163" selected, Triangle produces a conforming Delaunay triangulation, so\n");
2165" that every triangle is Delaunay, and not just constrained Delaunay.\n");
2168" The third section lists holes (and concavities, if -c is selected) in\n");
2170" the triangulation. Holes are specified by identifying a point inside\n");
2172" each hole. After the triangulation is formed, Triangle creates holes\n");
2174" by eating triangles, spreading out from each hole point until its\n");
2176" progress is blocked by segments in the PSLG. You must be careful to\n");
2178" enclose each hole in segments, or your whole triangulation might be\n");
2180" eaten away. If the two triangles abutting a segment are eaten, the\n");
2182" segment itself is also eaten. Do not place a hole directly on a\n");
2183 printf(
" segment; if you do, Triangle chooses one side of the segment\n");
2184 printf(
" arbitrarily.\n\n");
2186" The optional fourth section lists regional attributes (to be assigned\n");
2188" to all triangles in a region) and regional constraints on the maximum\n");
2190" triangle area. Triangle reads this section only if the -A switch is\n");
2192" used or the -a switch is used without a number following it, and the -r\n"
2195" switch is not used. Regional attributes and area constraints are\n");
2197" propagated in the same manner as holes: you specify a point for each\n");
2199" attribute and/or constraint, and the attribute and/or constraint\n");
2201" affects the whole region (bounded by segments) containing the point.\n");
2203" If two values are written on a line after the x and y coordinate, the\n");
2205" first such value is assumed to be a regional attribute (but is only\n");
2207" applied if the -A switch is selected), and the second value is assumed\n"
2210" to be a regional area constraint (but is only applied if the -a switch\n"
2213" is selected). You may specify just one value after the coordinates,\n");
2215" which can serve as both an attribute and an area constraint, depending\n"
2218" on the choice of switches. If you are using the -A and -a switches\n");
2220" simultaneously and wish to assign an attribute to some region without\n");
2221 printf(
" imposing an area constraint, use a negative maximum area.\n\n");
2223" When a triangulation is created from a .poly file, you must either\n");
2225" enclose the entire region to be triangulated in PSLG segments, or\n");
2227" use the -c switch, which automatically creates extra segments that\n");
2229" enclose the convex hull of the PSLG. If you do not use the -c switch,\n"
2232" Triangle eats all triangles that are not enclosed by segments; if you\n");
2234" are not careful, your whole triangulation may be eaten away. If you do\n"
2237" use the -c switch, you can still produce concavities by the appropriate\n"
2240" placement of holes just inside the boundary of the convex hull.\n");
2243" An ideal PSLG has no intersecting segments, nor any vertices that lie\n");
2245" upon segments (except, of course, the endpoints of each segment). You\n"
2248" aren't required to make your .poly files ideal, but you should be aware\n"
2251" of what can go wrong. Segment intersections are relatively safe--\n");
2253" Triangle calculates the intersection points for you and adds them to\n");
2255" the triangulation--as long as your machine's floating-point precision\n");
2257" doesn't become a problem. You are tempting the fates if you have three\n"
2260" segments that cross at the same location, and expect Triangle to figure\n"
2263" out where the intersection point is. Thanks to floating-point roundoff\n"
2266" error, Triangle will probably decide that the three segments intersect\n"
2269" at three different points, and you will find a minuscule triangle in\n");
2271" your output--unless Triangle tries to refine the tiny triangle, uses\n");
2273" up the last bit of machine precision, and fails to terminate at all.\n");
2275" You're better off putting the intersection point in the input files,\n");
2277" and manually breaking up each segment into two. Similarly, if you\n");
2279" place a vertex at the middle of a segment, and hope that Triangle will\n"
2282" break up the segment at that vertex, you might get lucky. On the other\n"
2285" hand, Triangle might decide that the vertex doesn't lie precisely on\n");
2287" the segment, and you'll have a needle-sharp triangle in your output--or\n"
2289 printf(
" a lot of tiny triangles if you're generating a quality mesh.\n");
2292" When Triangle reads a .poly file, it also writes a .poly file, which\n");
2294" includes all the subsegments--the edges that are parts of input\n");
2296" segments. If the -c switch is used, the output .poly file also\n");
2298" includes all of the edges on the convex hull. Hence, the output .poly\n"
2301" file is useful for finding edges associated with input segments and for\n"
2304" setting boundary conditions in finite element simulations. Moreover,\n");
2306" you will need the output .poly file if you plan to refine the output\n");
2308" mesh, and don't want segments to be missing in later triangulations.\n");
2310 printf(
" .area files:\n");
2311 printf(
" First line: <# of triangles>\n");
2312 printf(
" Following lines: <triangle #> <maximum area>\n");
2315" An .area file associates with each triangle a maximum area that is used\n"
2318" for mesh refinement. As with other file formats, every triangle must\n");
2320" be represented, and the triangles must be numbered consecutively. A\n");
2322" triangle may be left unconstrained by assigning it a negative maximum\n");
2323 printf(
" area.\n\n");
2324 printf(
" .edge files:\n");
2325 printf(
" First line: <# of edges> <# of boundary markers (0 or 1)>\n");
2327" Following lines: <edge #> <endpoint> <endpoint> [boundary marker]\n");
2330" Endpoints are indices into the corresponding .node file. Triangle can\n"
2333" produce .edge files (use the -e switch), but cannot read them. The\n");
2335" optional column of boundary markers is suppressed by the -B switch.\n");
2338" In Voronoi diagrams, one also finds a special kind of edge that is an\n");
2340" infinite ray with only one endpoint. For these edges, a different\n");
2341 printf(
" format is used:\n\n");
2342 printf(
" <edge #> <endpoint> -1 <direction x> <direction y>\n\n");
2344" The `direction' is a floating-point vector that indicates the direction\n"
2346 printf(
" of the infinite ray.\n\n");
2347 printf(
" .neigh files:\n");
2349" First line: <# of triangles> <# of neighbors per triangle (always 3)>\n"
2352" Following lines: <triangle #> <neighbor> <neighbor> <neighbor>\n");
2355" Neighbors are indices into the corresponding .ele file. An index of -1\n"
2358" indicates no neighbor (because the triangle is on an exterior\n");
2360" boundary). The first neighbor of triangle i is opposite the first\n");
2361 printf(
" corner of triangle i, and so on.\n\n");
2363" Triangle can produce .neigh files (use the -n switch), but cannot read\n"
2365 printf(
" them.\n\n");
2366 printf(
"Boundary Markers:\n\n");
2368" Boundary markers are tags used mainly to identify which output vertices\n");
2370" and edges are associated with which PSLG segment, and to identify which\n");
2372" vertices and edges occur on a boundary of the triangulation. A common\n");
2374" use is to determine where boundary conditions should be applied to a\n");
2376" finite element mesh. You can prevent boundary markers from being written\n"
2378 printf(
" into files produced by Triangle by using the -B switch.\n\n");
2380" The boundary marker associated with each segment in an output .poly file\n"
2382 printf(
" and each edge in an output .edge file is chosen as follows:\n");
2384" - If an output edge is part or all of a PSLG segment with a nonzero\n");
2386" boundary marker, then the edge is assigned the same marker.\n");
2388" - Otherwise, if the edge lies on a boundary of the triangulation\n");
2390" (even the boundary of a hole), then the edge is assigned the marker\n");
2391 printf(
" one (1).\n");
2392 printf(
" - Otherwise, the edge is assigned the marker zero (0).\n");
2394" The boundary marker associated with each vertex in an output .node file\n");
2395 printf(
" is chosen as follows:\n");
2397" - If a vertex is assigned a nonzero boundary marker in the input file,\n"
2400" then it is assigned the same marker in the output .node file.\n");
2402" - Otherwise, if the vertex lies on a PSLG segment (even if it is an\n");
2404" endpoint of the segment) with a nonzero boundary marker, then the\n");
2406" vertex is assigned the same marker. If the vertex lies on several\n");
2407 printf(
" such segments, one of the markers is chosen arbitrarily.\n");
2409" - Otherwise, if the vertex occurs on a boundary of the triangulation,\n");
2410 printf(
" then the vertex is assigned the marker one (1).\n");
2411 printf(
" - Otherwise, the vertex is assigned the marker zero (0).\n");
2414" If you want Triangle to determine for you which vertices and edges are on\n"
2417" the boundary, assign them the boundary marker zero (or use no markers at\n"
2420" all) in your input files. In the output files, all boundary vertices,\n");
2421 printf(
" edges, and segments will be assigned the value one.\n\n");
2422 printf(
"Triangulation Iteration Numbers:\n\n");
2424" Because Triangle can read and refine its own triangulations, input\n");
2426" and output files have iteration numbers. For instance, Triangle might\n");
2428" read the files mesh.3.node, mesh.3.ele, and mesh.3.poly, refine the\n");
2430" triangulation, and output the files mesh.4.node, mesh.4.ele, and\n");
2431 printf(
" mesh.4.poly. Files with no iteration number are treated as if\n");
2433" their iteration number is zero; hence, Triangle might read the file\n");
2435" points.node, triangulate it, and produce the files points.1.node and\n");
2436 printf(
" points.1.ele.\n\n");
2438" Iteration numbers allow you to create a sequence of successively finer\n");
2440" meshes suitable for multigrid methods. They also allow you to produce a\n"
2443" sequence of meshes using error estimate-driven mesh refinement.\n");
2446" If you're not using refinement or quality meshing, and you don't like\n");
2448" iteration numbers, use the -I switch to disable them. This switch also\n");
2450" disables output of .node and .poly files to prevent your input files from\n"
2453" being overwritten. (If the input is a .poly file that contains its own\n");
2455" points, a .node file is written. This can be quite convenient for\n");
2456 printf(
" computing CDTs or quality meshes.)\n\n");
2457 printf(
"Examples of How to Use Triangle:\n\n");
2459" `triangle dots' reads vertices from dots.node, and writes their Delaunay\n"
2462" triangulation to dots.1.node and dots.1.ele. (dots.1.node is identical\n");
2464" to dots.node.) `triangle -I dots' writes the triangulation to dots.ele\n");
2466" instead. (No additional .node file is needed, so none is written.)\n");
2469" `triangle -pe object.1' reads a PSLG from object.1.poly (and possibly\n");
2471" object.1.node, if the vertices are omitted from object.1.poly) and writes\n"
2474" its constrained Delaunay triangulation to object.2.node and object.2.ele.\n"
2477" The segments are copied to object.2.poly, and all edges are written to\n");
2478 printf(
" object.2.edge.\n\n");
2480" `triangle -pq31.5a.1 object' reads a PSLG from object.poly (and possibly\n"
2483" object.node), generates a mesh whose angles are all between 31.5 and 117\n"
2486" degrees and whose triangles all have areas of 0.1 or less, and writes the\n"
2489" mesh to object.1.node and object.1.ele. Each segment may be broken up\n");
2490 printf(
" into multiple subsegments; these are written to object.1.poly.\n");
2493" Here is a sample file `box.poly' describing a square with a square hole:\n"
2497" # A box with eight vertices in 2D, no attributes, one boundary marker.\n"
2499 printf(
" 8 2 0 1\n");
2500 printf(
" # Outer box has these vertices:\n");
2501 printf(
" 1 0 0 0\n");
2502 printf(
" 2 0 3 0\n");
2503 printf(
" 3 3 0 0\n");
2504 printf(
" 4 3 3 33 # A special marker for this vertex.\n");
2505 printf(
" # Inner square has these vertices:\n");
2506 printf(
" 5 1 1 0\n");
2507 printf(
" 6 1 2 0\n");
2508 printf(
" 7 2 1 0\n");
2509 printf(
" 8 2 2 0\n");
2510 printf(
" # Five segments with boundary markers.\n");
2512 printf(
" 1 1 2 5 # Left side of outer box.\n");
2513 printf(
" # Square hole has these segments:\n");
2514 printf(
" 2 5 7 0\n");
2515 printf(
" 3 7 8 0\n");
2516 printf(
" 4 8 6 10\n");
2517 printf(
" 5 6 5 0\n");
2518 printf(
" # One hole in the middle of the inner square.\n");
2520 printf(
" 1 1.5 1.5\n");
2523" Note that some segments are missing from the outer square, so you must\n");
2525" use the `-c' switch. After `triangle -pqc box.poly', here is the output\n"
2528" file `box.1.node', with twelve vertices. The last four vertices were\n");
2530" added to meet the angle constraint. Vertices 1, 2, and 9 have markers\n");
2532" from segment 1. Vertices 6 and 8 have markers from segment 4. All the\n");
2534" other vertices but 4 have been marked to indicate that they lie on a\n");
2535 printf(
" boundary.\n\n");
2536 printf(
" 12 2 0 1\n");
2537 printf(
" 1 0 0 5\n");
2538 printf(
" 2 0 3 5\n");
2539 printf(
" 3 3 0 1\n");
2540 printf(
" 4 3 3 33\n");
2541 printf(
" 5 1 1 1\n");
2542 printf(
" 6 1 2 10\n");
2543 printf(
" 7 2 1 1\n");
2544 printf(
" 8 2 2 10\n");
2545 printf(
" 9 0 1.5 5\n");
2546 printf(
" 10 1.5 0 1\n");
2547 printf(
" 11 3 1.5 1\n");
2548 printf(
" 12 1.5 3 1\n");
2549 printf(
" # Generated by triangle -pqc box.poly\n");
2551 printf(
" Here is the output file `box.1.ele', with twelve triangles.\n");
2553 printf(
" 12 3 0\n");
2554 printf(
" 1 5 6 9\n");
2555 printf(
" 2 10 3 7\n");
2556 printf(
" 3 6 8 12\n");
2557 printf(
" 4 9 1 5\n");
2558 printf(
" 5 6 2 9\n");
2559 printf(
" 6 7 3 11\n");
2560 printf(
" 7 11 4 8\n");
2561 printf(
" 8 7 5 10\n");
2562 printf(
" 9 12 2 6\n");
2563 printf(
" 10 8 7 11\n");
2564 printf(
" 11 5 1 10\n");
2565 printf(
" 12 8 4 12\n");
2566 printf(
" # Generated by triangle -pqc box.poly\n\n");
2568" Here is the output file `box.1.poly'. Note that segments have been added\n"
2571" to represent the convex hull, and some segments have been subdivided by\n");
2573" newly added vertices. Note also that <# of vertices> is set to zero to\n");
2574 printf(
" indicate that the vertices should be read from the .node file.\n");
2576 printf(
" 0 2 0 1\n");
2578 printf(
" 1 1 9 5\n");
2579 printf(
" 2 5 7 1\n");
2580 printf(
" 3 8 7 1\n");
2581 printf(
" 4 6 8 10\n");
2582 printf(
" 5 5 6 1\n");
2583 printf(
" 6 3 10 1\n");
2584 printf(
" 7 4 11 1\n");
2585 printf(
" 8 2 12 1\n");
2586 printf(
" 9 9 2 5\n");
2587 printf(
" 10 10 1 1\n");
2588 printf(
" 11 11 3 1\n");
2589 printf(
" 12 12 4 1\n");
2591 printf(
" 1 1.5 1.5\n");
2592 printf(
" # Generated by triangle -pqc box.poly\n");
2594 printf(
"Refinement and Area Constraints:\n");
2597" The -r switch causes a mesh (.node and .ele files) to be read and\n");
2599" refined. If the -p switch is also used, a .poly file is read and used to\n"
2602" specify edges that are constrained and cannot be eliminated (although\n");
2604" they can be subdivided into smaller edges) by the refinement process.\n");
2607" When you refine a mesh, you generally want to impose tighter constraints.\n"
2610" One way to accomplish this is to use -q with a larger angle, or -a\n");
2612" followed by a smaller area than you used to generate the mesh you are\n");
2614" refining. Another way to do this is to create an .area file, which\n");
2616" specifies a maximum area for each triangle, and use the -a switch\n");
2618" (without a number following). Each triangle's area constraint is applied\n"
2621" to that triangle. Area constraints tend to diffuse as the mesh is\n");
2623" refined, so if there are large variations in area constraint between\n");
2625" adjacent triangles, you may not get the results you want. In that case,\n"
2628" consider instead using the -u switch and writing a C procedure that\n");
2629 printf(
" determines which triangles are too large.\n\n");
2631" If you are refining a mesh composed of linear (three-node) elements, the\n"
2634" output mesh contains all the nodes present in the input mesh, in the same\n"
2637" order, with new nodes added at the end of the .node file. However, the\n");
2639" refinement is not hierarchical: there is no guarantee that each output\n");
2641" element is contained in a single input element. Often, an output element\n"
2644" can overlap two or three input elements, and some input edges are not\n");
2646" present in the output mesh. Hence, a sequence of refined meshes forms a\n"
2649" hierarchy of nodes, but not a hierarchy of elements. If you refine a\n");
2651" mesh of higher-order elements, the hierarchical property applies only to\n"
2654" the nodes at the corners of an element; the midpoint nodes on each edge\n");
2655 printf(
" are discarded before the mesh is refined.\n\n");
2657" Maximum area constraints in .poly files operate differently from those in\n"
2660" .area files. A maximum area in a .poly file applies to the whole\n");
2662" (segment-bounded) region in which a point falls, whereas a maximum area\n");
2664" in an .area file applies to only one triangle. Area constraints in .poly\n"
2667" files are used only when a mesh is first generated, whereas area\n");
2669" constraints in .area files are used only to refine an existing mesh, and\n"
2672" are typically based on a posteriori error estimates resulting from a\n");
2673 printf(
" finite element simulation on that mesh.\n\n");
2675" `triangle -rq25 object.1' reads object.1.node and object.1.ele, then\n");
2677" refines the triangulation to enforce a 25 degree minimum angle, and then\n"
2680" writes the refined triangulation to object.2.node and object.2.ele.\n");
2683" `triangle -rpaa6.2 z.3' reads z.3.node, z.3.ele, z.3.poly, and z.3.area.\n"
2686" After reconstructing the mesh and its subsegments, Triangle refines the\n");
2688" mesh so that no triangle has area greater than 6.2, and furthermore the\n");
2690" triangles satisfy the maximum area constraints in z.3.area. No angle\n");
2692" bound is imposed at all. The output is written to z.4.node, z.4.ele, and\n"
2694 printf(
" z.4.poly.\n\n");
2696" The sequence `triangle -qa1 x', `triangle -rqa.3 x.1', `triangle -rqa.1\n");
2698" x.2' creates a sequence of successively finer meshes x.1, x.2, and x.3,\n");
2699 printf(
" suitable for multigrid.\n\n");
2700 printf(
"Convex Hulls and Mesh Boundaries:\n\n");
2702" If the input is a vertex set (not a PSLG), Triangle produces its convex\n");
2704" hull as a by-product in the output .poly file if you use the -c switch.\n");
2706" There are faster algorithms for finding a two-dimensional convex hull\n");
2707 printf(
" than triangulation, of course, but this one comes for free.\n\n");
2709" If the input is an unconstrained mesh (you are using the -r switch but\n");
2711" not the -p switch), Triangle produces a list of its boundary edges\n");
2713" (including hole boundaries) as a by-product when you use the -c switch.\n");
2715" If you also use the -p switch, the output .poly file contains all the\n");
2716 printf(
" segments from the input .poly file as well.\n\n");
2717 printf(
"Voronoi Diagrams:\n\n");
2719" The -v switch produces a Voronoi diagram, in files suffixed .v.node and\n");
2721" .v.edge. For example, `triangle -v points' reads points.node, produces\n");
2723" its Delaunay triangulation in points.1.node and points.1.ele, and\n");
2725" produces its Voronoi diagram in points.1.v.node and points.1.v.edge. The\n"
2728" .v.node file contains a list of all Voronoi vertices, and the .v.edge\n");
2730" file contains a list of all Voronoi edges, some of which may be infinite\n"
2733" rays. (The choice of filenames makes it easy to run the set of Voronoi\n");
2734 printf(
" vertices through Triangle, if so desired.)\n\n");
2736" This implementation does not use exact arithmetic to compute the Voronoi\n"
2739" vertices, and does not check whether neighboring vertices are identical.\n"
2742" Be forewarned that if the Delaunay triangulation is degenerate or\n");
2744" near-degenerate, the Voronoi diagram may have duplicate vertices or\n");
2745 printf(
" crossing edges.\n\n");
2747" The result is a valid Voronoi diagram only if Triangle's output is a true\n"
2750" Delaunay triangulation. The Voronoi output is usually meaningless (and\n");
2752" may contain crossing edges and other pathology) if the output is a CDT or\n"
2755" CCDT, or if it has holes or concavities. If the triangulated domain is\n");
2757" convex and has no holes, you can use -D switch to force Triangle to\n");
2759" construct a conforming Delaunay triangulation instead of a CCDT, so the\n");
2760 printf(
" Voronoi diagram will be valid.\n\n");
2761 printf(
"Mesh Topology:\n\n");
2763" You may wish to know which triangles are adjacent to a certain Delaunay\n");
2765" edge in an .edge file, which Voronoi cells are adjacent to a certain\n");
2767" Voronoi edge in a .v.edge file, or which Voronoi cells are adjacent to\n");
2769" each other. All of this information can be found by cross-referencing\n");
2771" output files with the recollection that the Delaunay triangulation and\n");
2772 printf(
" the Voronoi diagram are planar duals.\n\n");
2774" Specifically, edge i of an .edge file is the dual of Voronoi edge i of\n");
2776" the corresponding .v.edge file, and is rotated 90 degrees counterclock-\n");
2778" wise from the Voronoi edge. Triangle j of an .ele file is the dual of\n");
2780" vertex j of the corresponding .v.node file. Voronoi cell k is the dual\n");
2781 printf(
" of vertex k of the corresponding .node file.\n\n");
2783" Hence, to find the triangles adjacent to a Delaunay edge, look at the\n");
2785" vertices of the corresponding Voronoi edge. If the endpoints of a\n");
2787" Voronoi edge are Voronoi vertices 2 and 6 respectively, then triangles 2\n"
2790" and 6 adjoin the left and right sides of the corresponding Delaunay edge,\n"
2793" respectively. To find the Voronoi cells adjacent to a Voronoi edge, look\n"
2796" at the endpoints of the corresponding Delaunay edge. If the endpoints of\n"
2799" a Delaunay edge are input vertices 7 and 12, then Voronoi cells 7 and 12\n"
2802" adjoin the right and left sides of the corresponding Voronoi edge,\n");
2804" respectively. To find which Voronoi cells are adjacent to each other,\n");
2805 printf(
" just read the list of Delaunay edges.\n\n");
2807" Triangle does not write a list of the edges adjoining each Voronoi cell,\n"
2810" but you can reconstructed it straightforwardly. For instance, to find\n");
2812" all the edges of Voronoi cell 1, search the output .edge file for every\n");
2814" edge that has input vertex 1 as an endpoint. The corresponding dual\n");
2816" edges in the output .v.edge file form the boundary of Voronoi cell 1.\n");
2819" For each Voronoi vertex, the .neigh file gives a list of the three\n");
2821" Voronoi vertices attached to it. You might find this more convenient\n");
2822 printf(
" than the .v.edge file.\n\n");
2823 printf(
"Quadratic Elements:\n\n");
2825" Triangle generates meshes with subparametric quadratic elements if the\n");
2827" -o2 switch is specified. Quadratic elements have six nodes per element,\n"
2830" rather than three. `Subparametric' means that the edges of the triangles\n"
2833" are always straight, so that subparametric quadratic elements are\n");
2835" geometrically identical to linear elements, even though they can be used\n"
2838" with quadratic interpolating functions. The three extra nodes of an\n");
2840" element fall at the midpoints of the three edges, with the fourth, fifth,\n"
2843" and sixth nodes appearing opposite the first, second, and third corners\n");
2844 printf(
" respectively.\n\n");
2845 printf(
"Domains with Small Angles:\n\n");
2847" If two input segments adjoin each other at a small angle, clearly the -q\n"
2850" switch cannot remove the small angle. Moreover, Triangle may have no\n");
2852" choice but to generate additional triangles whose smallest angles are\n");
2854" smaller than the specified bound. However, these triangles only appear\n");
2856" between input segments separated by small angles. Moreover, if you\n");
2858" request a minimum angle of theta degrees, Triangle will generally produce\n"
2861" no angle larger than 180 - 2 theta, even if it is forced to compromise on\n"
2863 printf(
" the minimum angle.\n\n");
2864 printf(
"Statistics:\n\n");
2866" After generating a mesh, Triangle prints a count of entities in the\n");
2868" output mesh, including the number of vertices, triangles, edges, exterior\n"
2871" boundary edges (i.e. subsegments on the boundary of the triangulation,\n");
2873" including hole boundaries), interior boundary edges (i.e. subsegments of\n"
2876" input segments not on the boundary), and total subsegments. If you've\n");
2878" forgotten the statistics for an existing mesh, run Triangle on that mesh\n"
2881" with the -rNEP switches to read the mesh and print the statistics without\n"
2884" writing any files. Use -rpNEP if you've got a .poly file for the mesh.\n");
2887" The -V switch produces extended statistics, including a rough estimate\n");
2889" of memory use, the number of calls to geometric predicates, and\n");
2891" histograms of the angles and the aspect ratios of the triangles in the\n");
2892 printf(
" mesh.\n\n");
2893 printf(
"Exact Arithmetic:\n\n");
2895" Triangle uses adaptive exact arithmetic to perform what computational\n");
2897" geometers call the `orientation' and `incircle' tests. If the floating-\n"
2900" point arithmetic of your machine conforms to the IEEE 754 standard (as\n");
2902" most workstations do), and does not use extended precision internal\n");
2904" floating-point registers, then your output is guaranteed to be an\n");
2906" absolutely true Delaunay or constrained Delaunay triangulation, roundoff\n"
2909" error notwithstanding. The word `adaptive' implies that these arithmetic\n"
2912" routines compute the result only to the precision necessary to guarantee\n"
2915" correctness, so they are usually nearly as fast as their approximate\n");
2916 printf(
" counterparts.\n\n");
2918" May CPUs, including Intel x86 processors, have extended precision\n");
2920" floating-point registers. These must be reconfigured so their precision\n"
2923" is reduced to memory precision. Triangle does this if it is compiled\n");
2924 printf(
" correctly. See the makefile for details.\n\n");
2926" The exact tests can be disabled with the -X switch. On most inputs, this\n"
2929" switch reduces the computation time by about eight percent--it's not\n");
2931" worth the risk. There are rare difficult inputs (having many collinear\n");
2933" and cocircular vertices), however, for which the difference in speed\n");
2935" could be a factor of two. Be forewarned that these are precisely the\n");
2937" inputs most likely to cause errors if you use the -X switch. Hence, the\n"
2939 printf(
" -X switch is not recommended.\n\n");
2941" Unfortunately, the exact tests don't solve every numerical problem.\n");
2943" Exact arithmetic is not used to compute the positions of new vertices,\n");
2945" because the bit complexity of vertex coordinates would grow without\n");
2947" bound. Hence, segment intersections aren't computed exactly; in very\n");
2949" unusual cases, roundoff error in computing an intersection point might\n");
2951" actually lead to an inverted triangle and an invalid triangulation.\n");
2953" (This is one reason to specify your own intersection points in your .poly\n"
2956" files.) Similarly, exact arithmetic is not used to compute the vertices\n"
2958 printf(
" of the Voronoi diagram.\n\n");
2960" Another pair of problems not solved by the exact arithmetic routines is\n");
2962" underflow and overflow. If Triangle is compiled for double precision\n");
2964" arithmetic, I believe that Triangle's geometric predicates work correctly\n"
2967" if the exponent of every input coordinate falls in the range [-148, 201].\n"
2970" Underflow can silently prevent the orientation and incircle tests from\n");
2972" being performed exactly, while overflow typically causes a floating\n");
2973 printf(
" exception.\n\n");
2974 printf(
"Calling Triangle from Another Program:\n\n");
2975 printf(
" Read the file triangle.h for details.\n\n");
2976 printf(
"Troubleshooting:\n\n");
2977 printf(
" Please read this section before mailing me bugs.\n\n");
2978 printf(
" `My output mesh has no triangles!'\n\n");
2980" If you're using a PSLG, you've probably failed to specify a proper set\n"
2983" of bounding segments, or forgotten to use the -c switch. Or you may\n");
2985" have placed a hole badly, thereby eating all your triangles. To test\n");
2986 printf(
" these possibilities, try again with the -c and -O switches.\n");
2988" Alternatively, all your input vertices may be collinear, in which case\n"
2990 printf(
" you can hardly expect to triangulate them.\n\n");
2991 printf(
" `Triangle doesn't terminate, or just crashes.'\n\n");
2993" Bad things can happen when triangles get so small that the distance\n");
2995" between their vertices isn't much larger than the precision of your\n");
2997" machine's arithmetic. If you've compiled Triangle for single-precision\n"
3000" arithmetic, you might do better by recompiling it for double-precision.\n"
3003" Then again, you might just have to settle for more lenient constraints\n"
3006" on the minimum angle and the maximum area than you had planned.\n");
3009" You can minimize precision problems by ensuring that the origin lies\n");
3011" inside your vertex set, or even inside the densest part of your\n");
3013" mesh. If you're triangulating an object whose x-coordinates all fall\n");
3015" between 6247133 and 6247134, you're not leaving much floating-point\n");
3016 printf(
" precision for Triangle to work with.\n\n");
3018" Precision problems can occur covertly if the input PSLG contains two\n");
3020" segments that meet (or intersect) at an extremely small angle, or if\n");
3022" such an angle is introduced by the -c switch. If you don't realize\n");
3024" that a tiny angle is being formed, you might never discover why\n");
3026" Triangle is crashing. To check for this possibility, use the -S switch\n"
3029" (with an appropriate limit on the number of Steiner points, found by\n");
3031" trial-and-error) to stop Triangle early, and view the output .poly file\n"
3034" with Show Me (described below). Look carefully for regions where dense\n"
3037" clusters of vertices are forming and for small angles between segments.\n"
3040" Zoom in closely, as such segments might look like a single segment from\n"
3042 printf(
" a distance.\n\n");
3044" If some of the input values are too large, Triangle may suffer a\n");
3046" floating exception due to overflow when attempting to perform an\n");
3048" orientation or incircle test. (Read the section on exact arithmetic\n");
3050" above.) Again, I recommend compiling Triangle for double (rather\n");
3051 printf(
" than single) precision arithmetic.\n\n");
3053" Unexpected problems can arise if you use quality meshing (-q, -a, or\n");
3055" -u) with an input that is not segment-bounded--that is, if your input\n");
3057" is a vertex set, or you're using the -c switch. If the convex hull of\n"
3060" your input vertices has collinear vertices on its boundary, an input\n");
3062" vertex that you think lies on the convex hull might actually lie just\n");
3064" inside the convex hull. If so, the vertex and the nearby convex hull\n");
3066" edge form an extremely thin triangle. When Triangle tries to refine\n");
3068" the mesh to enforce angle and area constraints, Triangle might generate\n"
3071" extremely tiny triangles, or it might fail because of insufficient\n");
3072 printf(
" floating-point precision.\n\n");
3074" `The numbering of the output vertices doesn't match the input vertices.'\n"
3078" You may have had duplicate input vertices, or you may have eaten some\n");
3080" of your input vertices with a hole, or by placing them outside the area\n"
3083" enclosed by segments. In any case, you can solve the problem by not\n");
3084 printf(
" using the -j switch.\n\n");
3086" `Triangle executes without incident, but when I look at the resulting\n");
3088" mesh, it has overlapping triangles or other geometric inconsistencies.'\n");
3091" If you select the -X switch, Triangle occasionally makes mistakes due\n");
3093" to floating-point roundoff error. Although these errors are rare,\n");
3095" don't use the -X switch. If you still have problems, please report the\n"
3097 printf(
" bug.\n\n");
3099" `Triangle executes without incident, but when I look at the resulting\n");
3100 printf(
" Voronoi diagram, it has overlapping edges or other geometric\n");
3101 printf(
" inconsistencies.'\n");
3104" If your input is a PSLG (-p), you can only expect a meaningful Voronoi\n"
3107" diagram if the domain you are triangulating is convex and free of\n");
3109" holes, and you use the -D switch to construct a conforming Delaunay\n");
3110 printf(
" triangulation (instead of a CDT or CCDT).\n\n");
3112" Strange things can happen if you've taken liberties with your PSLG. Do\n");
3114" you have a vertex lying in the middle of a segment? Triangle sometimes\n");
3116" copes poorly with that sort of thing. Do you want to lay out a collinear\n"
3119" row of evenly spaced, segment-connected vertices? Have you simply\n");
3121" defined one long segment connecting the leftmost vertex to the rightmost\n"
3124" vertex, and a bunch of vertices lying along it? This method occasionally\n"
3127" works, especially with horizontal and vertical lines, but often it\n");
3129" doesn't, and you'll have to connect each adjacent pair of vertices with a\n"
3131 printf(
" separate segment. If you don't like it, tough.\n\n");
3133" Furthermore, if you have segments that intersect other than at their\n");
3135" endpoints, try not to let the intersections fall extremely close to PSLG\n"
3137 printf(
" vertices or each other.\n\n");
3139" If you have problems refining a triangulation not produced by Triangle:\n");
3141" Are you sure the triangulation is geometrically valid? Is it formatted\n");
3143" correctly for Triangle? Are the triangles all listed so the first three\n"
3146" vertices are their corners in counterclockwise order? Are all of the\n");
3148" triangles constrained Delaunay? Triangle's Delaunay refinement algorithm\n"
3150 printf(
" assumes that it starts with a CDT.\n\n");
3151 printf(
"Show Me:\n\n");
3153" Triangle comes with a separate program named `Show Me', whose primary\n");
3155" purpose is to draw meshes on your screen or in PostScript. Its secondary\n"
3158" purpose is to check the validity of your input files, and do so more\n");
3160" thoroughly than Triangle does. Unlike Triangle, Show Me requires that\n");
3162" you have the X Windows system. Sorry, Microsoft Windows users.\n");
3164 printf(
"Triangle on the Web:\n");
3166 printf(
" To see an illustrated version of these instructions, check out\n");
3168 printf(
" http://www.cs.cmu.edu/~quake/triangle.html\n");
3170 printf(
"A Brief Plea:\n");
3173" If you use Triangle, and especially if you use it to accomplish real\n");
3175" work, I would like very much to hear from you. A short letter or email\n");
3177" (to jrs@cs.berkeley.edu) describing how you use Triangle will mean a lot\n"
3180" to me. The more people I know are using this program, the more easily I\n"
3183" can justify spending time on improvements, which in turn will benefit\n");
3185" you. Also, I can put you on a list to receive email whenever a new\n");
3186 printf(
" version of Triangle is available.\n\n");
3188" If you use a mesh generated by Triangle in a publication, please include\n"
3191" an acknowledgment as well. And please spell Triangle with a capital `T'!\n"
3194" If you want to include a citation, use `Jonathan Richard Shewchuk,\n");
3196" ``Triangle: Engineering a 2D Quality Mesh Generator and Delaunay\n");
3198" Triangulator,'' in Applied Computational Geometry: Towards Geometric\n");
3200" Engineering (Ming C. Lin and Dinesh Manocha, editors), volume 1148 of\n");
3202" Lecture Notes in Computer Science, pages 203-222, Springer-Verlag,\n");
3204" Berlin, May 1996. (From the First ACM Workshop on Applied Computational\n"
3206 printf(
" Geometry.)'\n\n");
3207 printf(
"Research credit:\n\n");
3209" Of course, I can take credit for only a fraction of the ideas that made\n");
3211" this mesh generator possible. Triangle owes its existence to the efforts\n"
3214" of many fine computational geometers and other researchers, including\n");
3216" Marshall Bern, L. Paul Chew, Kenneth L. Clarkson, Boris Delaunay, Rex A.\n"
3219" Dwyer, David Eppstein, Steven Fortune, Leonidas J. Guibas, Donald E.\n");
3221" Knuth, Charles L. Lawson, Der-Tsai Lee, Gary L. Miller, Ernst P. Mucke,\n");
3223" Steven E. Pav, Douglas M. Priest, Jim Ruppert, Isaac Saias, Bruce J.\n");
3225" Schachter, Micha Sharir, Peter W. Shor, Daniel D. Sleator, Jorge Stolfi,\n"
3227 printf(
" Robert E. Tarjan, Alper Ungor, Christopher J. Van Wyk, Noel J.\n");
3229" Walkington, and Binhai Zhu. See the comments at the beginning of the\n");
3230 printf(
" source code for references.\n\n");
3244 printf(
" Please report this bug to jrs@cs.berkeley.edu\n");
3245 printf(
" Include the message above, your input data set, and the exact\n");
3246 printf(
" command line you used to run Triangle.\n");
3257#ifdef ANSI_DECLARATORS
3278 b->
poly =
b->refine =
b->quality = 0;
3279 b->vararea =
b->fixedarea =
b->usertest = 0;
3280 b->regionattrib =
b->convex =
b->weighted =
b->jettison = 0;
3282 b->edgesout =
b->voronoi =
b->neighbors =
b->geomview = 0;
3283 b->nobound =
b->nopolywritten =
b->nonodewritten =
b->noelewritten = 0;
3284 b->noiterationnum = 0;
3285 b->noholes =
b->noexact = 0;
3286 b->incremental =
b->sweepline = 0;
3296 b->quiet =
b->verbose = 0;
3298 b->innodefilename[0] =
'\0';
3303 if (argv[i][0] ==
'-') {
3305 for (j =
STARTINDEX; argv[i][j] !=
'\0'; j++) {
3306 if (argv[i][j] ==
'p') {
3310 if (argv[i][j] ==
'r') {
3313 if (argv[i][j] ==
'q') {
3315 if (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
3316 (argv[i][j + 1] ==
'.')) {
3318 while (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
3319 (argv[i][j + 1] ==
'.')) {
3321 workstring[k] = argv[i][j];
3324 workstring[k] =
'\0';
3325 b->minangle = (
REAL) strtod(workstring, (
char **) NULL);
3330 if (argv[i][j] ==
'a') {
3332 if (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
3333 (argv[i][j + 1] ==
'.')) {
3336 while (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
3337 (argv[i][j + 1] ==
'.')) {
3339 workstring[k] = argv[i][j];
3342 workstring[k] =
'\0';
3343 b->maxarea = (
REAL) strtod(workstring, (
char **) NULL);
3344 if (
b->maxarea <= 0.0) {
3345 printf(
"Error: Maximum area must be greater than zero.\n");
3352 if (argv[i][j] ==
'u') {
3357 if (argv[i][j] ==
'A') {
3358 b->regionattrib = 1;
3360 if (argv[i][j] ==
'c') {
3363 if (argv[i][j] ==
'w') {
3366 if (argv[i][j] ==
'W') {
3369 if (argv[i][j] ==
'j') {
3372 if (argv[i][j] ==
'z') {
3375 if (argv[i][j] ==
'e') {
3378 if (argv[i][j] ==
'v') {
3381 if (argv[i][j] ==
'n') {
3384 if (argv[i][j] ==
'g') {
3387 if (argv[i][j] ==
'B') {
3390 if (argv[i][j] ==
'P') {
3391 b->nopolywritten = 1;
3393 if (argv[i][j] ==
'N') {
3394 b->nonodewritten = 1;
3396 if (argv[i][j] ==
'E') {
3397 b->noelewritten = 1;
3400 if (argv[i][j] ==
'I') {
3401 b->noiterationnum = 1;
3404 if (argv[i][j] ==
'O') {
3407 if (argv[i][j] ==
'X') {
3410 if (argv[i][j] ==
'o') {
3411 if (argv[i][j + 1] ==
'2') {
3417 if (argv[i][j] ==
'Y') {
3420 if (argv[i][j] ==
'S') {
3422 while ((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) {
3424 b->steiner =
b->steiner * 10 + (
int) (argv[i][j] -
'0');
3429 if (argv[i][j] ==
'i') {
3432 if (argv[i][j] ==
'F') {
3436 if (argv[i][j] ==
'l') {
3441 if (argv[i][j] ==
's') {
3444 if ((argv[i][j] ==
'D') || (argv[i][j] ==
'L')) {
3449 if (argv[i][j] ==
'C') {
3453 if (argv[i][j] ==
'Q') {
3456 if (argv[i][j] ==
'V') {
3460 if ((argv[i][j] ==
'h') || (argv[i][j] ==
'H') ||
3461 (argv[i][j] ==
'?')) {
3474 if (
b->innodefilename[0] ==
'\0') {
3477 if (!strcmp(&
b->innodefilename[strlen(
b->innodefilename) - 5],
".node")) {
3478 b->innodefilename[strlen(
b->innodefilename) - 5] =
'\0';
3480 if (!strcmp(&
b->innodefilename[strlen(
b->innodefilename) - 5],
".poly")) {
3481 b->innodefilename[strlen(
b->innodefilename) - 5] =
'\0';
3485 if (!strcmp(&
b->innodefilename[strlen(
b->innodefilename) - 4],
".ele")) {
3486 b->innodefilename[strlen(
b->innodefilename) - 4] =
'\0';
3489 if (!strcmp(&
b->innodefilename[strlen(
b->innodefilename) - 5],
".area")) {
3490 b->innodefilename[strlen(
b->innodefilename) - 5] =
'\0';
3497 b->usesegments =
b->poly ||
b->refine ||
b->quality ||
b->convex;
3498 b->goodangle = cos(
b->minangle *
PI / 180.0);
3499 if (
b->goodangle == 1.0) {
3500 b->offconstant = 0.0;
3502 b->offconstant = 0.475 * sqrt((1.0 +
b->goodangle) / (1.0 -
b->goodangle));
3504 b->goodangle *=
b->goodangle;
3505 if (
b->refine &&
b->noiterationnum) {
3507 "Error: You cannot use the -I switch when refining a triangulation.\n");
3512 if (!
b->refine && !
b->poly) {
3517 if (
b->refine || !
b->poly) {
3518 b->regionattrib = 0;
3522 if (
b->weighted && (
b->poly ||
b->quality)) {
3525 printf(
"Warning: weighted triangulations (-w, -W) are incompatible\n");
3526 printf(
" with PSLGs (-p) and meshing (-q, -a, -u). Weights ignored.\n"
3530 if (
b->jettison &&
b->nonodewritten && !
b->quiet) {
3531 printf(
"Warning: -j and -N switches are somewhat incompatible.\n");
3532 printf(
" If any vertices are jettisoned, you will need the output\n");
3533 printf(
" .node file to reconstruct the new node indices.");
3537 strcpy(
b->inpolyfilename,
b->innodefilename);
3538 strcpy(
b->inelefilename,
b->innodefilename);
3539 strcpy(
b->areafilename,
b->innodefilename);
3541 strcpy(workstring,
b->innodefilename);
3543 while (workstring[j] !=
'\0') {
3544 if ((workstring[j] ==
'.') && (workstring[j + 1] !=
'\0')) {
3550 if (increment > 0) {
3553 if ((workstring[j] >=
'0') && (workstring[j] <=
'9')) {
3554 meshnumber = meshnumber * 10 + (
int) (workstring[j] -
'0');
3559 }
while (workstring[j] !=
'\0');
3561 if (
b->noiterationnum) {
3562 strcpy(
b->outnodefilename,
b->innodefilename);
3563 strcpy(
b->outelefilename,
b->innodefilename);
3564 strcpy(
b->edgefilename,
b->innodefilename);
3565 strcpy(
b->vnodefilename,
b->innodefilename);
3566 strcpy(
b->vedgefilename,
b->innodefilename);
3567 strcpy(
b->neighborfilename,
b->innodefilename);
3568 strcpy(
b->offfilename,
b->innodefilename);
3569 strcat(
b->outnodefilename,
".node");
3570 strcat(
b->outelefilename,
".ele");
3571 strcat(
b->edgefilename,
".edge");
3572 strcat(
b->vnodefilename,
".v.node");
3573 strcat(
b->vedgefilename,
".v.edge");
3574 strcat(
b->neighborfilename,
".neigh");
3575 strcat(
b->offfilename,
".off");
3576 }
else if (increment == 0) {
3577 strcpy(
b->outnodefilename,
b->innodefilename);
3578 strcpy(
b->outpolyfilename,
b->innodefilename);
3579 strcpy(
b->outelefilename,
b->innodefilename);
3580 strcpy(
b->edgefilename,
b->innodefilename);
3581 strcpy(
b->vnodefilename,
b->innodefilename);
3582 strcpy(
b->vedgefilename,
b->innodefilename);
3583 strcpy(
b->neighborfilename,
b->innodefilename);
3584 strcpy(
b->offfilename,
b->innodefilename);
3585 strcat(
b->outnodefilename,
".1.node");
3586 strcat(
b->outpolyfilename,
".1.poly");
3587 strcat(
b->outelefilename,
".1.ele");
3588 strcat(
b->edgefilename,
".1.edge");
3589 strcat(
b->vnodefilename,
".1.v.node");
3590 strcat(
b->vedgefilename,
".1.v.edge");
3591 strcat(
b->neighborfilename,
".1.neigh");
3592 strcat(
b->offfilename,
".1.off");
3594 workstring[increment] =
'%';
3595 workstring[increment + 1] =
'd';
3596 workstring[increment + 2] =
'\0';
3597 sprintf(
b->outnodefilename, workstring, meshnumber + 1);
3598 strcpy(
b->outpolyfilename,
b->outnodefilename);
3599 strcpy(
b->outelefilename,
b->outnodefilename);
3600 strcpy(
b->edgefilename,
b->outnodefilename);
3601 strcpy(
b->vnodefilename,
b->outnodefilename);
3602 strcpy(
b->vedgefilename,
b->outnodefilename);
3603 strcpy(
b->neighborfilename,
b->outnodefilename);
3604 strcpy(
b->offfilename,
b->outnodefilename);
3605 strcat(
b->outnodefilename,
".node");
3606 strcat(
b->outpolyfilename,
".poly");
3607 strcat(
b->outelefilename,
".ele");
3608 strcat(
b->edgefilename,
".edge");
3609 strcat(
b->vnodefilename,
".v.node");
3610 strcat(
b->vedgefilename,
".v.edge");
3611 strcat(
b->neighborfilename,
".neigh");
3612 strcat(
b->offfilename,
".off");
3614 strcat(
b->innodefilename,
".node");
3615 strcat(
b->inpolyfilename,
".poly");
3616 strcat(
b->inelefilename,
".ele");
3617 strcat(
b->areafilename,
".area");
3640#ifdef ANSI_DECLARATORS
3650 struct otri printtri;
3651 struct osub printsh;
3654 printf(
"triangle x%zx with orientation %d:\n", (
size_t) t->
tri,
3657 if (printtri.
tri ==
m->dummytri) {
3658 printf(
" [0] = Outer space\n");
3660 printf(
" [0] = x%zx %d\n", (
size_t) printtri.
tri,
3664 if (printtri.
tri ==
m->dummytri) {
3665 printf(
" [1] = Outer space\n");
3667 printf(
" [1] = x%zx %d\n", (
size_t) printtri.
tri,
3671 if (printtri.
tri ==
m->dummytri) {
3672 printf(
" [2] = Outer space\n");
3674 printf(
" [2] = x%zx %d\n", (
size_t) printtri.
tri,
3678 org(*t, printvertex);
3679 if (printvertex == (
vertex) NULL)
3680 printf(
" Origin[%d] = NULL\n", (t->
orient + 1) % 3 + 3);
3682 printf(
" Origin[%d] = x%zx (%.12g, %.12g)\n",
3683 (t->
orient + 1) % 3 + 3, (
size_t) printvertex,
3684 printvertex[0], printvertex[1]);
3685 dest(*t, printvertex);
3686 if (printvertex == (
vertex) NULL)
3687 printf(
" Dest [%d] = NULL\n", (t->
orient + 2) % 3 + 3);
3689 printf(
" Dest [%d] = x%zx (%.12g, %.12g)\n",
3690 (t->
orient + 2) % 3 + 3, (
size_t) printvertex,
3691 printvertex[0], printvertex[1]);
3692 apex(*t, printvertex);
3693 if (printvertex == (
vertex) NULL)
3694 printf(
" Apex [%d] = NULL\n", t->
orient + 3);
3696 printf(
" Apex [%d] = x%zx (%.12g, %.12g)\n",
3697 t->
orient + 3, (
size_t) printvertex,
3698 printvertex[0], printvertex[1]);
3700 if (
b->usesegments) {
3702 if (printsh.
ss !=
m->dummysub) {
3703 printf(
" [6] = x%zx %d\n", (
size_t) printsh.
ss,
3707 if (printsh.
ss !=
m->dummysub) {
3708 printf(
" [7] = x%zx %d\n", (
size_t) printsh.
ss,
3712 if (printsh.
ss !=
m->dummysub) {
3713 printf(
" [8] = x%zx %d\n", (
size_t) printsh.
ss,
3719 printf(
" Area constraint: %.4g\n",
areabound(*t));
3734#ifdef ANSI_DECLARATORS
3744 struct osub printsh;
3745 struct otri printtri;
3750 printf(
"subsegment x%zx with orientation %d and mark %d:\n",
3753 if (printsh.
ss ==
m->dummysub) {
3754 printf(
" [0] = No subsegment\n");
3756 printf(
" [0] = x%zx %d\n", (
size_t) printsh.
ss,
3760 if (printsh.
ss ==
m->dummysub) {
3761 printf(
" [1] = No subsegment\n");
3763 printf(
" [1] = x%zx %d\n", (
size_t) printsh.
ss,
3767 sorg(*s, printvertex);
3768 if (printvertex == (
vertex) NULL)
3769 printf(
" Origin[%d] = NULL\n", 2 + s->
ssorient);
3771 printf(
" Origin[%d] = x%zx (%.12g, %.12g)\n",
3772 2 + s->
ssorient, (
size_t) printvertex,
3773 printvertex[0], printvertex[1]);
3774 sdest(*s, printvertex);
3775 if (printvertex == (
vertex) NULL)
3776 printf(
" Dest [%d] = NULL\n", 3 - s->
ssorient);
3778 printf(
" Dest [%d] = x%zx (%.12g, %.12g)\n",
3779 3 - s->
ssorient, (
size_t) printvertex,
3780 printvertex[0], printvertex[1]);
3783 if (printtri.
tri ==
m->dummytri) {
3784 printf(
" [6] = Outer space\n");
3786 printf(
" [6] = x%zx %d\n", (
size_t) printtri.
tri,
3790 if (printtri.
tri ==
m->dummytri) {
3791 printf(
" [7] = Outer space\n");
3793 printf(
" [7] = x%zx %d\n", (
size_t) printtri.
tri,
3798 if (printvertex == (
vertex) NULL)
3799 printf(
" Segment origin[%d] = NULL\n", 4 + s->
ssorient);
3801 printf(
" Segment origin[%d] = x%zx (%.12g, %.12g)\n",
3802 4 + s->
ssorient, (
size_t) printvertex,
3803 printvertex[0], printvertex[1]);
3805 if (printvertex == (
vertex) NULL)
3806 printf(
" Segment dest [%d] = NULL\n", 5 - s->
ssorient);
3808 printf(
" Segment dest [%d] = x%zx (%.12g, %.12g)\n",
3809 5 - s->
ssorient, (
size_t) printvertex,
3810 printvertex[0], printvertex[1]);
3830#ifdef ANSI_DECLARATORS
3864#ifdef ANSI_DECLARATORS
3880 alignptr = (uintptr_t) (pool->
nowblock + 1);
3910#ifdef ANSI_DECLARATORS
3935 if (firstitemcount == 0) {
3958#ifdef ANSI_DECLARATORS
3979#ifdef ANSI_DECLARATORS
4003 (
int)
sizeof(
VOID *) +
4007 *newblock = (
VOID *) NULL;
4014 alignptr = (uintptr_t) (pool->
nowblock + 1);
4042#ifdef ANSI_DECLARATORS
4065#ifdef ANSI_DECLARATORS
4078 alignptr = (uintptr_t) (pool->
pathblock + 1);
4101#ifdef ANSI_DECLARATORS
4114 return (
VOID *) NULL;
4122 alignptr = (uintptr_t) (pool->
pathblock + 1);
4166#ifdef ANSI_DECLARATORS
4182 m->triangles.alignbytes);
4184 alignptr = (uintptr_t)
m->dummytribase;
4186 (alignptr + (uintptr_t)
m->triangles.alignbytes -
4187 (alignptr % (uintptr_t)
m->triangles.alignbytes));
4200 if (
b->usesegments) {
4205 m->subsegs.alignbytes);
4207 alignptr = (uintptr_t)
m->dummysubbase;
4209 (alignptr + (uintptr_t)
m->subsegs.alignbytes -
4210 (alignptr % (uintptr_t)
m->subsegs.alignbytes));
4215 m->dummysub[0] = (
subseg)
m->dummysub;
4216 m->dummysub[1] = (
subseg)
m->dummysub;
4218 m->dummysub[2] = (
subseg) NULL;
4219 m->dummysub[3] = (
subseg) NULL;
4220 m->dummysub[4] = (
subseg) NULL;
4221 m->dummysub[5] = (
subseg) NULL;
4223 m->dummysub[6] = (
subseg)
m->dummytri;
4224 m->dummysub[7] = (
subseg)
m->dummytri;
4226 * (
int *) (
m->dummysub + 8) = 0;
4246#ifdef ANSI_DECLARATORS
4260 m->vertexmarkindex = ((
m->mesh_dim +
m->nextras) *
sizeof(
REAL) +
4263 vertexsize = (
m->vertexmarkindex + 2) *
sizeof(
int);
4267 m->vertex2triindex = (vertexsize +
sizeof(
triangle) - 1) /
4269 vertexsize = (
m->vertex2triindex + 1) *
sizeof(
triangle);
4289#ifdef ANSI_DECLARATORS
4304 m->highorderindex = 6 + (
b->usesegments * 3);
4306 trisize = ((
b->order + 1) * (
b->order + 2) / 2 + (
m->highorderindex - 3)) *
4310 m->elemattribindex = (trisize +
sizeof(
REAL) - 1) /
sizeof(
REAL);
4314 m->areaboundindex =
m->elemattribindex +
m->eextras +
b->regionattrib;
4318 trisize = (
m->areaboundindex + 1) *
sizeof(
REAL);
4319 }
else if (
m->eextras +
b->regionattrib > 0) {
4320 trisize =
m->areaboundindex *
sizeof(
REAL);
4326 if ((
b->voronoi ||
b->neighbors) &&
4328 trisize = 6 *
sizeof(
triangle) +
sizeof(
int);
4333 (2 *
m->invertices - 2) >
TRIPERBLOCK ? (2 *
m->invertices - 2) :
4336 if (
b->usesegments) {
4343 dummyinit(
m,
b,
m->triangles.itembytes,
m->subsegs.itembytes);
4356#ifdef ANSI_DECLARATORS
4377#ifdef ANSI_DECLARATORS
4389 if (newtriangle == (
triangle *) NULL) {
4392 }
while (
deadtri(newtriangle));
4402#ifdef ANSI_DECLARATORS
4423#ifdef ANSI_DECLARATORS
4435 if (newsubseg == (
subseg *) NULL) {
4448#ifdef ANSI_DECLARATORS
4469#ifdef ANSI_DECLARATORS
4481 if (newvertex == (
vertex) NULL) {
4497#ifdef ANSI_DECLARATORS
4500void badsubsegdealloc(
m, dyingseg)
4522#ifdef ANSI_DECLARATORS
4534 if (newseg == (
struct badsubseg *) NULL) {
4555#ifdef ANSI_DECLARATORS
4570 getblock =
m->vertices.firstblock;
4571 current =
b->firstnumber;
4574 if (current +
m->vertices.itemsfirstblock <= number) {
4575 getblock = (
VOID **) *getblock;
4576 current +=
m->vertices.itemsfirstblock;
4577 while (current +
m->vertices.itemsperblock <= number) {
4578 getblock = (
VOID **) *getblock;
4579 current +=
m->vertices.itemsperblock;
4584 alignptr = (uintptr_t) (getblock + 1);
4585 foundvertex = (
char *) (alignptr + (uintptr_t)
m->vertices.alignbytes -
4586 (alignptr % (uintptr_t)
m->vertices.alignbytes));
4587 return (
vertex) (foundvertex +
m->vertices.itembytes * (number - current));
4596#ifdef ANSI_DECLARATORS
4607 if (
b->usesegments) {
4615 if ((
b->minangle > 0.0) ||
b->vararea ||
b->fixedarea ||
b->usertest) {
4637#ifdef ANSI_DECLARATORS
4643struct otri *newotri;
4658 if (
b->usesegments) {
4665 for (i = 0; i <
m->eextras; i++) {
4681#ifdef ANSI_DECLARATORS
4686struct osub *newsubseg;
4693 newsubseg->
ss[0] = (
subseg)
m->dummysub;
4694 newsubseg->
ss[1] = (
subseg)
m->dummysub;
4701 newsubseg->
ss[6] = (
subseg)
m->dummytri;
4702 newsubseg->
ss[7] = (
subseg)
m->dummytri;
4730#define Absolute(a) ((a) >= 0.0 ? (a) : -(a))
4746#define Fast_Two_Sum_Tail(a, b, x, y) \
4750#define Fast_Two_Sum(a, b, x, y) \
4751 x = (REAL) (a + b); \
4752 Fast_Two_Sum_Tail(a, b, x, y)
4754#define Two_Sum_Tail(a, b, x, y) \
4755 bvirt = (REAL) (x - a); \
4756 avirt = x - bvirt; \
4757 bround = b - bvirt; \
4758 around = a - avirt; \
4761#define Two_Sum(a, b, x, y) \
4762 x = (REAL) (a + b); \
4763 Two_Sum_Tail(a, b, x, y)
4765#define Two_Diff_Tail(a, b, x, y) \
4766 bvirt = (REAL) (a - x); \
4767 avirt = x + bvirt; \
4768 bround = bvirt - b; \
4769 around = a - avirt; \
4772#define Two_Diff(a, b, x, y) \
4773 x = (REAL) (a - b); \
4774 Two_Diff_Tail(a, b, x, y)
4776#define Split(a, ahi, aLo) \
4777 c = (REAL) (splitter * a); \
4778 abig = (REAL) (c - a); \
4782#define Two_Product_Tail(a, b, x, y) \
4783 Split(a, ahi, aLo); \
4784 Split(b, bhi, blo); \
4785 err1 = x - (ahi * bhi); \
4786 err2 = err1 - (aLo * bhi); \
4787 err3 = err2 - (ahi * blo); \
4788 y = (aLo * blo) - err3
4790#define Two_Product(a, b, x, y) \
4791 x = (REAL) (a * b); \
4792 Two_Product_Tail(a, b, x, y)
4797#define Two_Product_Presplit(a, b, bhi, blo, x, y) \
4798 x = (REAL) (a * b); \
4799 Split(a, ahi, aLo); \
4800 err1 = x - (ahi * bhi); \
4801 err2 = err1 - (aLo * bhi); \
4802 err3 = err2 - (ahi * blo); \
4803 y = (aLo * blo) - err3
4807#define Square_Tail(a, x, y) \
4808 Split(a, ahi, aLo); \
4809 err1 = x - (ahi * ahi); \
4810 err3 = err1 - ((ahi + ahi) * aLo); \
4811 y = (aLo * aLo) - err3
4813#define Square(a, x, y) \
4814 x = (REAL) (a * a); \
4815 Square_Tail(a, x, y)
4820#define Two_One_Sum(a1, a0, b, x2, x1, x0) \
4821 Two_Sum(a0, b , _i, x0); \
4822 Two_Sum(a1, _i, x2, x1)
4824#define Two_One_Diff(a1, a0, b, x2, x1, x0) \
4825 Two_Diff(a0, b , _i, x0); \
4826 Two_Sum( a1, _i, x2, x1)
4828#define Two_Two_Sum(a1, a0, b1, b0, x3, x2, x1, x0) \
4829 Two_One_Sum(a1, a0, b0, _j, _0, x0); \
4830 Two_One_Sum(_j, _0, b1, x3, x2, x1)
4832#define Two_Two_Diff(a1, a0, b1, b0, x3, x2, x1, x0) \
4833 Two_One_Diff(a1, a0, b0, _j, _0, x0); \
4834 Two_One_Diff(_j, _0, b1, x3, x2, x1)
4838#define Two_One_Product(a1, a0, b, x3, x2, x1, x0) \
4839 Split(b, bhi, blo); \
4840 Two_Product_Presplit(a0, b, bhi, blo, _i, x0); \
4841 Two_Product_Presplit(a1, b, bhi, blo, _j, _0); \
4842 Two_Sum(_i, _0, _k, x1); \
4843 Fast_Two_Sum(_j, _k, x3, x2)
4867 REAL check, lastcheck;
4875 _control87(_PC_24, _MCW_PC);
4877 _control87(_PC_53, _MCW_PC);
4883 fenv.__control_word = 4210;
4886 fenv.__control_word = 4722;
4906 every_other = !every_other;
4908 }
while ((check != 1.0) && (check != lastcheck));
4937#ifdef ANSI_DECLARATORS
4953 REAL avirt, bround, around;
4954 int eindex, findex, hindex;
4959 eindex = findex = 0;
4960 if ((fnow > enow) == (fnow > -enow)) {
4968 if ((eindex < elen) && (findex < flen)) {
4969 if ((fnow > enow) == (fnow > -enow)) {
4980 while ((eindex < elen) && (findex < flen)) {
4981 if ((fnow > enow) == (fnow > -enow)) {
4994 while (eindex < elen) {
5002 while (findex < flen) {
5010 if ((Q != 0.0) || (hindex == 0)) {
5031#ifdef ANSI_DECLARATORS
5049 REAL avirt, bround, around;
5052 REAL ahi, aLo, bhi, blo;
5053 REAL err1, err2, err3;
5061 for (eindex = 1; eindex < elen; eindex++) {
5073 if ((Q != 0.0) || (hindex == 0)) {
5087#ifdef ANSI_DECLARATORS
5100 for (eindex = 1; eindex < elen; eindex++) {
5126#ifdef ANSI_DECLARATORS
5138 REAL acxtail, acytail, bcxtail, bcytail;
5140 REAL detlefttail, detrighttail;
5142 REAL B[4], C1[8], C2[12], D[16];
5144 int C1length, C2length, Dlength;
5151 REAL avirt, bround, around;
5154 REAL ahi, aLo, bhi, blo;
5155 REAL err1, err2, err3;
5159 acx = (
REAL) (pa[0] - pc[0]);
5160 bcx = (
REAL) (pb[0] - pc[0]);
5161 acy = (
REAL) (pa[1] - pc[1]);
5162 bcy = (
REAL) (pb[1] - pc[1]);
5167 Two_Two_Diff(detleft, detlefttail, detright, detrighttail,
5168 B3, B[2], B[1], B[0]);
5173 if ((det >= errbound) || (-det >= errbound)) {
5182 if ((acxtail == 0.0) && (acytail == 0.0)
5183 && (bcxtail == 0.0) && (bcytail == 0.0)) {
5188 det += (acx * bcytail + bcy * acxtail)
5189 - (acy * bcxtail + bcx * acytail);
5190 if ((det >= errbound) || (-det >= errbound)) {
5212 return(D[Dlength - 1]);
5215#ifdef ANSI_DECLARATORS
5228 REAL detleft, detright, det;
5229 REAL detsum, errbound;
5231 m->counterclockcount++;
5233 detleft = (pa[0] - pc[0]) * (pb[1] - pc[1]);
5234 detright = (pa[1] - pc[1]) * (pb[0] - pc[0]);
5235 det = detleft - detright;
5241 if (detleft > 0.0) {
5242 if (detright <= 0.0) {
5245 detsum = detleft + detright;
5247 }
else if (detleft < 0.0) {
5248 if (detright >= 0.0) {
5251 detsum = -detleft - detright;
5258 if ((det >= errbound) || (-det >= errbound)) {
5284#ifdef ANSI_DECLARATORS
5299 INEXACT REAL bdxcdy1, cdxbdy1, cdxady1, adxcdy1, adxbdy1, bdxady1;
5300 REAL bdxcdy0, cdxbdy0, cdxady0, adxcdy0, adxbdy0, bdxady0;
5301 REAL bc[4], ca[4], ab[4];
5303 REAL axbc[8], axxbc[16], aybc[8], ayybc[16], adet[32];
5304 int axbcLen, axxbcLen, aybcLen, ayybcLen, alen;
5305 REAL bxca[8], bxxca[16], byca[8], byyca[16], bdet[32];
5306 int bxcalen, bxxcalen, bycalen, byycalen, blen;
5307 REAL cxab[8], cxxab[16], cyab[8], cyyab[16], cdet[32];
5308 int cxablen, cxxablen, cyablen, cyyablen, cLen;
5311 REAL fin1[1152], fin2[1152];
5312 REAL *finnow, *finother, *finswap;
5315 REAL adxtail, bdxtail, cdxtail, adytail, bdytail, cdytail;
5316 INEXACT REAL adxadx1, adyady1, bdxbdx1, bdybdy1, cdxcdx1, cdycdy1;
5317 REAL adxadx0, adyady0, bdxbdx0, bdybdy0, cdxcdx0, cdycdy0;
5318 REAL aa[4], bb[4], cc[4];
5324 REAL temp8[8], temp16a[16], temp16b[16], temp16c[16];
5325 REAL temp32a[32], temp32b[32], temp48[48], temp64[64];
5326 int temp8len, temp16alen, temp16blen, temp16cLen;
5327 int temp32alen, temp32blen, temp48len, temp64len;
5328 REAL axtbb[8], axtcc[8], aytbb[8], aytcc[8];
5329 int axtbblen, axtccLen, aytbblen, aytccLen;
5330 REAL bxtaa[8], bxtcc[8], bytaa[8], bytcc[8];
5331 int bxtaalen, bxtccLen, bytaalen, bytccLen;
5332 REAL cxtaa[8], cxtbb[8], cytaa[8], cytbb[8];
5333 int cxtaalen, cxtbblen, cytaalen, cytbblen;
5334 REAL axtbc[8], aytbc[8], bxtca[8], bytca[8], cxtab[8], cytab[8];
5335 int axtbcLen, aytbcLen, bxtcalen, bytcalen, cxtablen, cytablen;
5336 REAL axtbct[16], aytbct[16], bxtcat[16], bytcat[16], cxtabt[16], cytabt[16];
5337 int axtbctlen, aytbctlen, bxtcatlen, bytcatlen, cxtabtlen, cytabtlen;
5338 REAL axtbctt[8], aytbctt[8], bxtcatt[8];
5339 REAL bytcatt[8], cxtabtt[8], cytabtt[8];
5340 int axtbcttlen, aytbcttlen, bxtcattlen, bytcattlen, cxtabttlen, cytabttlen;
5341 REAL abt[8], bct[8], cat[8];
5342 int abtlen, bctlen, catlen;
5343 REAL abtt[4], bctt[4], catt[4];
5344 int abttlen, bcttlen, cattlen;
5349 REAL avirt, bround, around;
5352 REAL ahi, aLo, bhi, blo;
5353 REAL err1, err2, err3;
5357 adx = (
REAL) (pa[0] - pd[0]);
5358 bdx = (
REAL) (pb[0] - pd[0]);
5359 cdx = (
REAL) (pc[0] - pd[0]);
5360 ady = (
REAL) (pa[1] - pd[1]);
5361 bdy = (
REAL) (pb[1] - pd[1]);
5362 cdy = (
REAL) (pc[1] - pd[1]);
5366 Two_Two_Diff(bdxcdy1, bdxcdy0, cdxbdy1, cdxbdy0, bc3, bc[2], bc[1], bc[0]);
5376 Two_Two_Diff(cdxady1, cdxady0, adxcdy1, adxcdy0, ca3, ca[2], ca[1], ca[0]);
5386 Two_Two_Diff(adxbdy1, adxbdy0, bdxady1, bdxady0, ab3, ab[2], ab[1], ab[0]);
5399 if ((det >= errbound) || (-det >= errbound)) {
5409 if ((adxtail == 0.0) && (bdxtail == 0.0) && (cdxtail == 0.0)
5410 && (adytail == 0.0) && (bdytail == 0.0) && (cdytail == 0.0)) {
5415 det += ((adx * adx + ady * ady) * ((bdx * cdytail + cdy * bdxtail)
5416 - (bdy * cdxtail + cdx * bdytail))
5417 + 2.0 * (adx * adxtail + ady * adytail) * (bdx * cdy - bdy * cdx))
5418 + ((bdx * bdx + bdy * bdy) * ((cdx * adytail + ady * cdxtail)
5419 - (cdy * adxtail + adx * cdytail))
5420 + 2.0 * (bdx * bdxtail + bdy * bdytail) * (cdx * ady - cdy * adx))
5421 + ((cdx * cdx + cdy * cdy) * ((adx * bdytail + bdy * adxtail)
5422 - (ady * bdxtail + bdx * adytail))
5423 + 2.0 * (cdx * cdxtail + cdy * cdytail) * (adx * bdy - ady * bdx));
5424 if ((det >= errbound) || (-det >= errbound)) {
5431 if ((bdxtail != 0.0) || (bdytail != 0.0)
5432 || (cdxtail != 0.0) || (cdytail != 0.0)) {
5433 Square(adx, adxadx1, adxadx0);
5434 Square(ady, adyady1, adyady0);
5435 Two_Two_Sum(adxadx1, adxadx0, adyady1, adyady0, aa3, aa[2], aa[1], aa[0]);
5438 if ((cdxtail != 0.0) || (cdytail != 0.0)
5439 || (adxtail != 0.0) || (adytail != 0.0)) {
5440 Square(bdx, bdxbdx1, bdxbdx0);
5441 Square(bdy, bdybdy1, bdybdy0);
5442 Two_Two_Sum(bdxbdx1, bdxbdx0, bdybdy1, bdybdy0, bb3, bb[2], bb[1], bb[0]);
5445 if ((adxtail != 0.0) || (adytail != 0.0)
5446 || (bdxtail != 0.0) || (bdytail != 0.0)) {
5447 Square(cdx, cdxcdx1, cdxcdx0);
5448 Square(cdy, cdycdy1, cdycdy0);
5449 Two_Two_Sum(cdxcdx1, cdxcdx0, cdycdy1, cdycdy0, cc3, cc[2], cc[1], cc[0]);
5453 if (adxtail != 0.0) {
5465 temp16blen, temp16b, temp32a);
5467 temp32alen, temp32a, temp48);
5470 finswap = finnow; finnow = finother; finother = finswap;
5472 if (adytail != 0.0) {
5484 temp16blen, temp16b, temp32a);
5486 temp32alen, temp32a, temp48);
5489 finswap = finnow; finnow = finother; finother = finswap;
5491 if (bdxtail != 0.0) {
5503 temp16blen, temp16b, temp32a);
5505 temp32alen, temp32a, temp48);
5508 finswap = finnow; finnow = finother; finother = finswap;
5510 if (bdytail != 0.0) {
5522 temp16blen, temp16b, temp32a);
5524 temp32alen, temp32a, temp48);
5527 finswap = finnow; finnow = finother; finother = finswap;
5529 if (cdxtail != 0.0) {
5541 temp16blen, temp16b, temp32a);
5543 temp32alen, temp32a, temp48);
5546 finswap = finnow; finnow = finother; finother = finswap;
5548 if (cdytail != 0.0) {
5560 temp16blen, temp16b, temp32a);
5562 temp32alen, temp32a, temp48);
5565 finswap = finnow; finnow = finother; finother = finswap;
5568 if ((adxtail != 0.0) || (adytail != 0.0)) {
5569 if ((bdxtail != 0.0) || (bdytail != 0.0)
5570 || (cdxtail != 0.0) || (cdytail != 0.0)) {
5573 Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]);
5585 Two_Two_Diff(ti1, ti0, tj1, tj0, bctt3, bctt[2], bctt[1], bctt[0]);
5595 if (adxtail != 0.0) {
5601 temp32alen, temp32a, temp48);
5604 finswap = finnow; finnow = finother; finother = finswap;
5605 if (bdytail != 0.0) {
5611 finswap = finnow; finnow = finother; finother = finswap;
5613 if (cdytail != 0.0) {
5619 finswap = finnow; finnow = finother; finother = finswap;
5630 temp16blen, temp16b, temp32b);
5632 temp32blen, temp32b, temp64);
5635 finswap = finnow; finnow = finother; finother = finswap;
5637 if (adytail != 0.0) {
5643 temp32alen, temp32a, temp48);
5646 finswap = finnow; finnow = finother; finother = finswap;
5657 temp16blen, temp16b, temp32b);
5659 temp32blen, temp32b, temp64);
5662 finswap = finnow; finnow = finother; finother = finswap;
5665 if ((bdxtail != 0.0) || (bdytail != 0.0)) {
5666 if ((cdxtail != 0.0) || (cdytail != 0.0)
5667 || (adxtail != 0.0) || (adytail != 0.0)) {
5670 Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]);
5682 Two_Two_Diff(ti1, ti0, tj1, tj0, catt3, catt[2], catt[1], catt[0]);
5692 if (bdxtail != 0.0) {
5698 temp32alen, temp32a, temp48);
5701 finswap = finnow; finnow = finother; finother = finswap;
5702 if (cdytail != 0.0) {
5708 finswap = finnow; finnow = finother; finother = finswap;
5710 if (adytail != 0.0) {
5716 finswap = finnow; finnow = finother; finother = finswap;
5727 temp16blen, temp16b, temp32b);
5729 temp32blen, temp32b, temp64);
5732 finswap = finnow; finnow = finother; finother = finswap;
5734 if (bdytail != 0.0) {
5740 temp32alen, temp32a, temp48);
5743 finswap = finnow; finnow = finother; finother = finswap;
5754 temp16blen, temp16b, temp32b);
5756 temp32blen, temp32b, temp64);
5759 finswap = finnow; finnow = finother; finother = finswap;
5762 if ((cdxtail != 0.0) || (cdytail != 0.0)) {
5763 if ((adxtail != 0.0) || (adytail != 0.0)
5764 || (bdxtail != 0.0) || (bdytail != 0.0)) {
5767 Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]);
5779 Two_Two_Diff(ti1, ti0, tj1, tj0, abtt3, abtt[2], abtt[1], abtt[0]);
5789 if (cdxtail != 0.0) {
5795 temp32alen, temp32a, temp48);
5798 finswap = finnow; finnow = finother; finother = finswap;
5799 if (adytail != 0.0) {
5805 finswap = finnow; finnow = finother; finother = finswap;
5807 if (bdytail != 0.0) {
5813 finswap = finnow; finnow = finother; finother = finswap;
5824 temp16blen, temp16b, temp32b);
5826 temp32blen, temp32b, temp64);
5829 finswap = finnow; finnow = finother; finother = finswap;
5831 if (cdytail != 0.0) {
5837 temp32alen, temp32a, temp48);
5840 finswap = finnow; finnow = finother; finother = finswap;
5851 temp16blen, temp16b, temp32b);
5853 temp32blen, temp32b, temp64);
5856 finswap = finnow; finnow = finother; finother = finswap;
5860 return finnow[finlength - 1];
5863#ifdef ANSI_DECLARATORS
5877 REAL adx, bdx, cdx, ady, bdy, cdy;
5878 REAL bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady;
5879 REAL alift, blift, clift;
5881 REAL permanent, errbound;
5885 adx = pa[0] - pd[0];
5886 bdx = pb[0] - pd[0];
5887 cdx = pc[0] - pd[0];
5888 ady = pa[1] - pd[1];
5889 bdy = pb[1] - pd[1];
5890 cdy = pc[1] - pd[1];
5894 alift = adx * adx + ady * ady;
5898 blift = bdx * bdx + bdy * bdy;
5902 clift = cdx * cdx + cdy * cdy;
5904 det = alift * (bdxcdy - cdxbdy)
5905 + blift * (cdxady - adxcdy)
5906 + clift * (adxbdy - bdxady);
5916 if ((det > errbound) || (-det > errbound)) {
5945#ifdef ANSI_DECLARATORS
5951 aheight, bheight, cheight, dheight, permanent)
5964 INEXACT REAL adx, bdx, cdx, ady, bdy, cdy, adheight, bdheight, cdheight;
5967 INEXACT REAL bdxcdy1, cdxbdy1, cdxady1, adxcdy1, adxbdy1, bdxady1;
5968 REAL bdxcdy0, cdxbdy0, cdxady0, adxcdy0, adxbdy0, bdxady0;
5969 REAL bc[4], ca[4], ab[4];
5971 REAL adet[8], bdet[8], cdet[8];
5972 int alen, blen, cLen;
5975 REAL *finnow, *finother, *finswap;
5976 REAL fin1[192], fin2[192];
5979 REAL adxtail, bdxtail, cdxtail;
5980 REAL adytail, bdytail, cdytail;
5981 REAL adheighttail, bdheighttail, cdheighttail;
5985 REAL at_b[4], at_c[4], bt_c[4], bt_a[4], ct_a[4], ct_b[4];
5986 int at_blen, at_cLen, bt_cLen, bt_alen, ct_alen, ct_blen;
5989 REAL bdxt_cdy0, cdxt_bdy0, cdxt_ady0;
5990 REAL adxt_cdy0, adxt_bdy0, bdxt_ady0;
5993 REAL bdyt_cdx0, cdyt_bdx0, cdyt_adx0;
5994 REAL adyt_cdx0, adyt_bdx0, bdyt_adx0;
5995 REAL bct[8], cat[8], abt[8];
5996 int bctlen, catlen, abtlen;
5999 REAL bdxt_cdyt0, cdxt_bdyt0, cdxt_adyt0;
6000 REAL adxt_cdyt0, adxt_bdyt0, bdxt_adyt0;
6003 int vlength, wlength;
6007 REAL avirt, bround, around;
6010 REAL ahi, aLo, bhi, blo;
6011 REAL err1, err2, err3;
6015 adx = (
REAL) (pa[0] - pd[0]);
6016 bdx = (
REAL) (pb[0] - pd[0]);
6017 cdx = (
REAL) (pc[0] - pd[0]);
6018 ady = (
REAL) (pa[1] - pd[1]);
6019 bdy = (
REAL) (pb[1] - pd[1]);
6020 cdy = (
REAL) (pc[1] - pd[1]);
6021 adheight = (
REAL) (aheight - dheight);
6022 bdheight = (
REAL) (bheight - dheight);
6023 cdheight = (
REAL) (cheight - dheight);
6027 Two_Two_Diff(bdxcdy1, bdxcdy0, cdxbdy1, cdxbdy0, bc3, bc[2], bc[1], bc[0]);
6033 Two_Two_Diff(cdxady1, cdxady0, adxcdy1, adxcdy0, ca3, ca[2], ca[1], ca[0]);
6039 Two_Two_Diff(adxbdy1, adxbdy0, bdxady1, bdxady0, ab3, ab[2], ab[1], ab[0]);
6048 if ((det >= errbound) || (-det >= errbound)) {
6062 if ((adxtail == 0.0) && (bdxtail == 0.0) && (cdxtail == 0.0) &&
6063 (adytail == 0.0) && (bdytail == 0.0) && (cdytail == 0.0) &&
6064 (adheighttail == 0.0) &&
6065 (bdheighttail == 0.0) &&
6066 (cdheighttail == 0.0)) {
6071 det += (adheight * ((bdx * cdytail + cdy * bdxtail) -
6072 (bdy * cdxtail + cdx * bdytail)) +
6073 adheighttail * (bdx * cdy - bdy * cdx)) +
6074 (bdheight * ((cdx * adytail + ady * cdxtail) -
6075 (cdy * adxtail + adx * cdytail)) +
6076 bdheighttail * (cdx * ady - cdy * adx)) +
6077 (cdheight * ((adx * bdytail + bdy * adxtail) -
6078 (ady * bdxtail + bdx * adytail)) +
6079 cdheighttail * (adx * bdy - ady * bdx));
6080 if ((det >= errbound) || (-det >= errbound)) {
6087 if (adxtail == 0.0) {
6088 if (adytail == 0.0) {
6096 at_b[1] = at_blarge;
6099 at_c[1] = at_clarge;
6103 if (adytail == 0.0) {
6105 at_b[1] = at_blarge;
6109 at_c[1] = at_clarge;
6114 Two_Two_Diff(adxt_bdy1, adxt_bdy0, adyt_bdx1, adyt_bdx0,
6115 at_blarge, at_b[2], at_b[1], at_b[0]);
6116 at_b[3] = at_blarge;
6120 Two_Two_Diff(adyt_cdx1, adyt_cdx0, adxt_cdy1, adxt_cdy0,
6121 at_clarge, at_c[2], at_c[1], at_c[0]);
6122 at_c[3] = at_clarge;
6126 if (bdxtail == 0.0) {
6127 if (bdytail == 0.0) {
6135 bt_c[1] = bt_clarge;
6138 bt_a[1] = bt_alarge;
6142 if (bdytail == 0.0) {
6144 bt_c[1] = bt_clarge;
6148 bt_a[1] = bt_alarge;
6153 Two_Two_Diff(bdxt_cdy1, bdxt_cdy0, bdyt_cdx1, bdyt_cdx0,
6154 bt_clarge, bt_c[2], bt_c[1], bt_c[0]);
6155 bt_c[3] = bt_clarge;
6159 Two_Two_Diff(bdyt_adx1, bdyt_adx0, bdxt_ady1, bdxt_ady0,
6160 bt_alarge, bt_a[2], bt_a[1], bt_a[0]);
6161 bt_a[3] = bt_alarge;
6165 if (cdxtail == 0.0) {
6166 if (cdytail == 0.0) {
6174 ct_a[1] = ct_alarge;
6177 ct_b[1] = ct_blarge;
6181 if (cdytail == 0.0) {
6183 ct_a[1] = ct_alarge;
6187 ct_b[1] = ct_blarge;
6192 Two_Two_Diff(cdxt_ady1, cdxt_ady0, cdyt_adx1, cdyt_adx0,
6193 ct_alarge, ct_a[2], ct_a[1], ct_a[0]);
6194 ct_a[3] = ct_alarge;
6198 Two_Two_Diff(cdyt_bdx1, cdyt_bdx0, cdxt_bdy1, cdxt_bdy0,
6199 ct_blarge, ct_b[2], ct_b[1], ct_b[0]);
6200 ct_b[3] = ct_blarge;
6209 finswap = finnow; finnow = finother; finother = finswap;
6215 finswap = finnow; finnow = finother; finother = finswap;
6221 finswap = finnow; finnow = finother; finother = finswap;
6223 if (adheighttail != 0.0) {
6227 finswap = finnow; finnow = finother; finother = finswap;
6229 if (bdheighttail != 0.0) {
6233 finswap = finnow; finnow = finother; finother = finswap;
6235 if (cdheighttail != 0.0) {
6239 finswap = finnow; finnow = finother; finother = finswap;
6242 if (adxtail != 0.0) {
6243 if (bdytail != 0.0) {
6244 Two_Product(adxtail, bdytail, adxt_bdyt1, adxt_bdyt0);
6245 Two_One_Product(adxt_bdyt1, adxt_bdyt0, cdheight, u3, u[2], u[1], u[0]);
6249 finswap = finnow; finnow = finother; finother = finswap;
6250 if (cdheighttail != 0.0) {
6252 u3, u[2], u[1], u[0]);
6256 finswap = finnow; finnow = finother; finother = finswap;
6259 if (cdytail != 0.0) {
6261 Two_Product(negate, cdytail, adxt_cdyt1, adxt_cdyt0);
6262 Two_One_Product(adxt_cdyt1, adxt_cdyt0, bdheight, u3, u[2], u[1], u[0]);
6266 finswap = finnow; finnow = finother; finother = finswap;
6267 if (bdheighttail != 0.0) {
6269 u3, u[2], u[1], u[0]);
6273 finswap = finnow; finnow = finother; finother = finswap;
6277 if (bdxtail != 0.0) {
6278 if (cdytail != 0.0) {
6279 Two_Product(bdxtail, cdytail, bdxt_cdyt1, bdxt_cdyt0);
6280 Two_One_Product(bdxt_cdyt1, bdxt_cdyt0, adheight, u3, u[2], u[1], u[0]);
6284 finswap = finnow; finnow = finother; finother = finswap;
6285 if (adheighttail != 0.0) {
6287 u3, u[2], u[1], u[0]);
6291 finswap = finnow; finnow = finother; finother = finswap;
6294 if (adytail != 0.0) {
6296 Two_Product(negate, adytail, bdxt_adyt1, bdxt_adyt0);
6297 Two_One_Product(bdxt_adyt1, bdxt_adyt0, cdheight, u3, u[2], u[1], u[0]);
6301 finswap = finnow; finnow = finother; finother = finswap;
6302 if (cdheighttail != 0.0) {
6304 u3, u[2], u[1], u[0]);
6308 finswap = finnow; finnow = finother; finother = finswap;
6312 if (cdxtail != 0.0) {
6313 if (adytail != 0.0) {
6314 Two_Product(cdxtail, adytail, cdxt_adyt1, cdxt_adyt0);
6315 Two_One_Product(cdxt_adyt1, cdxt_adyt0, bdheight, u3, u[2], u[1], u[0]);
6319 finswap = finnow; finnow = finother; finother = finswap;
6320 if (bdheighttail != 0.0) {
6322 u3, u[2], u[1], u[0]);
6326 finswap = finnow; finnow = finother; finother = finswap;
6329 if (bdytail != 0.0) {
6331 Two_Product(negate, bdytail, cdxt_bdyt1, cdxt_bdyt0);
6332 Two_One_Product(cdxt_bdyt1, cdxt_bdyt0, adheight, u3, u[2], u[1], u[0]);
6336 finswap = finnow; finnow = finother; finother = finswap;
6337 if (adheighttail != 0.0) {
6339 u3, u[2], u[1], u[0]);
6343 finswap = finnow; finnow = finother; finother = finswap;
6348 if (adheighttail != 0.0) {
6352 finswap = finnow; finnow = finother; finother = finswap;
6354 if (bdheighttail != 0.0) {
6358 finswap = finnow; finnow = finother; finother = finswap;
6360 if (cdheighttail != 0.0) {
6364 finswap = finnow; finnow = finother; finother = finswap;
6367 return finnow[finlength - 1];
6370#ifdef ANSI_DECLARATORS
6375REAL orient3d(
m,
b, pa, pb, pc, pd, aheight, bheight, cheight, dheight)
6389 REAL adx, bdx, cdx, ady, bdy, cdy, adheight, bdheight, cdheight;
6390 REAL bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady;
6392 REAL permanent, errbound;
6396 adx = pa[0] - pd[0];
6397 bdx = pb[0] - pd[0];
6398 cdx = pc[0] - pd[0];
6399 ady = pa[1] - pd[1];
6400 bdy = pb[1] - pd[1];
6401 cdy = pc[1] - pd[1];
6402 adheight = aheight - dheight;
6403 bdheight = bheight - dheight;
6404 cdheight = cheight - dheight;
6415 det = adheight * (bdxcdy - cdxbdy)
6416 + bdheight * (cdxady - adxcdy)
6417 + cdheight * (adxbdy - bdxady);
6427 if ((det > errbound) || (-det > errbound)) {
6431 return orient3dadapt(pa, pb, pc, pd, aheight, bheight, cheight, dheight,
6453#ifdef ANSI_DECLARATORS
6467 if (
b->weighted == 0) {
6469 }
else if (
b->weighted == 1) {
6471 pa[0] * pa[0] + pa[1] * pa[1] - pa[2],
6472 pb[0] * pb[0] + pb[1] * pb[1] - pb[2],
6473 pc[0] * pc[0] + pc[1] * pc[1] - pc[2],
6474 pd[0] * pd[0] + pd[1] * pd[1] - pd[2]);
6476 return orient3d(
m,
b, pa, pb, pc, pd, pa[2], pb[2], pc[2], pd[2]);
6494#ifdef ANSI_DECLARATORS
6513 REAL xdo, ydo, xao, yao;
6514 REAL dodist, aodist, dadist;
6516 REAL dx, dy, dxoff, dyoff;
6518 m->circumcentercount++;
6521 xdo = tdest[0] - torg[0];
6522 ydo = tdest[1] - torg[1];
6523 xao = tapex[0] - torg[0];
6524 yao = tapex[1] - torg[1];
6525 dodist = xdo * xdo + ydo * ydo;
6526 aodist = xao * xao + yao * yao;
6527 dadist = (tdest[0] - tapex[0]) * (tdest[0] - tapex[0]) +
6528 (tdest[1] - tapex[1]) * (tdest[1] - tapex[1]);
6530 denominator = 0.5 / (xdo * yao - xao * ydo);
6537 m->counterclockcount--;
6539 dx = (yao * dodist - ydo * aodist) * denominator;
6540 dy = (xdo * aodist - xao * dodist) * denominator;
6547 if ((dodist < aodist) && (dodist < dadist)) {
6548 if (offcenter && (
b->offconstant > 0.0)) {
6550 dxoff = 0.5 * xdo -
b->offconstant * ydo;
6551 dyoff = 0.5 * ydo +
b->offconstant * xdo;
6554 if (dxoff * dxoff + dyoff * dyoff < dx * dx + dy * dy) {
6559 }
else if (aodist < dadist) {
6560 if (offcenter && (
b->offconstant > 0.0)) {
6561 dxoff = 0.5 * xao +
b->offconstant * yao;
6562 dyoff = 0.5 * yao -
b->offconstant * xao;
6565 if (dxoff * dxoff + dyoff * dyoff < dx * dx + dy * dy) {
6571 if (offcenter && (
b->offconstant > 0.0)) {
6572 dxoff = 0.5 * (tapex[0] - tdest[0]) -
6573 b->offconstant * (tapex[1] - tdest[1]);
6574 dyoff = 0.5 * (tapex[1] - tdest[1]) +
6575 b->offconstant * (tapex[0] - tdest[0]);
6578 if (dxoff * dxoff + dyoff * dyoff <
6579 (dx - xdo) * (dx - xdo) + (dy - ydo) * (dy - ydo)) {
6586 circumcenter[0] = torg[0] + dx;
6587 circumcenter[1] = torg[1] + dy;
6594 *xi = (yao * dx - xao * dy) * (2.0 * denominator);
6595 *eta = (xdo * dy - ydo * dx) * (2.0 * denominator);
6608#ifdef ANSI_DECLARATORS
6628 m->checksegments = 0;
6629 m->checkquality = 0;
6630 m->incirclecount =
m->counterclockcount =
m->orient3dcount = 0;
6631 m->hyperbolacount =
m->circletopcount =
m->circumcentercount = 0;
6647#ifdef ANSI_DECLARATORS
6651unsigned int choices;
6671#ifdef ANSI_DECLARATORS
6680 struct otri triangleloop;
6681 struct otri oppotri, oppooppotri;
6682 vertex triorg, tridest, triapex;
6683 vertex oppoorg, oppodest;
6689 saveexact =
b->noexact;
6692 printf(
" Checking consistency of mesh...\n");
6698 while (triangleloop.tri != (
triangle *) NULL) {
6700 for (triangleloop.orient = 0; triangleloop.orient < 3;
6701 triangleloop.orient++) {
6702 org(triangleloop, triorg);
6703 dest(triangleloop, tridest);
6704 if (triangleloop.orient == 0) {
6706 apex(triangleloop, triapex);
6708 printf(
" !! !! Inverted ");
6714 sym(triangleloop, oppotri);
6715 if (oppotri.tri !=
m->dummytri) {
6717 sym(oppotri, oppooppotri);
6718 if ((triangleloop.tri != oppooppotri.tri)
6719 || (triangleloop.orient != oppooppotri.orient)) {
6720 printf(
" !! !! Asymmetric triangle-triangle bond:\n");
6721 if (triangleloop.tri == oppooppotri.tri) {
6722 printf(
" (Right triangle, wrong orientation)\n");
6726 printf(
" Second (nonreciprocating) ");
6732 org(oppotri, oppoorg);
6733 dest(oppotri, oppodest);
6734 if ((triorg != oppodest) || (tridest != oppoorg)) {
6735 printf(
" !! !! Mismatched edge coordinates between two triangles:\n"
6737 printf(
" First mismatched ");
6739 printf(
" Second mismatched ");
6749 printf(
" In my studied opinion, the mesh appears to be consistent.\n");
6751 }
else if (horrors == 1) {
6752 printf(
" !! !! !! !! Precisely one festering wound discovered.\n");
6754 printf(
" !! !! !! !! %d abominations witnessed.\n", horrors);
6757 b->noexact = saveexact;
6770#ifdef ANSI_DECLARATORS
6773void checkdelaunay(
m,
b)
6779 struct otri triangleloop;
6780 struct otri oppotri;
6781 struct osub opposubseg;
6782 vertex triorg, tridest, triapex;
6784 int shouldbedelaunay;
6791 saveexact =
b->noexact;
6794 printf(
" Checking Delaunay property of mesh...\n");
6800 while (triangleloop.tri != (
triangle *) NULL) {
6802 for (triangleloop.orient = 0; triangleloop.orient < 3;
6803 triangleloop.orient++) {
6804 org(triangleloop, triorg);
6805 dest(triangleloop, tridest);
6806 apex(triangleloop, triapex);
6807 sym(triangleloop, oppotri);
6808 apex(oppotri, oppoapex);
6812 shouldbedelaunay = (oppotri.tri !=
m->dummytri) &&
6813 !
deadtri(oppotri.tri) && (triangleloop.tri < oppotri.tri) &&
6814 (triorg !=
m->infvertex1) && (triorg !=
m->infvertex2) &&
6815 (triorg !=
m->infvertex3) &&
6816 (tridest !=
m->infvertex1) && (tridest !=
m->infvertex2) &&
6817 (tridest !=
m->infvertex3) &&
6818 (triapex !=
m->infvertex1) && (triapex !=
m->infvertex2) &&
6819 (triapex !=
m->infvertex3) &&
6820 (oppoapex !=
m->infvertex1) && (oppoapex !=
m->infvertex2) &&
6821 (oppoapex !=
m->infvertex3);
6822 if (
m->checksegments && shouldbedelaunay) {
6825 tspivot(triangleloop, opposubseg);
6826 if (opposubseg.ss !=
m->dummysub){
6827 shouldbedelaunay = 0;
6830 if (shouldbedelaunay) {
6831 if (
nonregular(
m,
b, triorg, tridest, triapex, oppoapex) > 0.0) {
6833 printf(
" !! !! Non-Delaunay pair of triangles:\n");
6834 printf(
" First non-Delaunay ");
6836 printf(
" Second non-Delaunay ");
6838 printf(
" !! !! Non-regular pair of triangles:\n");
6839 printf(
" First non-regular ");
6841 printf(
" Second non-regular ");
6853 " By virtue of my perceptive intelligence, I declare the mesh Delaunay.\n");
6855 }
else if (horrors == 1) {
6857 " !! !! !! !! Precisely one terrifying transgression identified.\n");
6859 printf(
" !! !! !! !! %d obscenities viewed with horror.\n", horrors);
6862 b->noexact = saveexact;
6880#ifdef ANSI_DECLARATORS
6884void enqueuebadtriang(
m,
b, badtri)
6892 int exponent, expincrement;
6897 if (
b->verbose > 2) {
6898 printf(
" Queueing bad triangle:\n");
6899 printf(
" (%.12g, %.12g) (%.12g, %.12g) (%.12g, %.12g)\n",
6907 if (badtri->
key >= 1.0) {
6923 while (
length * multiplier * multiplier > 1.0) {
6925 multiplier *= multiplier;
6928 exponent += expincrement;
6937 queuenumber = 2047 - exponent;
6939 queuenumber = 2048 + exponent;
6943 if (
m->queuefront[queuenumber] == (
struct badtriang *) NULL) {
6946 if (queuenumber >
m->firstnonemptyq) {
6948 m->nextnonemptyq[queuenumber] =
m->firstnonemptyq;
6949 m->firstnonemptyq = queuenumber;
6953 i = queuenumber + 1;
6954 while (
m->queuefront[i] == (
struct badtriang *) NULL) {
6958 m->nextnonemptyq[queuenumber] =
m->nextnonemptyq[i];
6959 m->nextnonemptyq[i] = queuenumber;
6962 m->queuefront[queuenumber] = badtri;
6965 m->queuetail[queuenumber]->nexttriang = badtri;
6968 m->queuetail[queuenumber] = badtri;
6986#ifdef ANSI_DECLARATORS
6990void enqueuebadtri(
m,
b, enqtri, minedge, enqapex, enqorg, enqdest)
7006 newbad->
key = minedge;
7010 enqueuebadtriang(
m,
b, newbad);
7023#ifdef ANSI_DECLARATORS
7034 if (
m->firstnonemptyq < 0) {
7038 result =
m->queuefront[
m->firstnonemptyq];
7040 m->queuefront[
m->firstnonemptyq] =
result->nexttriang;
7043 if (
result ==
m->queuetail[
m->firstnonemptyq]) {
7044 m->firstnonemptyq =
m->nextnonemptyq[
m->firstnonemptyq];
7076#ifdef ANSI_DECLARATORS
7078 struct osub *testsubseg)
7080int checkseg4encroach(
m,
b, testsubseg)
7083struct osub *testsubseg;
7087 struct otri neighbortri;
7088 struct osub testsym;
7093 vertex eorg, edest, eapex;
7099 sorg(*testsubseg, eorg);
7100 sdest(*testsubseg, edest);
7102 stpivot(*testsubseg, neighbortri);
7104 if (neighbortri.tri !=
m->dummytri) {
7107 apex(neighbortri, eapex);
7113 dotproduct = (eorg[0] - eapex[0]) * (edest[0] - eapex[0]) +
7114 (eorg[1] - eapex[1]) * (edest[1] - eapex[1]);
7115 if (dotproduct < 0.0) {
7116 if (
b->conformdel ||
7117 (dotproduct * dotproduct >=
7118 (2.0 *
b->goodangle - 1.0) * (2.0 *
b->goodangle - 1.0) *
7119 ((eorg[0] - eapex[0]) * (eorg[0] - eapex[0]) +
7120 (eorg[1] - eapex[1]) * (eorg[1] - eapex[1])) *
7121 ((edest[0] - eapex[0]) * (edest[0] - eapex[0]) +
7122 (edest[1] - eapex[1]) * (edest[1] - eapex[1])))) {
7128 ssym(*testsubseg, testsym);
7129 stpivot(testsym, neighbortri);
7131 if (neighbortri.tri !=
m->dummytri) {
7134 apex(neighbortri, eapex);
7137 dotproduct = (eorg[0] - eapex[0]) * (edest[0] - eapex[0]) +
7138 (eorg[1] - eapex[1]) * (edest[1] - eapex[1]);
7139 if (dotproduct < 0.0) {
7140 if (
b->conformdel ||
7141 (dotproduct * dotproduct >=
7142 (2.0 *
b->goodangle - 1.0) * (2.0 *
b->goodangle - 1.0) *
7143 ((eorg[0] - eapex[0]) * (eorg[0] - eapex[0]) +
7144 (eorg[1] - eapex[1]) * (eorg[1] - eapex[1])) *
7145 ((edest[0] - eapex[0]) * (edest[0] - eapex[0]) +
7146 (edest[1] - eapex[1]) * (edest[1] - eapex[1])))) {
7152 if (encroached && (!
b->nobisect || ((
b->nobisect == 1) && (sides == 2)))) {
7153 if (
b->verbose > 2) {
7155 " Queueing encroached subsegment (%.12g, %.12g) (%.12g, %.12g).\n",
7156 eorg[0], eorg[1], edest[0], edest[1]);
7161 if (encroached == 1) {
7189#ifdef ANSI_DECLARATORS
7192void testtriangle(
m,
b, testtri)
7195struct otri *testtri;
7199 struct otri tri1, tri2;
7200 struct osub testsub;
7201 vertex torg, tdest, tapex;
7203 vertex org1, dest1, org2, dest2;
7205 REAL dxod, dyod, dxda, dyda, dxao, dyao;
7206 REAL dxod2, dyod2, dxda2, dyda2, dxao2, dyao2;
7207 REAL apexlen, orglen, destlen, minedge;
7214 org(*testtri, torg);
7215 dest(*testtri, tdest);
7216 apex(*testtri, tapex);
7217 dxod = torg[0] - tdest[0];
7218 dyod = torg[1] - tdest[1];
7219 dxda = tdest[0] - tapex[0];
7220 dyda = tdest[1] - tapex[1];
7221 dxao = tapex[0] - torg[0];
7222 dyao = tapex[1] - torg[1];
7223 dxod2 = dxod * dxod;
7224 dyod2 = dyod * dyod;
7225 dxda2 = dxda * dxda;
7226 dyda2 = dyda * dyda;
7227 dxao2 = dxao * dxao;
7228 dyao2 = dyao * dyao;
7230 apexlen = dxod2 + dyod2;
7231 orglen = dxda2 + dyda2;
7232 destlen = dxao2 + dyao2;
7234 if ((apexlen < orglen) && (apexlen < destlen)) {
7238 angle = dxda * dxao + dyda * dyao;
7243 }
else if (orglen < destlen) {
7247 angle = dxod * dxao + dyod * dyao;
7251 lnext(*testtri, tri1);
7256 angle = dxod * dxda + dyod * dyda;
7260 lprev(*testtri, tri1);
7263 if (
b->vararea ||
b->fixedarea ||
b->usertest) {
7265 area = 0.5 * (dxod * dyda - dyod * dxda);
7266 if (
b->fixedarea && (area >
b->maxarea)) {
7268 enqueuebadtri(
m,
b, testtri, minedge, tapex, torg, tdest);
7273 if ((
b->vararea) && (area >
areabound(*testtri)) &&
7276 enqueuebadtri(
m,
b, testtri, minedge, tapex, torg, tdest);
7283 enqueuebadtri(
m,
b, testtri, minedge, tapex, torg, tdest);
7290 if (
angle >
b->goodangle) {
7305 if (testsub.ss ==
m->dummysub) {
7311 }
while (testsub.ss ==
m->dummysub);
7319 }
while (testsub.ss ==
m->dummysub);
7324 joinvertex = (
vertex) NULL;
7325 if ((dest1[0] == org2[0]) && (dest1[1] == org2[1])) {
7327 }
else if ((org1[0] == dest2[0]) && (org1[1] == dest2[1])) {
7330 if (joinvertex != (
vertex) NULL) {
7333 dist1 = ((base1[0] - joinvertex[0]) * (base1[0] - joinvertex[0]) +
7334 (base1[1] - joinvertex[1]) * (base1[1] - joinvertex[1]));
7335 dist2 = ((base2[0] - joinvertex[0]) * (base2[0] - joinvertex[0]) +
7336 (base2[1] - joinvertex[1]) * (base2[1] - joinvertex[1]));
7338 if ((dist1 < 1.001 * dist2) && (dist1 > 0.999 * dist2)) {
7347 enqueuebadtri(
m,
b, testtri, minedge, tapex, torg, tdest);
7375#ifdef ANSI_DECLARATORS
7384 struct otri triangleloop;
7388 printf(
" Constructing mapping from vertices to triangles.\n");
7396 org(triangleloop, triorg);
7470#ifdef ANSI_DECLARATORS
7473 int stopatsubsegment)
7479struct otri *searchtri;
7480int stopatsubsegment;
7484 struct otri backtracktri;
7485 struct osub checkedge;
7486 vertex forg, fdest, fapex;
7487 REAL orgorient, destorient;
7492 if (
b->verbose > 2) {
7493 printf(
" Searching for point (%.12g, %.12g).\n",
7494 searchpoint[0], searchpoint[1]);
7497 org(*searchtri, forg);
7498 dest(*searchtri, fdest);
7499 apex(*searchtri, fapex);
7501 if (
b->verbose > 2) {
7502 printf(
" At (%.12g, %.12g) (%.12g, %.12g) (%.12g, %.12g)\n",
7503 forg[0], forg[1], fdest[0], fdest[1], fapex[0], fapex[1]);
7506 if ((fapex[0] == searchpoint[0]) && (fapex[1] == searchpoint[1])) {
7516 if (destorient > 0.0) {
7517 if (orgorient > 0.0) {
7523 moveleft = (fapex[0] - searchpoint[0]) * (fdest[0] - forg[0]) +
7524 (fapex[1] - searchpoint[1]) * (fdest[1] - forg[1]) > 0.0;
7529 if (orgorient > 0.0) {
7534 if (destorient == 0.0) {
7538 if (orgorient == 0.0) {
7550 lprev(*searchtri, backtracktri);
7553 lnext(*searchtri, backtracktri);
7556 sym(backtracktri, *searchtri);
7558 if (
m->checksegments && stopatsubsegment) {
7560 tspivot(backtracktri, checkedge);
7561 if (checkedge.
ss !=
m->dummysub) {
7563 otricopy(backtracktri, *searchtri);
7568 if (searchtri->
tri ==
m->dummytri) {
7570 otricopy(backtracktri, *searchtri);
7574 apex(*searchtri, fapex);
7614#ifdef ANSI_DECLARATORS
7622struct otri *searchtri;
7628 struct otri sampletri;
7631 REAL searchdist, dist;
7633 long samplesperblock, totalsamplesleft, samplesleft;
7634 long population, totalpopulation;
7637 if (
b->verbose > 2) {
7638 printf(
" Randomly sampling for a triangle near point (%.12g, %.12g).\n",
7639 searchpoint[0], searchpoint[1]);
7643 org(*searchtri, torg);
7644 searchdist = (searchpoint[0] - torg[0]) * (searchpoint[0] - torg[0]) +
7645 (searchpoint[1] - torg[1]) * (searchpoint[1] - torg[1]);
7646 if (
b->verbose > 2) {
7647 printf(
" Boundary triangle has origin (%.12g, %.12g).\n",
7653 if (
m->recenttri.tri != (
triangle *) NULL) {
7655 org(
m->recenttri, torg);
7656 if ((torg[0] == searchpoint[0]) && (torg[1] == searchpoint[1])) {
7660 dist = (searchpoint[0] - torg[0]) * (searchpoint[0] - torg[0]) +
7661 (searchpoint[1] - torg[1]) * (searchpoint[1] - torg[1]);
7662 if (dist < searchdist) {
7665 if (
b->verbose > 2) {
7666 printf(
" Choosing recent triangle with origin (%.12g, %.12g).\n",
7678 m->triangles.items) {
7686 samplesperblock = (
m->samples *
TRIPERBLOCK - 1) /
m->triangles.maxitems + 1;
7689 samplesleft = (
m->samples *
m->triangles.itemsfirstblock - 1) /
7690 m->triangles.maxitems + 1;
7691 totalsamplesleft =
m->samples;
7692 population =
m->triangles.itemsfirstblock;
7693 totalpopulation =
m->triangles.maxitems;
7694 sampleblock =
m->triangles.firstblock;
7696 while (totalsamplesleft > 0) {
7698 if (population > totalpopulation) {
7699 population = totalpopulation;
7702 alignptr = (uintptr_t) (sampleblock + 1);
7703 firsttri = (
char *) (alignptr +
7704 (uintptr_t)
m->triangles.alignbytes -
7706 (uintptr_t)
m->triangles.alignbytes));
7712 m->triangles.itembytes));
7714 org(sampletri, torg);
7715 dist = (searchpoint[0] - torg[0]) * (searchpoint[0] - torg[0]) +
7716 (searchpoint[1] - torg[1]) * (searchpoint[1] - torg[1]);
7717 if (dist < searchdist) {
7720 if (
b->verbose > 2) {
7721 printf(
" Choosing triangle with origin (%.12g, %.12g).\n",
7729 }
while ((samplesleft > 0) && (totalsamplesleft > 0));
7731 if (totalsamplesleft > 0) {
7732 sampleblock = (
VOID **) *sampleblock;
7733 samplesleft = samplesperblock;
7734 totalpopulation -= population;
7740 org(*searchtri, torg);
7741 dest(*searchtri, tdest);
7743 if ((torg[0] == searchpoint[0]) && (torg[1] == searchpoint[1])) {
7746 if ((tdest[0] == searchpoint[0]) && (tdest[1] == searchpoint[1])) {
7756 }
else if (ahead == 0.0) {
7758 if (((torg[0] < searchpoint[0]) == (searchpoint[0] < tdest[0])) &&
7759 ((torg[1] < searchpoint[1]) == (searchpoint[1] < tdest[1]))) {
7785#ifdef ANSI_DECLARATORS
7797 struct otri oppotri;
7798 struct osub newsubseg;
7804 dest(*tri, tridest);
7814 if (newsubseg.
ss ==
m->dummysub) {
7828 tsbond(oppotri, newsubseg);
7829 setmark(newsubseg, subsegmark);
7830 if (
b->verbose > 2) {
7831 printf(
" Inserting new ");
7835 if (
mark(newsubseg) == 0) {
7836 setmark(newsubseg, subsegmark);
7889#ifdef ANSI_DECLARATORS
7895struct otri *flipedge;
7899 struct otri botleft, botright;
7900 struct otri topleft, topright;
7902 struct otri botlcasing, botrcasing;
7903 struct otri toplcasing, toprcasing;
7904 struct osub botlsubseg, botrsubseg;
7905 struct osub toplsubseg, toprsubseg;
7906 vertex leftvertex, rightvertex, botvertex;
7912 org(*flipedge, rightvertex);
7913 dest(*flipedge, leftvertex);
7914 apex(*flipedge, botvertex);
7915 sym(*flipedge, top);
7917 if (top.
tri ==
m->dummytri) {
7918 printf(
"Internal error in flip(): Attempt to flip on boundary.\n");
7922 if (
m->checksegments) {
7923 tspivot(*flipedge, toplsubseg);
7924 if (toplsubseg.
ss !=
m->dummysub) {
7925 printf(
"Internal error in flip(): Attempt to flip a segment.\n");
7931 apex(top, farvertex);
7934 lprev(top, topleft);
7935 sym(topleft, toplcasing);
7936 lnext(top, topright);
7937 sym(topright, toprcasing);
7938 lnext(*flipedge, botleft);
7939 sym(botleft, botlcasing);
7940 lprev(*flipedge, botright);
7941 sym(botright, botrcasing);
7943 bond(topleft, botlcasing);
7944 bond(botleft, botrcasing);
7945 bond(botright, toprcasing);
7946 bond(topright, toplcasing);
7948 if (
m->checksegments) {
7952 tspivot(botright, botrsubseg);
7953 tspivot(topright, toprsubseg);
7954 if (toplsubseg.
ss ==
m->dummysub) {
7957 tsbond(topright, toplsubseg);
7959 if (botlsubseg.
ss ==
m->dummysub) {
7962 tsbond(topleft, botlsubseg);
7964 if (botrsubseg.
ss ==
m->dummysub) {
7967 tsbond(botleft, botrsubseg);
7969 if (toprsubseg.
ss ==
m->dummysub) {
7972 tsbond(botright, toprsubseg);
7977 setorg(*flipedge, farvertex);
7978 setdest(*flipedge, botvertex);
7979 setapex(*flipedge, rightvertex);
7983 if (
b->verbose > 2) {
7984 printf(
" Edge flip results in left ");
7986 printf(
" and right ");
8024#ifdef ANSI_DECLARATORS
8030struct otri *flipedge;
8034 struct otri botleft, botright;
8035 struct otri topleft, topright;
8037 struct otri botlcasing, botrcasing;
8038 struct otri toplcasing, toprcasing;
8039 struct osub botlsubseg, botrsubseg;
8040 struct osub toplsubseg, toprsubseg;
8041 vertex leftvertex, rightvertex, botvertex;
8047 org(*flipedge, rightvertex);
8048 dest(*flipedge, leftvertex);
8049 apex(*flipedge, botvertex);
8050 sym(*flipedge, top);
8052 if (top.
tri ==
m->dummytri) {
8053 printf(
"Internal error in unflip(): Attempt to flip on boundary.\n");
8057 if (
m->checksegments) {
8058 tspivot(*flipedge, toplsubseg);
8059 if (toplsubseg.
ss !=
m->dummysub) {
8060 printf(
"Internal error in unflip(): Attempt to flip a subsegment.\n");
8066 apex(top, farvertex);
8069 lprev(top, topleft);
8070 sym(topleft, toplcasing);
8071 lnext(top, topright);
8072 sym(topright, toprcasing);
8073 lnext(*flipedge, botleft);
8074 sym(botleft, botlcasing);
8075 lprev(*flipedge, botright);
8076 sym(botright, botrcasing);
8078 bond(topleft, toprcasing);
8079 bond(botleft, toplcasing);
8080 bond(botright, botlcasing);
8081 bond(topright, botrcasing);
8083 if (
m->checksegments) {
8087 tspivot(botright, botrsubseg);
8088 tspivot(topright, toprsubseg);
8089 if (toplsubseg.
ss ==
m->dummysub) {
8092 tsbond(botleft, toplsubseg);
8094 if (botlsubseg.
ss ==
m->dummysub) {
8097 tsbond(botright, botlsubseg);
8099 if (botrsubseg.
ss ==
m->dummysub) {
8102 tsbond(topright, botrsubseg);
8104 if (toprsubseg.
ss ==
m->dummysub) {
8107 tsbond(topleft, toprsubseg);
8112 setorg(*flipedge, botvertex);
8113 setdest(*flipedge, farvertex);
8114 setapex(*flipedge, leftvertex);
8118 if (
b->verbose > 2) {
8119 printf(
" Edge unflip results in left ");
8121 printf(
" and right ");
8173#ifdef ANSI_DECLARATORS
8176 struct osub *splitseg,
8177 int segmentflaws,
int triflaws )
8180 segmentflaws, triflaws)
8184struct otri *searchtri;
8185struct osub *splitseg;
8194 struct otri botleft, botright;
8195 struct otri topleft, topright;
8196 struct otri newbotleft, newbotright;
8197 struct otri newtopright;
8198 struct otri botlcasing, botrcasing;
8199 struct otri toplcasing, toprcasing;
8200 struct otri testtri;
8201 struct osub botlsubseg, botrsubseg;
8202 struct osub toplsubseg, toprsubseg;
8203 struct osub brokensubseg;
8204 struct osub checksubseg;
8205 struct osub rightsubseg;
8206 struct osub newsubseg;
8210 vertex leftvertex, rightvertex, botvertex, topvertex, farvertex;
8211 vertex segmentorg, segmentdest;
8226 if (
b->verbose > 1) {
8227 printf(
" Inserting (%.12g, %.12g).\n", newvertex[0], newvertex[1]);
8230 if (splitseg == (
struct osub *) NULL) {
8233 if (searchtri->
tri ==
m->dummytri) {
8235 horiz.
tri =
m->dummytri;
8239 intersect =
locate(
m,
b, newvertex, &horiz);
8261 if (
m->checksegments && (splitseg == (
struct osub *) NULL)) {
8264 if (brokensubseg.
ss !=
m->dummysub) {
8267 enq =
b->nobisect != 2;
8268 if (enq && (
b->nobisect == 1)) {
8271 sym(horiz, testtri);
8272 enq = testtri.
tri !=
m->dummytri;
8280 if (
b->verbose > 2) {
8282 " Queueing encroached subsegment (%.12g, %.12g) (%.12g, %.12g).\n",
8298 lprev(horiz, botright);
8299 sym(botright, botrcasing);
8300 sym(horiz, topright);
8302 mirrorflag = topright.
tri !=
m->dummytri;
8305 sym(topright, toprcasing);
8314 org(horiz, rightvertex);
8315 dest(horiz, leftvertex);
8316 apex(horiz, botvertex);
8317 setorg(newbotright, botvertex);
8318 setdest(newbotright, rightvertex);
8319 setapex(newbotright, newvertex);
8320 setorg(horiz, newvertex);
8321 for (i = 0; i <
m->eextras; i++) {
8330 dest(topright, topvertex);
8331 setorg(newtopright, rightvertex);
8332 setdest(newtopright, topvertex);
8333 setapex(newtopright, newvertex);
8334 setorg(topright, newvertex);
8335 for (i = 0; i <
m->eextras; i++) {
8347 if (
m->checksegments) {
8348 tspivot(botright, botrsubseg);
8349 if (botrsubseg.
ss !=
m->dummysub) {
8351 tsbond(newbotright, botrsubseg);
8354 tspivot(topright, toprsubseg);
8355 if (toprsubseg.
ss !=
m->dummysub) {
8357 tsbond(newtopright, toprsubseg);
8363 bond(newbotright, botrcasing);
8365 bond(newbotright, botright);
8368 bond(newtopright, toprcasing);
8370 bond(newtopright, topright);
8372 bond(newtopright, newbotright);
8375 if (splitseg != (
struct osub *) NULL) {
8378 segorg(*splitseg, segmentorg);
8379 segdest(*splitseg, segmentdest);
8381 spivot(*splitseg, rightsubseg);
8383 tspivot(newbotright, newsubseg);
8386 sbond(*splitseg, newsubseg);
8388 sbond(newsubseg, rightsubseg);
8397 if (
m->checkquality) {
8400 m->lastflip->flippedtri =
encode(horiz);
8406 printf(
"Internal error in insertvertex():\n");
8408 " Clockwise triangle prior to edge vertex insertion (bottom).\n");
8412 printf(
"Internal error in insertvertex():\n");
8413 printf(
" Clockwise triangle prior to edge vertex insertion (top).\n");
8416 printf(
"Internal error in insertvertex():\n");
8418 " Clockwise triangle after edge vertex insertion (top right).\n");
8421 printf(
"Internal error in insertvertex():\n");
8423 " Clockwise triangle after edge vertex insertion (top left).\n");
8427 printf(
"Internal error in insertvertex():\n");
8429 " Clockwise triangle after edge vertex insertion (bottom left).\n");
8432 printf(
"Internal error in insertvertex():\n");
8434 " Clockwise triangle after edge vertex insertion (bottom right).\n");
8437 if (
b->verbose > 2) {
8438 printf(
" Updating bottom left ");
8441 printf(
" Updating top left ");
8443 printf(
" Creating top right ");
8446 printf(
" Creating bottom right ");
8455 lnext(horiz, botleft);
8456 lprev(horiz, botright);
8457 sym(botleft, botlcasing);
8458 sym(botright, botrcasing);
8463 org(horiz, rightvertex);
8464 dest(horiz, leftvertex);
8465 apex(horiz, botvertex);
8466 setorg(newbotleft, leftvertex);
8467 setdest(newbotleft, botvertex);
8468 setapex(newbotleft, newvertex);
8469 setorg(newbotright, botvertex);
8470 setdest(newbotright, rightvertex);
8471 setapex(newbotright, newvertex);
8473 for (i = 0; i <
m->eextras; i++) {
8488 if (
m->checksegments) {
8490 if (botlsubseg.
ss !=
m->dummysub) {
8492 tsbond(newbotleft, botlsubseg);
8494 tspivot(botright, botrsubseg);
8495 if (botrsubseg.
ss !=
m->dummysub) {
8497 tsbond(newbotright, botrsubseg);
8502 bond(newbotleft, botlcasing);
8503 bond(newbotright, botrcasing);
8506 bond(newbotleft, newbotright);
8508 bond(botleft, newbotleft);
8510 bond(botright, newbotright);
8512 if (
m->checkquality) {
8515 m->lastflip->flippedtri =
encode(horiz);
8521 printf(
"Internal error in insertvertex():\n");
8522 printf(
" Clockwise triangle prior to vertex insertion.\n");
8525 printf(
"Internal error in insertvertex():\n");
8526 printf(
" Clockwise triangle after vertex insertion (top).\n");
8529 printf(
"Internal error in insertvertex():\n");
8530 printf(
" Clockwise triangle after vertex insertion (left).\n");
8533 printf(
"Internal error in insertvertex():\n");
8534 printf(
" Clockwise triangle after vertex insertion (right).\n");
8537 if (
b->verbose > 2) {
8538 printf(
" Updating top ");
8540 printf(
" Creating left ");
8542 printf(
" Creating right ");
8555 rightvertex =
first;
8556 dest(horiz, leftvertex);
8562 if (
m->checksegments) {
8565 if (checksubseg.
ss !=
m->dummysub) {
8571 if (checkseg4encroach(
m,
b, &checksubseg)) {
8582 if (top.
tri ==
m->dummytri) {
8587 apex(top, farvertex);
8593 if ((leftvertex ==
m->infvertex1) || (leftvertex ==
m->infvertex2) ||
8594 (leftvertex ==
m->infvertex3)) {
8601 }
else if ((rightvertex ==
m->infvertex1) ||
8602 (rightvertex ==
m->infvertex2) ||
8603 (rightvertex ==
m->infvertex3)) {
8610 }
else if ((farvertex ==
m->infvertex1) ||
8611 (farvertex ==
m->infvertex2) ||
8612 (farvertex ==
m->infvertex3)) {
8618 doflip =
incircle(
m,
b, leftvertex, newvertex, rightvertex,
8625 lprev(top, topleft);
8626 sym(topleft, toplcasing);
8627 lnext(top, topright);
8628 sym(topright, toprcasing);
8629 lnext(horiz, botleft);
8630 sym(botleft, botlcasing);
8631 lprev(horiz, botright);
8632 sym(botright, botrcasing);
8634 bond(topleft, botlcasing);
8635 bond(botleft, botrcasing);
8636 bond(botright, toprcasing);
8637 bond(topright, toplcasing);
8638 if (
m->checksegments) {
8642 tspivot(botright, botrsubseg);
8643 tspivot(topright, toprsubseg);
8644 if (toplsubseg.
ss ==
m->dummysub) {
8647 tsbond(topright, toplsubseg);
8649 if (botlsubseg.
ss ==
m->dummysub) {
8652 tsbond(topleft, botlsubseg);
8654 if (botrsubseg.
ss ==
m->dummysub) {
8657 tsbond(botleft, botrsubseg);
8659 if (toprsubseg.
ss ==
m->dummysub) {
8662 tsbond(botright, toprsubseg);
8666 setorg(horiz, farvertex);
8672 for (i = 0; i <
m->eextras; i++) {
8691 if (
m->checkquality) {
8695 m->lastflip = newflip;
8699 if (newvertex != (
vertex) NULL) {
8702 printf(
"Internal error in insertvertex():\n");
8703 printf(
" Clockwise triangle prior to edge flip (bottom).\n");
8717 printf(
"Internal error in insertvertex():\n");
8718 printf(
" Clockwise triangle after edge flip (left).\n");
8722 printf(
"Internal error in insertvertex():\n");
8723 printf(
" Clockwise triangle after edge flip (right).\n");
8727 if (
b->verbose > 2) {
8728 printf(
" Edge flip results in left ");
8731 printf(
" and right ");
8738 leftvertex = farvertex;
8747 testtriangle(
m,
b, &horiz);
8752 sym(horiz, testtri);
8756 if ((leftvertex ==
first) || (testtri.
tri ==
m->dummytri)) {
8758 lnext(horiz, *searchtri);
8759 lnext(horiz,
m->recenttri);
8763 lnext(testtri, horiz);
8764 rightvertex = leftvertex;
8765 dest(horiz, leftvertex);
8834#ifdef ANSI_DECLARATORS
8836 struct otri *firstedge,
struct otri *lastedge,
8837 int edgecount,
int doflip,
int triflaws)
8842struct otri *firstedge;
8843struct otri *lastedge;
8850 struct otri testtri;
8851 struct otri besttri;
8852 struct otri tempedge;
8853 vertex leftbasevertex, rightbasevertex;
8861 apex(*lastedge, leftbasevertex);
8862 dest(*firstedge, rightbasevertex);
8863 if (
b->verbose > 2) {
8864 printf(
" Triangulating interior polygon at edge\n");
8865 printf(
" (%.12g, %.12g) (%.12g, %.12g)\n", leftbasevertex[0],
8866 leftbasevertex[1], rightbasevertex[0], rightbasevertex[1]);
8869 onext(*firstedge, besttri);
8870 dest(besttri, bestvertex);
8873 for (i = 2; i <= edgecount - 2; i++) {
8875 dest(testtri, testvertex);
8877 if (
incircle(
m,
b, leftbasevertex, rightbasevertex, bestvertex,
8878 testvertex) > 0.0) {
8880 bestvertex = testvertex;
8884 if (
b->verbose > 2) {
8885 printf(
" Connecting edge to (%.12g, %.12g)\n", bestvertex[0],
8888 if (bestnumber > 1) {
8890 oprev(besttri, tempedge);
8894 if (bestnumber < edgecount - 2) {
8896 sym(besttri, tempedge);
8900 sym(tempedge, besttri);
8908 sym(besttri, testtri);
8909 testtriangle(
m,
b, &testtri);
8933#ifdef ANSI_DECLARATORS
8936void deletevertex(
m,
b, deltri)
8943 struct otri countingtri;
8944 struct otri firstedge, lastedge;
8945 struct otri deltriright;
8946 struct otri lefttri, righttri;
8947 struct otri leftcasing, rightcasing;
8948 struct osub leftsubseg, rightsubseg;
8955 org(*deltri, delvertex);
8956 if (
b->verbose > 1) {
8957 printf(
" Deleting (%.12g, %.12g).\n", delvertex[0], delvertex[1]);
8962 onext(*deltri, countingtri);
8964 while (!
otriequal(*deltri, countingtri)) {
8966 if (countingtri.tri ==
m->dummytri) {
8967 printf(
"Internal error in deletevertex():\n");
8968 printf(
" Attempt to delete boundary vertex.\n");
8977 if (edgecount < 3) {
8978 printf(
"Internal error in deletevertex():\n Vertex has degree %d.\n",
8983 if (edgecount > 3) {
8987 onext(*deltri, firstedge);
8988 oprev(*deltri, lastedge);
8993 lprev(*deltri, deltriright);
8994 dnext(*deltri, lefttri);
8995 sym(lefttri, leftcasing);
8996 oprev(deltriright, righttri);
8997 sym(righttri, rightcasing);
8998 bond(*deltri, leftcasing);
8999 bond(deltriright, rightcasing);
9001 if (leftsubseg.ss !=
m->dummysub) {
9002 tsbond(*deltri, leftsubseg);
9004 tspivot(righttri, rightsubseg);
9005 if (rightsubseg.ss !=
m->dummysub) {
9006 tsbond(deltriright, rightsubseg);
9010 org(lefttri, neworg);
9013 testtriangle(
m,
b, deltri);
9036#ifdef ANSI_DECLARATORS
9039void undovertex(
m,
b)
9045 struct otri fliptri;
9046 struct otri botleft, botright, topright;
9047 struct otri botlcasing, botrcasing, toprcasing;
9048 struct otri gluetri;
9049 struct osub botlsubseg, botrsubseg, toprsubseg;
9050 vertex botvertex, rightvertex;
9058 decode(
m->lastflip->flippedtri, fliptri);
9064 if (
m->lastflip->prevflip == (
struct flipstacker *) NULL) {
9067 dprev(fliptri, botleft);
9069 onext(fliptri, botright);
9071 sym(botleft, botlcasing);
9072 sym(botright, botrcasing);
9073 dest(botleft, botvertex);
9077 bond(fliptri, botlcasing);
9079 tsbond(fliptri, botlsubseg);
9081 bond(fliptri, botrcasing);
9082 tspivot(botright, botrsubseg);
9083 tsbond(fliptri, botrsubseg);
9091 lprev(fliptri, gluetri);
9092 sym(gluetri, botright);
9094 sym(botright, botrcasing);
9095 dest(botright, rightvertex);
9097 setorg(fliptri, rightvertex);
9098 bond(gluetri, botrcasing);
9099 tspivot(botright, botrsubseg);
9100 tsbond(gluetri, botrsubseg);
9105 sym(fliptri, gluetri);
9106 if (gluetri.tri !=
m->dummytri) {
9108 dnext(gluetri, topright);
9109 sym(topright, toprcasing);
9111 setorg(gluetri, rightvertex);
9112 bond(gluetri, toprcasing);
9113 tspivot(topright, toprsubseg);
9114 tsbond(gluetri, toprsubseg);
9128 m->lastflip =
m->lastflip->prevflip;
9183#ifdef ANSI_DECLARATORS
9194 REAL pivotx, pivoty;
9197 if (arraysize == 2) {
9199 if ((sortarray[0][0] > sortarray[1][0]) ||
9200 ((sortarray[0][0] == sortarray[1][0]) &&
9201 (sortarray[0][1] > sortarray[1][1]))) {
9202 temp = sortarray[1];
9203 sortarray[1] = sortarray[0];
9204 sortarray[0] = temp;
9210 pivotx = sortarray[pivot][0];
9211 pivoty = sortarray[pivot][1];
9215 while (left < right) {
9219 }
while ((left <= right) && ((sortarray[left][0] < pivotx) ||
9220 ((sortarray[left][0] == pivotx) &&
9221 (sortarray[left][1] < pivoty))));
9225 }
while ((left <= right) && ((sortarray[right][0] > pivotx) ||
9226 ((sortarray[right][0] == pivotx) &&
9227 (sortarray[right][1] > pivoty))));
9230 temp = sortarray[left];
9231 sortarray[left] = sortarray[right];
9232 sortarray[right] = temp;
9239 if (right < arraysize - 2) {
9241 vertexsort(&sortarray[right + 1], arraysize - right - 1);
9257#ifdef ANSI_DECLARATORS
9270 REAL pivot1, pivot2;
9273 if (arraysize == 2) {
9275 if ((sortarray[0][axis] > sortarray[1][axis]) ||
9276 ((sortarray[0][axis] == sortarray[1][axis]) &&
9277 (sortarray[0][1 - axis] > sortarray[1][1 - axis]))) {
9278 temp = sortarray[1];
9279 sortarray[1] = sortarray[0];
9280 sortarray[0] = temp;
9286 pivot1 = sortarray[pivot][axis];
9287 pivot2 = sortarray[pivot][1 - axis];
9291 while (left < right) {
9295 }
while ((left <= right) && ((sortarray[left][axis] < pivot1) ||
9296 ((sortarray[left][axis] == pivot1) &&
9297 (sortarray[left][1 - axis] < pivot2))));
9301 }
while ((left <= right) && ((sortarray[right][axis] > pivot1) ||
9302 ((sortarray[right][axis] == pivot1) &&
9303 (sortarray[right][1 - axis] > pivot2))));
9306 temp = sortarray[left];
9307 sortarray[left] = sortarray[right];
9308 sortarray[right] = temp;
9313 if (left > median) {
9317 if (right < median - 1) {
9319 vertexmedian(&sortarray[right + 1], arraysize - right - 1,
9320 median - right - 1, axis);
9335#ifdef ANSI_DECLARATORS
9347 divider = arraysize >> 1;
9348 if (arraysize <= 3) {
9356 if (arraysize - divider >= 2) {
9360 alternateaxes(&sortarray[divider], arraysize - divider, 1 - axis);
9399#ifdef ANSI_DECLARATORS
9401 struct otri *innerleft,
struct otri *innerright,
9402 struct otri *farright,
int axis)
9404void mergehulls(
m,
b, farleft, innerleft, innerright, farright, axis)
9407struct otri *farleft;
9408struct otri *innerleft;
9409struct otri *innerright;
9410struct otri *farright;
9415 struct otri leftcand, rightcand;
9416 struct otri baseedge;
9417 struct otri nextedge;
9418 struct otri sidecasing, topcasing, outercasing;
9419 struct otri checkedge;
9422 vertex innerleftapex, innerrightapex;
9423 vertex farleftpt, farrightpt;
9424 vertex farleftapex, farrightapex;
9425 vertex lowerleft, lowerright;
9426 vertex upperleft, upperright;
9431 int leftfinished, rightfinished;
9434 dest(*innerleft, innerleftdest);
9435 apex(*innerleft, innerleftapex);
9436 org(*innerright, innerrightorg);
9437 apex(*innerright, innerrightapex);
9439 if (
b->dwyer && (axis == 1)) {
9440 org(*farleft, farleftpt);
9441 apex(*farleft, farleftapex);
9442 dest(*farright, farrightpt);
9443 apex(*farright, farrightapex);
9447 while (farleftapex[1] < farleftpt[1]) {
9450 farleftpt = farleftapex;
9451 apex(*farleft, farleftapex);
9453 sym(*innerleft, checkedge);
9454 apex(checkedge, checkvertex);
9455 while (checkvertex[1] > innerleftdest[1]) {
9456 lnext(checkedge, *innerleft);
9457 innerleftapex = innerleftdest;
9458 innerleftdest = checkvertex;
9459 sym(*innerleft, checkedge);
9460 apex(checkedge, checkvertex);
9462 while (innerrightapex[1] < innerrightorg[1]) {
9465 innerrightorg = innerrightapex;
9466 apex(*innerright, innerrightapex);
9468 sym(*farright, checkedge);
9469 apex(checkedge, checkvertex);
9470 while (checkvertex[1] > farrightpt[1]) {
9471 lnext(checkedge, *farright);
9472 farrightapex = farrightpt;
9473 farrightpt = checkvertex;
9474 sym(*farright, checkedge);
9475 apex(checkedge, checkvertex);
9486 innerleftdest = innerleftapex;
9487 apex(*innerleft, innerleftapex);
9495 innerrightorg = innerrightapex;
9496 apex(*innerright, innerrightapex);
9499 }
while (changemade);
9501 sym(*innerleft, leftcand);
9502 sym(*innerright, rightcand);
9506 bond(baseedge, *innerleft);
9508 bond(baseedge, *innerright);
9510 setorg(baseedge, innerrightorg);
9511 setdest(baseedge, innerleftdest);
9513 if (
b->verbose > 2) {
9514 printf(
" Creating base bounding ");
9518 org(*farleft, farleftpt);
9519 if (innerleftdest == farleftpt) {
9520 lnext(baseedge, *farleft);
9522 dest(*farright, farrightpt);
9523 if (innerrightorg == farrightpt) {
9524 lprev(baseedge, *farright);
9527 lowerleft = innerleftdest;
9528 lowerright = innerrightorg;
9530 apex(leftcand, upperleft);
9531 apex(rightcand, upperright);
9542 if (leftfinished && rightfinished) {
9545 setorg(nextedge, lowerleft);
9546 setdest(nextedge, lowerright);
9549 bond(nextedge, baseedge);
9551 bond(nextedge, rightcand);
9553 bond(nextedge, leftcand);
9554 if (
b->verbose > 2) {
9555 printf(
" Creating top bounding ");
9559 if (
b->dwyer && (axis == 1)) {
9560 org(*farleft, farleftpt);
9561 apex(*farleft, farleftapex);
9562 dest(*farright, farrightpt);
9563 apex(*farright, farrightapex);
9564 sym(*farleft, checkedge);
9565 apex(checkedge, checkvertex);
9569 while (checkvertex[0] < farleftpt[0]) {
9570 lprev(checkedge, *farleft);
9571 farleftapex = farleftpt;
9572 farleftpt = checkvertex;
9573 sym(*farleft, checkedge);
9574 apex(checkedge, checkvertex);
9576 while (farrightapex[0] > farrightpt[0]) {
9579 farrightpt = farrightapex;
9580 apex(*farright, farrightapex);
9586 if (!leftfinished) {
9588 lprev(leftcand, nextedge);
9590 apex(nextedge, nextapex);
9593 if (nextapex != (
vertex) NULL) {
9595 badedge =
incircle(
m,
b, lowerleft, lowerright, upperleft, nextapex) >
9601 sym(nextedge, topcasing);
9603 sym(nextedge, sidecasing);
9604 bond(nextedge, topcasing);
9605 bond(leftcand, sidecasing);
9607 sym(leftcand, outercasing);
9609 bond(nextedge, outercasing);
9611 setorg(leftcand, lowerleft);
9618 upperleft = nextapex;
9621 apex(nextedge, nextapex);
9622 if (nextapex != (
vertex) NULL) {
9624 badedge =
incircle(
m,
b, lowerleft, lowerright, upperleft,
9634 if (!rightfinished) {
9636 lnext(rightcand, nextedge);
9638 apex(nextedge, nextapex);
9641 if (nextapex != (
vertex) NULL) {
9643 badedge =
incircle(
m,
b, lowerleft, lowerright, upperright, nextapex) >
9649 sym(nextedge, topcasing);
9651 sym(nextedge, sidecasing);
9652 bond(nextedge, topcasing);
9653 bond(rightcand, sidecasing);
9655 sym(rightcand, outercasing);
9657 bond(nextedge, outercasing);
9660 setdest(rightcand, lowerright);
9662 setorg(nextedge, upperright);
9666 upperright = nextapex;
9669 apex(nextedge, nextapex);
9670 if (nextapex != (
vertex) NULL) {
9672 badedge =
incircle(
m,
b, lowerleft, lowerright, upperright,
9681 if (leftfinished || (!rightfinished &&
9682 (
incircle(
m,
b, upperleft, lowerleft, lowerright, upperright) >
9686 bond(baseedge, rightcand);
9687 lprev(rightcand, baseedge);
9689 lowerright = upperright;
9690 sym(baseedge, rightcand);
9691 apex(rightcand, upperright);
9695 bond(baseedge, leftcand);
9696 lnext(leftcand, baseedge);
9697 setorg(baseedge, lowerright);
9698 lowerleft = upperleft;
9699 sym(baseedge, leftcand);
9700 apex(leftcand, upperleft);
9702 if (
b->verbose > 2) {
9703 printf(
" Connecting ");
9726#ifdef ANSI_DECLARATORS
9728 int vertices,
int axis,
9729 struct otri *farleft,
struct otri *farright)
9737struct otri *farleft;
9738struct otri *farright;
9742 struct otri midtri, tri1, tri2, tri3;
9743 struct otri innerleft, innerright;
9747 if (
b->verbose > 2) {
9748 printf(
" Triangulating %d vertices.\n", vertices);
9750 if (vertices == 2) {
9754 setorg(*farleft, sortarray[0]);
9755 setdest(*farleft, sortarray[1]);
9758 setorg(*farright, sortarray[1]);
9759 setdest(*farright, sortarray[0]);
9761 bond(*farleft, *farright);
9764 bond(*farleft, *farright);
9767 bond(*farleft, *farright);
9768 if (
b->verbose > 2) {
9769 printf(
" Creating ");
9771 printf(
" Creating ");
9775 lprev(*farright, *farleft);
9777 }
else if (vertices == 3) {
9788 setorg(midtri, sortarray[0]);
9789 setdest(midtri, sortarray[1]);
9790 setorg(tri1, sortarray[1]);
9792 setorg(tri2, sortarray[2]);
9794 setorg(tri3, sortarray[1]);
9818 setorg(midtri, sortarray[0]);
9820 setorg(tri3, sortarray[0]);
9824 setdest(midtri, sortarray[1]);
9825 setorg(tri1, sortarray[1]);
9827 setapex(midtri, sortarray[2]);
9828 setorg(tri2, sortarray[2]);
9832 setdest(midtri, sortarray[2]);
9833 setorg(tri1, sortarray[2]);
9835 setapex(midtri, sortarray[1]);
9836 setorg(tri2, sortarray[1]);
9860 lnext(*farleft, *farright);
9863 if (
b->verbose > 2) {
9864 printf(
" Creating ");
9866 printf(
" Creating ");
9868 printf(
" Creating ");
9870 printf(
" Creating ");
9876 divider = vertices >> 1;
9880 &innerright, farright);
9881 if (
b->verbose > 1) {
9882 printf(
" Joining triangulations with %d and %d vertices.\n", divider,
9883 vertices - divider);
9886 mergehulls(
m,
b, farleft, &innerleft, &innerright, farright, axis);
9890#ifdef ANSI_DECLARATORS
9896struct otri *startghost;
9900 struct otri searchedge;
9901 struct otri dissolveedge;
9902 struct otri deadtriangle;
9908 printf(
" Removing ghost triangles.\n");
9911 lprev(*startghost, searchedge);
9913 m->dummytri[0] =
encode(searchedge);
9915 otricopy(*startghost, dissolveedge);
9919 lnext(dissolveedge, deadtriangle);
9926 if (dissolveedge.
tri !=
m->dummytri) {
9927 org(dissolveedge, markorg);
9936 sym(deadtriangle, dissolveedge);
9939 }
while (!
otriequal(dissolveedge, *startghost));
9953#ifdef ANSI_DECLARATORS
9963 struct otri hullleft, hullright;
9968 printf(
" Sorting vertices.\n");
9974 for (i = 0; i <
m->invertices; i++) {
9981 for (j = 1; j <
m->invertices; j++) {
9982 if ((sortarray[i][0] == sortarray[j][0])
9983 && (sortarray[i][1] == sortarray[j][1])) {
9986"Warning: A duplicate vertex at (%.12g, %.12g) appeared and was ignored.\n",
9987 sortarray[j][0], sortarray[j][1]);
9993 sortarray[i] = sortarray[j];
10000 if (i - divider >= 2) {
10001 if (divider >= 2) {
10009 printf(
" Forming triangulation.\n");
10040#ifdef ANSI_DECLARATORS
10043void boundingbox(
m,
b)
10049 struct otri inftri;
10053 printf(
" Creating triangular bounding box.\n");
10057 if (
m->ymax -
m->ymin >
width) {
10060 if (
width == 0.0) {
10067 m->infvertex1[0] =
m->xmin - 50.0 *
width;
10068 m->infvertex1[1] =
m->ymin - 40.0 *
width;
10069 m->infvertex2[0] =
m->xmax + 50.0 *
width;
10070 m->infvertex2[1] =
m->ymin - 40.0 *
width;
10071 m->infvertex3[0] = 0.5 * (
m->xmin +
m->xmax);
10072 m->infvertex3[1] =
m->ymax + 60.0 *
width;
10076 setorg(inftri,
m->infvertex1);
10081 m->dummytri[0] = (
triangle) inftri.tri;
10082 if (
b->verbose > 2) {
10083 printf(
" Creating ");
10106#ifdef ANSI_DECLARATORS
10109long removebox(
m,
b)
10115 struct otri deadtriangle;
10116 struct otri searchedge;
10117 struct otri checkedge;
10118 struct otri nextedge, finaledge, dissolveedge;
10124 printf(
" Removing triangular bounding box.\n");
10127 nextedge.tri =
m->dummytri;
10128 nextedge.orient = 0;
10131 lprev(nextedge, finaledge);
10136 lprev(nextedge, searchedge);
10140 lnext(nextedge, checkedge);
10142 if (checkedge.tri ==
m->dummytri) {
10151 m->dummytri[0] =
encode(searchedge);
10153 while (!
otriequal(nextedge, finaledge)) {
10155 lprev(nextedge, dissolveedge);
10164 if (dissolveedge.tri !=
m->dummytri) {
10165 org(dissolveedge, markorg);
10173 lnext(nextedge, deadtriangle);
10174 sym(deadtriangle, nextedge);
10178 if (nextedge.tri ==
m->dummytri) {
10205#ifdef ANSI_DECLARATORS
10208long incrementaldelaunay(
m,
b)
10214 struct otri starttri;
10220 printf(
" Incrementally inserting vertices.\n");
10224 while (vertexloop != (
vertex) NULL) {
10225 starttri.tri =
m->dummytri;
10230"Warning: A duplicate vertex at (%.12g, %.12g) appeared and was ignored.\n",
10231 vertexloop[0], vertexloop[1]);
10239 return removebox(
m,
b);
10254#ifdef ANSI_DECLARATORS
10255void eventheapinsert(
struct event **heap,
int heapsize,
struct event *newevent)
10257void eventheapinsert(heap, heapsize, newevent)
10258struct event **heap;
10260struct event *newevent;
10264 REAL eventx, eventy;
10269 eventx = newevent->
xkey;
10270 eventy = newevent->
ykey;
10271 eventnum = heapsize;
10272 notdone = eventnum > 0;
10274 parent = (eventnum - 1) >> 1;
10275 if ((heap[parent]->
ykey < eventy) ||
10276 ((heap[parent]->
ykey == eventy)
10277 && (heap[parent]->
xkey <= eventx))) {
10280 heap[eventnum] = heap[parent];
10284 notdone = eventnum > 0;
10287 heap[eventnum] = newevent;
10295#ifdef ANSI_DECLARATORS
10296void eventheapify(
struct event **heap,
int heapsize,
int eventnum)
10298void eventheapify(heap, heapsize, eventnum)
10299struct event **heap;
10305 struct event *thisevent;
10306 REAL eventx, eventy;
10307 int leftchild, rightchild;
10311 thisevent = heap[eventnum];
10312 eventx = thisevent->
xkey;
10313 eventy = thisevent->
ykey;
10314 leftchild = 2 * eventnum + 1;
10315 notdone = leftchild < heapsize;
10317 if ((heap[leftchild]->
ykey < eventy) ||
10318 ((heap[leftchild]->
ykey == eventy)
10319 && (heap[leftchild]->
xkey < eventx))) {
10320 smallest = leftchild;
10322 smallest = eventnum;
10324 rightchild = leftchild + 1;
10325 if (rightchild < heapsize) {
10326 if ((heap[rightchild]->
ykey < heap[smallest]->
ykey) ||
10327 ((heap[rightchild]->
ykey == heap[smallest]->
ykey)
10328 && (heap[rightchild]->
xkey < heap[smallest]->
xkey))) {
10329 smallest = rightchild;
10332 if (smallest == eventnum) {
10335 heap[eventnum] = heap[smallest];
10337 heap[smallest] = thisevent;
10340 eventnum = smallest;
10341 leftchild = 2 * eventnum + 1;
10342 notdone = leftchild < heapsize;
10351#ifdef ANSI_DECLARATORS
10352void eventheapdelete(
struct event **heap,
int heapsize,
int eventnum)
10354void eventheapdelete(heap, heapsize, eventnum)
10355struct event **heap;
10361 struct event *moveevent;
10362 REAL eventx, eventy;
10366 moveevent = heap[heapsize - 1];
10367 if (eventnum > 0) {
10368 eventx = moveevent->
xkey;
10369 eventy = moveevent->
ykey;
10371 parent = (eventnum - 1) >> 1;
10372 if ((heap[parent]->
ykey < eventy) ||
10373 ((heap[parent]->
ykey == eventy)
10374 && (heap[parent]->
xkey <= eventx))) {
10377 heap[eventnum] = heap[parent];
10381 notdone = eventnum > 0;
10385 heap[eventnum] = moveevent;
10387 eventheapify(heap, heapsize - 1, eventnum);
10394#ifdef ANSI_DECLARATORS
10395void createeventheap(
struct mesh *
m,
struct event ***eventheap,
10396 struct event **events,
struct event **freeevents)
10398void createeventheap(
m, eventheap, events, freeevents)
10400struct event ***eventheap;
10401struct event **events;
10402struct event **freeevents;
10410 maxevents = (3 *
m->invertices) / 2;
10412 (
int)
sizeof(
struct event *));
10415 for (i = 0; i <
m->invertices; i++) {
10417 (*events)[i].eventptr = (
VOID *) thisvertex;
10418 (*events)[i].xkey = thisvertex[0];
10419 (*events)[i].ykey = thisvertex[1];
10420 eventheapinsert(*eventheap, i, *events + i);
10422 *freeevents = (
struct event *) NULL;
10423 for (i = maxevents - 1; i >=
m->invertices; i--) {
10425 *freeevents = *events + i;
10433#ifdef ANSI_DECLARATORS
10434int rightofhyperbola(
struct mesh *
m,
struct otri *fronttri,
vertex newsite)
10436int rightofhyperbola(
m, fronttri, newsite)
10438struct otri *fronttri;
10443 vertex leftvertex, rightvertex;
10444 REAL dxa, dya, dxb, dyb;
10446 m->hyperbolacount++;
10448 dest(*fronttri, leftvertex);
10449 apex(*fronttri, rightvertex);
10450 if ((leftvertex[1] < rightvertex[1]) ||
10451 ((leftvertex[1] == rightvertex[1]) &&
10452 (leftvertex[0] < rightvertex[0]))) {
10453 if (newsite[0] >= rightvertex[0]) {
10457 if (newsite[0] <= leftvertex[0]) {
10461 dxa = leftvertex[0] - newsite[0];
10462 dya = leftvertex[1] - newsite[1];
10463 dxb = rightvertex[0] - newsite[0];
10464 dyb = rightvertex[1] - newsite[1];
10465 return dya * (dxb * dxb + dyb * dyb) > dyb * (dxa * dxa + dya * dya);
10472#ifdef ANSI_DECLARATORS
10475REAL circletop(
m, pa, pb, pc, ccwabc)
10484 REAL xac, yac, xbc, ybc, xab, yab;
10485 REAL acLen2, bcLen2, ablen2;
10487 m->circletopcount++;
10489 xac = pa[0] - pc[0];
10490 yac = pa[1] - pc[1];
10491 xbc = pb[0] - pc[0];
10492 ybc = pb[1] - pc[1];
10493 xab = pa[0] - pb[0];
10494 yab = pa[1] - pb[1];
10495 acLen2 = xac * xac + yac * yac;
10496 bcLen2 = xbc * xbc + ybc * ybc;
10497 ablen2 = xab * xab + yab * yab;
10498 return pc[1] + (xac * bcLen2 - xbc * acLen2 +
sqrt(acLen2 * bcLen2 * ablen2))
10506#ifdef ANSI_DECLARATORS
10507void check4deadevent(
struct otri *checktri,
struct event **freeevents,
10508 struct event **eventheap,
int *heapsize)
10510void check4deadevent(checktri, freeevents, eventheap, heapsize)
10511struct otri *checktri;
10512struct event **freeevents;
10513struct event **eventheap;
10518 struct event *deadevent;
10522 org(*checktri, eventvertex);
10523 if (eventvertex != (
vertex) NULL) {
10524 deadevent = (
struct event *) eventvertex;
10527 *freeevents = deadevent;
10528 eventheapdelete(eventheap, *heapsize, eventnum);
10530 setorg(*checktri, NULL);
10538#ifdef ANSI_DECLARATORS
10542struct splaynode *splay(
m, splaytree, searchpoint, searchtri)
10546struct otri *searchtri;
10551 struct splaynode *lefttree, *righttree;
10554 int rightofroot, rightofchild;
10556 if (splaytree == (
struct splaynode *) NULL) {
10560 if (checkvertex == splaytree->
keydest) {
10561 rightofroot = rightofhyperbola(
m, &splaytree->
keyedge, searchpoint);
10572 if (checkvertex !=
child->keydest) {
10583 rightofchild = rightofhyperbola(
m, &
child->keyedge, searchpoint);
10584 if (rightofchild) {
10586 grandchild = splay(
m,
child->rchild, searchpoint, searchtri);
10587 child->rchild = grandchild;
10589 grandchild = splay(
m,
child->lchild, searchpoint, searchtri);
10590 child->lchild = grandchild;
10592 if (grandchild == (
struct splaynode *) NULL) {
10595 child->lchild = splaytree;
10598 child->rchild = splaytree;
10602 if (rightofchild) {
10605 child->lchild = splaytree;
10608 grandchild->
rchild = splaytree;
10615 grandchild->
lchild = splaytree;
10618 child->rchild = splaytree;
10625 lefttree = splay(
m, splaytree->
lchild, searchpoint, searchtri);
10626 righttree = splay(
m, splaytree->
rchild, searchpoint, searchtri);
10629 if (lefttree == (
struct splaynode *) NULL) {
10631 }
else if (righttree == (
struct splaynode *) NULL) {
10635 righttree->
lchild = lefttree;
10639 lefttree->
rchild = righttree;
10643 leftright = lefttree->
rchild;
10645 leftright = leftright->
rchild;
10647 leftright->
rchild = righttree;
10657#ifdef ANSI_DECLARATORS
10661struct splaynode *splayinsert(
m, splayroot, newkey, searchpoint)
10664struct otri *newkey;
10674 if (splayroot == (
struct splaynode *) NULL) {
10677 }
else if (rightofhyperbola(
m, &splayroot->
keyedge, searchpoint)) {
10678 newsplaynode->
lchild = splayroot;
10683 newsplaynode->
rchild = splayroot;
10686 return newsplaynode;
10693#ifdef ANSI_DECLARATORS
10696 struct otri *newkey,
10699struct splaynode *circletopinsert(
m,
b, splayroot, newkey, pa, pb, pc, topy)
10703struct otri *newkey;
10712 REAL xac, yac, xbc, ybc;
10713 REAL acLen2, bcLen2;
10714 REAL searchpoint[2];
10715 struct otri dummytri;
10718 xac = pa[0] - pc[0];
10719 yac = pa[1] - pc[1];
10720 xbc = pb[0] - pc[0];
10721 ybc = pb[1] - pc[1];
10722 acLen2 = xac * xac + yac * yac;
10723 bcLen2 = xbc * xbc + ybc * ybc;
10724 searchpoint[0] = pc[0] - (yac * bcLen2 - ybc * acLen2) / (2.0 * ccwabc);
10725 searchpoint[1] = topy;
10726 return splayinsert(
m, splay(
m, splayroot, (
vertex) searchpoint, &dummytri),
10727 newkey, (
vertex) searchpoint);
10734#ifdef ANSI_DECLARATORS
10736 struct otri *bottommost,
vertex searchvertex,
10737 struct otri *searchtri,
int *farright)
10739struct splaynode *frontlocate(
m, splayroot, bottommost, searchvertex,
10740 searchtri, farright)
10743struct otri *bottommost;
10745struct otri *searchtri;
10753 otricopy(*bottommost, *searchtri);
10754 splayroot = splay(
m, splayroot, searchvertex, searchtri);
10757 while (!farrightflag && rightofhyperbola(
m, searchtri, searchvertex)) {
10759 farrightflag =
otriequal(*searchtri, *bottommost);
10761 *farright = farrightflag;
10769#ifdef ANSI_DECLARATORS
10772long sweeplinedelaunay(
m,
b)
10778 struct event **eventheap;
10779 struct event *events;
10780 struct event *freeevents;
10781 struct event *nextevent;
10782 struct event *newevent;
10784 struct otri bottommost;
10785 struct otri searchtri;
10786 struct otri fliptri;
10787 struct otri lefttri, righttri, farlefttri, farrighttri;
10788 struct otri inserttri;
10789 vertex firstvertex, secondvertex;
10790 vertex nextvertex, lastvertex;
10792 vertex leftvertex, midvertex, rightvertex;
10793 REAL lefttest, righttest;
10795 int check4events, farrightflag;
10803 printf(
" Placing vertices in event heap.\n");
10805 createeventheap(
m, &eventheap, &events, &freeevents);
10806 heapsize =
m->invertices;
10809 printf(
" Forming triangulation.\n");
10813 bond(lefttri, righttri);
10816 bond(lefttri, righttri);
10819 bond(lefttri, righttri);
10820 firstvertex = (
vertex) eventheap[0]->eventptr;
10822 freeevents = eventheap[0];
10823 eventheapdelete(eventheap, heapsize, 0);
10826 if (heapsize == 0) {
10827 printf(
"Error: Input vertices are all identical.\n");
10830 secondvertex = (
vertex) eventheap[0]->eventptr;
10832 freeevents = eventheap[0];
10833 eventheapdelete(eventheap, heapsize, 0);
10835 if ((firstvertex[0] == secondvertex[0]) &&
10836 (firstvertex[1] == secondvertex[1])) {
10839"Warning: A duplicate vertex at (%.12g, %.12g) appeared and was ignored.\n",
10840 secondvertex[0], secondvertex[1]);
10845 }
while ((firstvertex[0] == secondvertex[0]) &&
10846 (firstvertex[1] == secondvertex[1]));
10847 setorg(lefttri, firstvertex);
10848 setdest(lefttri, secondvertex);
10849 setorg(righttri, secondvertex);
10850 setdest(righttri, firstvertex);
10851 lprev(lefttri, bottommost);
10852 lastvertex = secondvertex;
10853 while (heapsize > 0) {
10854 nextevent = eventheap[0];
10855 eventheapdelete(eventheap, heapsize, 0);
10858 if (nextevent->
xkey <
m->xmin) {
10860 oprev(fliptri, farlefttri);
10861 check4deadevent(&farlefttri, &freeevents, eventheap, &heapsize);
10862 onext(fliptri, farrighttri);
10863 check4deadevent(&farrighttri, &freeevents, eventheap, &heapsize);
10865 if (
otriequal(farlefttri, bottommost)) {
10866 lprev(fliptri, bottommost);
10870 lprev(fliptri, lefttri);
10871 lnext(fliptri, righttri);
10872 sym(lefttri, farlefttri);
10876 dest(fliptri, leftvertex);
10877 apex(fliptri, midvertex);
10878 org(fliptri, rightvertex);
10879 splayroot = circletopinsert(
m,
b, splayroot, &lefttri, leftvertex,
10880 midvertex, rightvertex, nextevent->
ykey);
10884 if ((nextvertex[0] == lastvertex[0]) &&
10885 (nextvertex[1] == lastvertex[1])) {
10888"Warning: A duplicate vertex at (%.12g, %.12g) appeared and was ignored.\n",
10889 nextvertex[0], nextvertex[1]);
10895 lastvertex = nextvertex;
10897 splayroot = frontlocate(
m, splayroot, &bottommost, nextvertex,
10898 &searchtri, &farrightflag);
10908 check4deadevent(&searchtri, &freeevents, eventheap, &heapsize);
10911 sym(searchtri, farlefttri);
10914 dest(farrighttri, connectvertex);
10915 setorg(lefttri, connectvertex);
10916 setdest(lefttri, nextvertex);
10917 setorg(righttri, nextvertex);
10918 setdest(righttri, connectvertex);
10919 bond(lefttri, righttri);
10922 bond(lefttri, righttri);
10925 bond(lefttri, farlefttri);
10926 bond(righttri, farrighttri);
10927 if (!farrightflag &&
otriequal(farrighttri, bottommost)) {
10932 splayroot = splayinsert(
m, splayroot, &lefttri, nextvertex);
10934 lnext(righttri, inserttri);
10935 splayroot = splayinsert(
m, splayroot, &inserttri, nextvertex);
10940 freeevents = nextevent;
10942 if (check4events) {
10943 apex(farlefttri, leftvertex);
10944 dest(lefttri, midvertex);
10945 apex(lefttri, rightvertex);
10947 if (lefttest > 0.0) {
10948 newevent = freeevents;
10950 newevent->
xkey =
m->xminextreme;
10951 newevent->
ykey = circletop(
m, leftvertex, midvertex, rightvertex,
10954 eventheapinsert(eventheap, heapsize, newevent);
10956 setorg(lefttri, newevent);
10958 apex(righttri, leftvertex);
10959 org(righttri, midvertex);
10960 apex(farrighttri, rightvertex);
10962 if (righttest > 0.0) {
10963 newevent = freeevents;
10965 newevent->
xkey =
m->xminextreme;
10966 newevent->
ykey = circletop(
m, leftvertex, midvertex, rightvertex,
10969 eventheapinsert(eventheap, heapsize, newevent);
10971 setorg(farrighttri, newevent);
10997#ifdef ANSI_DECLARATORS
11014 "Constructing Delaunay triangulation by divide-and-conquer method.\n");
11019 printf(
"Constructing Delaunay triangulation ");
11020 if (
b->incremental) {
11021 printf(
"by incremental method.\n");
11022 }
else if (
b->sweepline) {
11023 printf(
"by sweepline method.\n");
11025 printf(
"by divide-and-conquer method.\n");
11028 if (
b->incremental) {
11029 hulledges = incrementaldelaunay(
m,
b);
11030 }
else if (
b->sweepline) {
11031 hulledges = sweeplinedelaunay(
m,
b);
11037 if (
m->triangles.items == 0) {
11074#ifdef ANSI_DECLARATORS
11075int reconstruct(
struct mesh *
m,
struct behavior *
b,
int *trianglelist,
11076 REAL *triangleattriblist,
REAL *trianglearealist,
11077 int elements,
int corners,
int attribs,
11078 int *segmentlist,
int *segmentmarkerlist,
int numberofsegments)
11080int reconstruct(
m,
b, trianglelist, triangleattriblist, trianglearealist,
11081 elements, corners, attribs, segmentlist, segmentmarkerlist,
11086REAL *triangleattriblist;
11087REAL *trianglearealist;
11092int *segmentmarkerlist;
11093int numberofsegments;
11098#ifdef ANSI_DECLARATORS
11099long reconstruct(
struct mesh *
m,
struct behavior *
b,
char *elefilename,
11100 char *areafilename,
char *polyfilename, FILE *polyfile)
11102long reconstruct(
m,
b, elefilename, areafilename, polyfilename, polyfile)
11124 struct otri triangleloop;
11125 struct otri triangleleft;
11126 struct otri checktri;
11127 struct otri checkleft;
11128 struct otri checkneighbor;
11129 struct osub subsegloop;
11134 vertex checkdest, checkapex;
11137 vertex segmentorg, segmentdest;
11141 int killvertexindex;
11143 int segmentmarkers;
11148 long elementnumber, segmentnumber;
11153 m->inelements = elements;
11154 incorners = corners;
11155 if (incorners < 3) {
11156 printf(
"Error: Triangles must have at least 3 vertices.\n");
11159 m->eextras = attribs;
11163 printf(
"Opening %s.\n", elefilename);
11165 elefile = fopen(elefilename,
"r");
11166 if (elefile == (FILE *) NULL) {
11167 printf(
" Error: Cannot access file %s.\n", elefilename);
11172 stringptr =
readline(inputline, elefile, elefilename);
11173 m->inelements = (
int) strtol(stringptr, &stringptr, 0);
11174 stringptr = findfield(stringptr);
11175 if (*stringptr ==
'\0') {
11178 incorners = (
int) strtol(stringptr, &stringptr, 0);
11179 if (incorners < 3) {
11180 printf(
"Error: Triangles in %s must have at least 3 vertices.\n",
11185 stringptr = findfield(stringptr);
11186 if (*stringptr ==
'\0') {
11189 m->eextras = (
int) strtol(stringptr, &stringptr, 0);
11196 for (elementnumber = 1; elementnumber <=
m->inelements; elementnumber++) {
11199 triangleloop.tri[3] = (
triangle) triangleloop.tri;
11202 segmentmarkers = 0;
11205 m->insegments = numberofsegments;
11206 segmentmarkers = segmentmarkerlist != (
int *) NULL;
11210 stringptr =
readline(inputline, polyfile,
b->inpolyfilename);
11211 m->insegments = (
int) strtol(stringptr, &stringptr, 0);
11212 stringptr = findfield(stringptr);
11213 if (*stringptr !=
'\0') {
11214 segmentmarkers = (
int) strtol(stringptr, &stringptr, 0);
11219 for (segmentnumber = 1; segmentnumber <=
m->insegments; segmentnumber++) {
11222 subsegloop.ss[2] = (
subseg) subsegloop.ss;
11233 printf(
"Opening %s.\n", areafilename);
11235 areafile = fopen(areafilename,
"r");
11236 if (areafile == (FILE *) NULL) {
11237 printf(
" Error: Cannot access file %s.\n", areafilename);
11240 stringptr =
readline(inputline, areafile, areafilename);
11241 areaelements = (
int) strtol(stringptr, &stringptr, 0);
11242 if (areaelements !=
m->inelements) {
11243 printf(
"Error: %s and %s disagree on number of triangles.\n",
11244 elefilename, areafilename);
11251 printf(
"Reconstructing mesh.\n");
11259 for (i = 0; i <
m->vertices.items; i++) {
11260 vertexarray[i] = (
triangle)
m->dummytri;
11264 printf(
" Assembling triangles.\n");
11270 elementnumber =
b->firstnumber;
11271 while (triangleloop.tri != (
triangle *) NULL) {
11274 for (j = 0; j < 3; j++) {
11275 corner[j] = trianglelist[vertexindex++];
11276 if ((corner[j] <
b->firstnumber) ||
11277 (corner[j] >=
b->firstnumber +
m->invertices)) {
11278 printf(
"Error: Triangle %ld has an invalid vertex index.\n",
11285 stringptr =
readline(inputline, elefile, elefilename);
11286 for (j = 0; j < 3; j++) {
11287 stringptr = findfield(stringptr);
11288 if (*stringptr ==
'\0') {
11289 printf(
"Error: Triangle %ld is missing vertex %d in %s.\n",
11290 elementnumber, j + 1, elefilename);
11293 corner[j] = (
int) strtol(stringptr, &stringptr, 0);
11294 if ((corner[j] <
b->firstnumber) ||
11295 (corner[j] >=
b->firstnumber +
m->invertices)) {
11296 printf(
"Error: Triangle %ld has an invalid vertex index.\n",
11305 for (j = 3; j < incorners; j++) {
11307 killvertexindex = trianglelist[vertexindex++];
11309 stringptr = findfield(stringptr);
11310 if (*stringptr !=
'\0') {
11311 killvertexindex = (
int) strtol(stringptr, &stringptr, 0);
11313 if ((killvertexindex >=
b->firstnumber) &&
11314 (killvertexindex <
b->firstnumber +
m->invertices)) {
11327 for (j = 0; j <
m->eextras; j++) {
11331 stringptr = findfield(stringptr);
11332 if (*stringptr ==
'\0') {
11336 (
REAL) strtod(stringptr, &stringptr));
11343 area = trianglearealist[elementnumber -
b->firstnumber];
11346 stringptr =
readline(inputline, areafile, areafilename);
11347 stringptr = findfield(stringptr);
11348 if (*stringptr ==
'\0') {
11351 area = (
REAL) strtod(stringptr, &stringptr);
11358 triangleloop.orient = 0;
11363 for (triangleloop.orient = 0; triangleloop.orient < 3;
11364 triangleloop.orient++) {
11366 aroundvertex = corner[triangleloop.orient];
11368 nexttri = vertexarray[aroundvertex -
b->firstnumber];
11370 triangleloop.tri[6 + triangleloop.orient] = nexttri;
11372 vertexarray[aroundvertex -
b->firstnumber] =
encode(triangleloop);
11373 decode(nexttri, checktri);
11374 if (checktri.
tri !=
m->dummytri) {
11375 dest(triangleloop, tdest);
11376 apex(triangleloop, tapex);
11379 dest(checktri, checkdest);
11380 apex(checktri, checkapex);
11381 if (tapex == checkdest) {
11383 lprev(triangleloop, triangleleft);
11384 bond(triangleleft, checktri);
11386 if (tdest == checkapex) {
11388 lprev(checktri, checkleft);
11389 bond(triangleloop, checkleft);
11392 nexttri = checktri.
tri[6 + checktri.
orient];
11393 decode(nexttri, checktri);
11394 }
while (checktri.
tri !=
m->dummytri);
11413 printf(
" Marking segments in triangulation.\n");
11420 segmentnumber =
b->firstnumber;
11421 while (subsegloop.ss != (
subseg *) NULL) {
11423 end[0] = segmentlist[vertexindex++];
11424 end[1] = segmentlist[vertexindex++];
11425 if (segmentmarkers) {
11426 boundmarker = segmentmarkerlist[segmentnumber -
b->firstnumber];
11430 stringptr =
readline(inputline, polyfile,
b->inpolyfilename);
11432 stringptr = findfield(stringptr);
11433 if (*stringptr ==
'\0') {
11434 printf(
"Error: Segment %ld has no endpoints in %s.\n", segmentnumber,
11438 end[0] = (
int) strtol(stringptr, &stringptr, 0);
11440 stringptr = findfield(stringptr);
11441 if (*stringptr ==
'\0') {
11442 printf(
"Error: Segment %ld is missing its second endpoint in %s.\n",
11443 segmentnumber, polyfilename);
11446 end[1] = (
int) strtol(stringptr, &stringptr, 0);
11448 if (segmentmarkers) {
11449 stringptr = findfield(stringptr);
11450 if (*stringptr ==
'\0') {
11453 boundmarker = (
int) strtol(stringptr, &stringptr, 0);
11457 for (j = 0; j < 2; j++) {
11458 if ((end[j] <
b->firstnumber) ||
11459 (
end[j] >=
b->firstnumber +
m->invertices)) {
11460 printf(
"Error: Segment %ld has an invalid vertex index.\n",
11467 subsegloop.ssorient = 0;
11470 setsorg(subsegloop, segmentorg);
11471 setsdest(subsegloop, segmentdest);
11474 setmark(subsegloop, boundmarker);
11476 for (subsegloop.ssorient = 0; subsegloop.ssorient < 2;
11477 subsegloop.ssorient++) {
11479 aroundvertex =
end[1 - subsegloop.ssorient];
11481 prevlink = &vertexarray[aroundvertex -
b->firstnumber];
11482 nexttri = vertexarray[aroundvertex -
b->firstnumber];
11483 decode(nexttri, checktri);
11484 sorg(subsegloop, shorg);
11493 while (notfound && (checktri.
tri !=
m->dummytri)) {
11494 dest(checktri, checkdest);
11495 if (shorg == checkdest) {
11497 *prevlink = checktri.
tri[6 + checktri.
orient];
11499 tsbond(checktri, subsegloop);
11501 sym(checktri, checkneighbor);
11502 if (checkneighbor.tri ==
m->dummytri) {
11512 prevlink = &checktri.
tri[6 + checktri.
orient];
11513 nexttri = checktri.
tri[6 + checktri.
orient];
11514 decode(nexttri, checktri);
11524 for (i = 0; i <
m->vertices.items; i++) {
11526 nexttri = vertexarray[i];
11527 decode(nexttri, checktri);
11528 while (checktri.
tri !=
m->dummytri) {
11531 nexttri = checktri.
tri[6 + checktri.
orient];
11534 sym(checktri, checkneighbor);
11535 if (checkneighbor.tri ==
m->dummytri) {
11539 decode(nexttri, checktri);
11574#ifdef ANSI_DECLARATORS
11576 struct otri *searchtri,
11582struct otri *searchtri;
11587 struct otri checktri;
11589 vertex leftvertex, rightvertex;
11590 REAL leftccw, rightccw;
11591 int leftflag, rightflag;
11594 org(*searchtri, startvertex);
11595 dest(*searchtri, rightvertex);
11596 apex(*searchtri, leftvertex);
11599 leftflag = leftccw > 0.0;
11602 rightflag = rightccw > 0.0;
11603 if (leftflag && rightflag) {
11606 onext(*searchtri, checktri);
11607 if (checktri.
tri ==
m->dummytri) {
11616 if (searchtri->
tri ==
m->dummytri) {
11617 printf(
"Internal error in finddirection(): Unable to find a\n");
11618 printf(
" triangle leading from (%.12g, %.12g) to", startvertex[0],
11620 printf(
" (%.12g, %.12g).\n", searchpoint[0], searchpoint[1]);
11623 apex(*searchtri, leftvertex);
11624 rightccw = leftccw;
11626 leftflag = leftccw > 0.0;
11628 while (rightflag) {
11631 if (searchtri->
tri ==
m->dummytri) {
11632 printf(
"Internal error in finddirection(): Unable to find a\n");
11633 printf(
" triangle leading from (%.12g, %.12g) to", startvertex[0],
11635 printf(
" (%.12g, %.12g).\n", searchpoint[0], searchpoint[1]);
11638 dest(*searchtri, rightvertex);
11639 leftccw = rightccw;
11641 rightflag = rightccw > 0.0;
11643 if (leftccw == 0.0) {
11645 }
else if (rightccw == 0.0) {
11669#ifdef ANSI_DECLARATORS
11671 struct otri *splittri,
struct osub *splitsubseg,
11677struct otri *splittri;
11678struct osub *splitsubseg;
11683 struct osub opposubseg;
11686 vertex leftvertex, rightvertex;
11699 apex(*splittri, endpoint1);
11700 org(*splittri, torg);
11701 dest(*splittri, tdest);
11703 tx = tdest[0] - torg[0];
11704 ty = tdest[1] - torg[1];
11705 ex = endpoint2[0] - endpoint1[0];
11706 ey = endpoint2[1] - endpoint1[1];
11707 etx = torg[0] - endpoint2[0];
11708 ety = torg[1] - endpoint2[1];
11709 denom = ty *
ex - tx *
ey;
11710 if (denom == 0.0) {
11711 printf(
"Internal error in segmentintersection():");
11712 printf(
" Attempt to find intersection of parallel segments.\n");
11715 split = (
ey * etx -
ex * ety) / denom;
11719 for (i = 0; i < 2 +
m->nextras; i++) {
11720 newvertex[i] = torg[i] + split * (tdest[i] - torg[i]);
11724 if (
b->verbose > 1) {
11726 " Splitting subsegment (%.12g, %.12g) (%.12g, %.12g) at (%.12g, %.12g).\n",
11727 torg[0], torg[1], tdest[0], tdest[1], newvertex[0], newvertex[1]);
11730 success =
insertvertex(
m,
b, newvertex, splittri, splitsubseg, 0, 0);
11732 printf(
"Internal error in segmentintersection():\n");
11733 printf(
" Failure to split a segment.\n");
11738 if (
m->steinerleft > 0) {
11744 spivot(*splitsubseg, opposubseg);
11750 }
while (splitsubseg->
ss !=
m->dummysub);
11754 }
while (opposubseg.
ss !=
m->dummysub);
11760 dest(*splittri, rightvertex);
11761 apex(*splittri, leftvertex);
11762 if ((leftvertex[0] == endpoint1[0]) && (leftvertex[1] == endpoint1[1])) {
11764 }
else if ((rightvertex[0] != endpoint1[0]) ||
11765 (rightvertex[1] != endpoint1[1])) {
11766 printf(
"Internal error in segmentintersection():\n");
11767 printf(
" Topological inconsistency after splitting a segment.\n");
11797#ifdef ANSI_DECLARATORS
11799 vertex endpoint2,
int newmark)
11804struct otri *searchtri;
11810 struct otri crosstri;
11811 struct osub crosssubseg;
11812 vertex leftvertex, rightvertex;
11817 dest(*searchtri, rightvertex);
11818 apex(*searchtri, leftvertex);
11819 if (((leftvertex[0] == endpoint2[0]) && (leftvertex[1] == endpoint2[1])) ||
11820 ((rightvertex[0] == endpoint2[0]) && (rightvertex[1] == endpoint2[1]))) {
11822 if ((leftvertex[0] == endpoint2[0]) && (leftvertex[1] == endpoint2[1])) {
11843 lnext(*searchtri, crosstri);
11844 tspivot(crosstri, crosssubseg);
11846 if (crosssubseg.
ss ==
m->dummysub) {
11881#ifdef ANSI_DECLARATORS
11885void conformingedge(
m,
b, endpoint1, endpoint2, newmark)
11894 struct otri searchtri1, searchtri2;
11895 struct osub brokensubseg;
11897 vertex midvertex1, midvertex2;
11902 if (
b->verbose > 2) {
11903 printf(
"Forcing segment into triangulation by recursive splitting:\n");
11904 printf(
" (%.12g, %.12g) (%.12g, %.12g)\n", endpoint1[0], endpoint1[1],
11905 endpoint2[0], endpoint2[1]);
11910 for (i = 0; i < 2 +
m->nextras; i++) {
11911 newvertex[i] = 0.5 * (endpoint1[i] + endpoint2[i]);
11916 searchtri1.tri =
m->dummytri;
11921 if (
b->verbose > 2) {
11922 printf(
" Segment intersects existing vertex (%.12g, %.12g).\n",
11923 newvertex[0], newvertex[1]);
11927 org(searchtri1, newvertex);
11930 if (
b->verbose > 2) {
11931 printf(
" Two segments intersect at (%.12g, %.12g).\n",
11932 newvertex[0], newvertex[1]);
11935 tspivot(searchtri1, brokensubseg);
11936 success =
insertvertex(
m,
b, newvertex, &searchtri1, &brokensubseg,
11939 printf(
"Internal error in conformingedge():\n");
11940 printf(
" Failure to split a segment.\n");
11945 if (
m->steinerleft > 0) {
11959 org(searchtri1, midvertex1);
11960 conformingedge(
m,
b, midvertex1, endpoint1, newmark);
11965 org(searchtri2, midvertex2);
11966 conformingedge(
m,
b, midvertex2, endpoint2, newmark);
12011#ifdef ANSI_DECLARATORS
12013 struct otri *fixuptri,
int leftside)
12018struct otri *fixuptri;
12023 struct otri neartri;
12024 struct otri fartri;
12025 struct osub faredge;
12026 vertex nearvertex, leftvertex, rightvertex, farvertex;
12030 lnext(*fixuptri, neartri);
12031 sym(neartri, fartri);
12033 if (fartri.
tri ==
m->dummytri) {
12037 if (faredge.
ss !=
m->dummysub) {
12041 apex(neartri, nearvertex);
12042 org(neartri, leftvertex);
12043 dest(neartri, rightvertex);
12044 apex(fartri, farvertex);
12064 if (
incircle(
m,
b, leftvertex, farvertex, rightvertex, nearvertex) <=
12131#ifdef ANSI_DECLARATORS
12133 struct otri *starttri,
vertex endpoint2,
int newmark)
12138struct otri *starttri;
12144 struct otri fixuptri, fixuptri2;
12145 struct osub crosssubseg;
12154 org(*starttri, endpoint1);
12155 lnext(*starttri, fixuptri);
12162 org(fixuptri, farvertex);
12165 if ((farvertex[0] == endpoint2[0]) && (farvertex[1] == endpoint2[1])) {
12166 oprev(fixuptri, fixuptri2);
12179 oprev(fixuptri, fixuptri2);
12186 oprev(fixuptri, fixuptri2);
12202 tspivot(fixuptri, crosssubseg);
12203 if (crosssubseg.
ss ==
m->dummysub) {
12233#ifdef ANSI_DECLARATORS
12246 struct otri searchtri1, searchtri2;
12251 if (
b->verbose > 1) {
12252 printf(
" Connecting (%.12g, %.12g) to (%.12g, %.12g).\n",
12253 endpoint1[0], endpoint1[1], endpoint2[0], endpoint2[1]);
12257 checkvertex = (
vertex) NULL;
12259 if (encodedtri != (
triangle) NULL) {
12260 decode(encodedtri, searchtri1);
12261 org(searchtri1, checkvertex);
12263 if (checkvertex != endpoint1) {
12265 searchtri1.
tri =
m->dummytri;
12271 "Internal error in insertsegment(): Unable to locate PSLG vertex\n");
12272 printf(
" (%.12g, %.12g) in triangulation.\n",
12273 endpoint1[0], endpoint1[1]);
12287 org(searchtri1, endpoint1);
12290 checkvertex = (
vertex) NULL;
12292 if (encodedtri != (
triangle) NULL) {
12293 decode(encodedtri, searchtri2);
12294 org(searchtri2, checkvertex);
12296 if (checkvertex != endpoint2) {
12298 searchtri2.
tri =
m->dummytri;
12304 "Internal error in insertsegment(): Unable to locate PSLG vertex\n");
12305 printf(
" (%.12g, %.12g) in triangulation.\n",
12306 endpoint2[0], endpoint2[1]);
12320 org(searchtri2, endpoint2);
12326 conformingedge(
m,
b, endpoint1, endpoint2, newmark);
12345#ifdef ANSI_DECLARATORS
12354 struct otri hulltri;
12355 struct otri nexttri;
12356 struct otri starttri;
12360 hulltri.
tri =
m->dummytri;
12371 oprev(hulltri, nexttri);
12372 while (nexttri.
tri !=
m->dummytri) {
12374 oprev(hulltri, nexttri);
12376 }
while (!
otriequal(hulltri, starttri));
12391#ifdef ANSI_DECLARATORS
12393 int *segmentmarkerlist,
int numberofsegments)
12395void formskeleton(
m,
b, segmentlist, segmentmarkerlist, numberofsegments)
12399int *segmentmarkerlist;
12400int numberofsegments;
12405#ifdef ANSI_DECLARATORS
12407 FILE *polyfile,
char *polyfilename)
12420 char polyfilename[6];
12426 vertex endpoint1, endpoint2;
12427 int segmentmarkers;
12434 printf(
"Recovering segments in Delaunay triangulation.\n");
12437 strcpy(polyfilename,
"input");
12438 m->insegments = numberofsegments;
12439 segmentmarkers = segmentmarkerlist != (
int *) NULL;
12444 stringptr =
readline(inputline, polyfile, polyfilename);
12445 m->insegments = (
int) strtol(stringptr, &stringptr, 0);
12446 stringptr = findfield(stringptr);
12447 if (*stringptr ==
'\0') {
12448 segmentmarkers = 0;
12450 segmentmarkers = (
int) strtol(stringptr, &stringptr, 0);
12455 if (
m->triangles.items == 0) {
12461 if (
m->insegments > 0) {
12464 printf(
" Recovering PSLG segments.\n");
12470 for (i = 0; i <
m->insegments; i++) {
12472 end1 = segmentlist[
index++];
12473 end2 = segmentlist[
index++];
12474 if (segmentmarkers) {
12475 boundmarker = segmentmarkerlist[i];
12478 stringptr =
readline(inputline, polyfile,
b->inpolyfilename);
12479 stringptr = findfield(stringptr);
12480 if (*stringptr ==
'\0') {
12481 printf(
"Error: Segment %d has no endpoints in %s.\n",
12482 b->firstnumber + i, polyfilename);
12485 end1 = (
int) strtol(stringptr, &stringptr, 0);
12487 stringptr = findfield(stringptr);
12488 if (*stringptr ==
'\0') {
12489 printf(
"Error: Segment %d is missing its second endpoint in %s.\n",
12490 b->firstnumber + i, polyfilename);
12493 end2 = (
int) strtol(stringptr, &stringptr, 0);
12495 if (segmentmarkers) {
12496 stringptr = findfield(stringptr);
12497 if (*stringptr ==
'\0') {
12500 boundmarker = (
int) strtol(stringptr, &stringptr, 0);
12505 (end1 >=
b->firstnumber +
m->invertices)) {
12507 printf(
"Warning: Invalid first endpoint of segment %d in %s.\n",
12508 b->firstnumber + i, polyfilename);
12511 (end2 >=
b->firstnumber +
m->invertices)) {
12513 printf(
"Warning: Invalid second endpoint of segment %d in %s.\n",
12514 b->firstnumber + i, polyfilename);
12520 if ((endpoint1[0] == endpoint2[0]) && (endpoint1[1] == endpoint2[1])) {
12522 printf(
"Warning: Endpoints of segment %d are coincident in %s.\n",
12523 b->firstnumber + i, polyfilename);
12533 if (
b->convex || !
b->poly) {
12536 printf(
" Enclosing convex hull with segments.\n");
12558#ifdef ANSI_DECLARATORS
12567 struct otri hulltri;
12568 struct otri nexttri;
12569 struct otri starttri;
12570 struct osub hullsubseg;
12577 printf(
" Marking concavities (external triangles) for elimination.\n");
12580 hulltri.
tri =
m->dummytri;
12590 tspivot(hulltri, hullsubseg);
12591 if (hullsubseg.
ss ==
m->dummysub) {
12596 *deadtriangle = hulltri.
tri;
12600 if (
mark(hullsubseg) == 0) {
12602 org(hulltri, horg);
12603 dest(hulltri, hdest);
12615 oprev(hulltri, nexttri);
12616 while (nexttri.
tri !=
m->dummytri) {
12618 oprev(hulltri, nexttri);
12620 }
while (!
otriequal(hulltri, starttri));
12640#ifdef ANSI_DECLARATORS
12649 struct otri testtri;
12650 struct otri neighbor;
12653 struct osub neighborsubseg;
12656 vertex deadorg, deaddest, deadapex;
12662 printf(
" Marking neighbors of marked triangles.\n");
12668 while (virusloop != (
triangle **) NULL) {
12669 testtri.
tri = *virusloop;
12675 if (
b->verbose > 2) {
12679 org(testtri, deadorg);
12680 dest(testtri, deaddest);
12681 apex(testtri, deadapex);
12682 printf(
" Checking (%.12g, %.12g) (%.12g, %.12g) (%.12g, %.12g)\n",
12683 deadorg[0], deadorg[1], deaddest[0], deaddest[1],
12684 deadapex[0], deadapex[1]);
12689 sym(testtri, neighbor);
12691 tspivot(testtri, neighborsubseg);
12693 if ((neighbor.
tri ==
m->dummytri) ||
infected(neighbor)) {
12694 if (neighborsubseg.
ss !=
m->dummysub) {
12699 if (neighbor.
tri !=
m->dummytri) {
12708 if (neighborsubseg.
ss ==
m->dummysub) {
12711 if (
b->verbose > 2) {
12712 org(neighbor, deadorg);
12713 dest(neighbor, deaddest);
12714 apex(neighbor, deadapex);
12716 " Marking (%.12g, %.12g) (%.12g, %.12g) (%.12g, %.12g)\n",
12717 deadorg[0], deadorg[1], deaddest[0], deaddest[1],
12718 deadapex[0], deadapex[1]);
12723 *deadtriangle = neighbor.
tri;
12728 if (
mark(neighborsubseg) == 0) {
12731 org(neighbor, norg);
12732 dest(neighbor, ndest);
12749 printf(
" Deleting marked triangles.\n");
12754 while (virusloop != (
triangle **) NULL) {
12755 testtri.
tri = *virusloop;
12761 org(testtri, testvertex);
12763 if (testvertex != (
vertex) NULL) {
12768 onext(testtri, neighbor);
12770 while ((neighbor.
tri !=
m->dummytri) &&
12783 if (neighbor.
tri ==
m->dummytri) {
12785 oprev(testtri, neighbor);
12787 while (neighbor.
tri !=
m->dummytri) {
12800 if (
b->verbose > 1) {
12801 printf(
" Deleting vertex (%.12g, %.12g)\n",
12802 testvertex[0], testvertex[1]);
12813 sym(testtri, neighbor);
12814 if (neighbor.
tri ==
m->dummytri) {
12850#ifdef ANSI_DECLARATORS
12862 struct otri testtri;
12863 struct otri neighbor;
12866 struct osub neighborsubseg;
12867 vertex regionorg, regiondest, regionapex;
12871 if (
b->verbose > 1) {
12872 printf(
" Marking neighbors of marked triangles.\n");
12879 while (virusloop != (
triangle **) NULL) {
12880 testtri.
tri = *virusloop;
12886 if (
b->regionattrib) {
12894 if (
b->verbose > 2) {
12898 org(testtri, regionorg);
12899 dest(testtri, regiondest);
12900 apex(testtri, regionapex);
12901 printf(
" Checking (%.12g, %.12g) (%.12g, %.12g) (%.12g, %.12g)\n",
12902 regionorg[0], regionorg[1], regiondest[0], regiondest[1],
12903 regionapex[0], regionapex[1]);
12908 sym(testtri, neighbor);
12910 tspivot(testtri, neighborsubseg);
12913 if ((neighbor.
tri !=
m->dummytri) && !
infected(neighbor)
12914 && (neighborsubseg.
ss ==
m->dummysub)) {
12915 if (
b->verbose > 2) {
12916 org(neighbor, regionorg);
12917 dest(neighbor, regiondest);
12918 apex(neighbor, regionapex);
12919 printf(
" Marking (%.12g, %.12g) (%.12g, %.12g) (%.12g, %.12g)\n",
12920 regionorg[0], regionorg[1], regiondest[0], regiondest[1],
12921 regionapex[0], regionapex[1]);
12927 *regiontri = neighbor.
tri;
12937 if (
b->verbose > 1) {
12938 printf(
" Unmarking marked triangles.\n");
12942 while (virusloop != (
triangle **) NULL) {
12943 testtri.
tri = *virusloop;
12963#ifdef ANSI_DECLARATORS
12965 REAL *regionlist,
int regions)
12967void carveholes(
m,
b, holelist, holes, regionlist, regions)
12977 struct otri searchtri;
12978 struct otri triangleloop;
12979 struct otri *regiontris;
12982 vertex searchorg, searchdest;
12987 if (!(
b->quiet || (
b->noholes &&
b->convex))) {
12988 printf(
"Removing unwanted triangles.\n");
12989 if (
b->verbose && (holes > 0)) {
12990 printf(
" Marking holes for elimination.\n");
12997 (
int)
sizeof(
struct otri));
12999 regiontris = (
struct otri *) NULL;
13002 if (((holes > 0) && !
b->noholes) || !
b->convex || (regions > 0)) {
13014 if ((holes > 0) && !
b->noholes) {
13016 for (i = 0; i < 2 * holes; i += 2) {
13018 if ((holelist[i] >=
m->xmin) && (holelist[i] <=
m->xmax)
13019 && (holelist[i + 1] >=
m->ymin) && (holelist[i + 1] <=
m->ymax)) {
13021 searchtri.
tri =
m->dummytri;
13027 org(searchtri, searchorg);
13028 dest(searchtri, searchdest);
13032 intersect =
locate(
m,
b, &holelist[i], &searchtri);
13038 *holetri = searchtri.
tri;
13053 for (i = 0; i < regions; i++) {
13054 regiontris[i].
tri =
m->dummytri;
13056 if ((regionlist[4 * i] >=
m->xmin) && (regionlist[4 * i] <=
m->xmax) &&
13057 (regionlist[4 * i + 1] >=
m->ymin) &&
13058 (regionlist[4 * i + 1] <=
m->ymax)) {
13060 searchtri.
tri =
m->dummytri;
13066 org(searchtri, searchorg);
13067 dest(searchtri, searchdest);
13071 intersect =
locate(
m,
b, ®ionlist[4 * i], &searchtri);
13075 otricopy(searchtri, regiontris[i]);
13082 if (
m->viri.items > 0) {
13090 if (
b->regionattrib) {
13092 printf(
"Spreading regional attributes and area constraints.\n");
13094 printf(
"Spreading regional attributes.\n");
13097 printf(
"Spreading regional area constraints.\n");
13100 if (
b->regionattrib && !
b->refine) {
13103 triangleloop.
orient = 0;
13110 for (i = 0; i < regions; i++) {
13111 if (regiontris[i].
tri !=
m->dummytri) {
13118 *regiontri = regiontris[i].
tri;
13120 regionplague(
m,
b, regionlist[4 * i + 2], regionlist[4 * i + 3]);
13125 if (
b->regionattrib && !
b->refine) {
13132 if (((holes > 0) && !
b->noholes) || !
b->convex || (regions > 0)) {
13157#ifdef ANSI_DECLARATORS
13160void tallyencs(
m,
b)
13166 struct osub subsegloop;
13170 subsegloop.ssorient = 0;
13172 while (subsegloop.ss != (
subseg *) NULL) {
13174 dummy = checkseg4encroach(
m,
b, &subsegloop);
13189void precisionerror()
13191 printf(
"Try increasing the area criterion and/or reducing the minimum\n");
13192 printf(
" allowable angle so that tiny triangles are not created.\n");
13194 printf(
"Alternatively, try recompiling me with double precision\n");
13195 printf(
" arithmetic (by removing \"#define SINGLE\" from the\n");
13196 printf(
" source file or \"-DSINGLE\" from the makefile).\n");
13218#ifdef ANSI_DECLARATORS
13219void splitencsegs(
struct mesh *
m,
struct behavior *
b,
int triflaws)
13221void splitencsegs(
m,
b, triflaws)
13228 struct otri enctri;
13229 struct otri testtri;
13230 struct osub testsh;
13231 struct osub currentenc;
13233 vertex eorg, edest, eapex;
13236 REAL segmentlength, nearestpoweroftwo;
13238 REAL multiplier, divisor;
13239 int acuteorg, acuteorg2, acutedest, acutedest2;
13247 while ((
m->badsubsegs.items > 0) && (
m->steinerleft != 0)) {
13249 encloop = badsubsegtraverse(
m);
13250 while ((encloop != (
struct badsubseg *) NULL) && (
m->steinerleft != 0)) {
13252 sorg(currentenc, eorg);
13253 sdest(currentenc, edest);
13277 lnext(enctri, testtri);
13279 acuteorg = testsh.ss !=
m->dummysub;
13283 acutedest = testsh.ss !=
m->dummysub;
13288 if (!
b->conformdel && !acuteorg && !acutedest) {
13289 apex(enctri, eapex);
13291 ((eorg[0] - eapex[0]) * (edest[0] - eapex[0]) +
13292 (eorg[1] - eapex[1]) * (edest[1] - eapex[1]) < 0.0)) {
13293 deletevertex(
m,
b, &testtri);
13295 apex(enctri, eapex);
13296 lprev(enctri, testtri);
13302 sym(enctri, testtri);
13303 if (testtri.
tri !=
m->dummytri) {
13307 acutedest2 = testsh.ss !=
m->dummysub;
13308 acutedest = acutedest || acutedest2;
13312 acuteorg2 = testsh.ss !=
m->dummysub;
13313 acuteorg = acuteorg || acuteorg2;
13316 if (!
b->conformdel && !acuteorg2 && !acutedest2) {
13317 org(testtri, eapex);
13319 ((eorg[0] - eapex[0]) * (edest[0] - eapex[0]) +
13320 (eorg[1] - eapex[1]) * (edest[1] - eapex[1]) < 0.0)) {
13321 deletevertex(
m,
b, &testtri);
13322 sym(enctri, testtri);
13323 apex(testtri, eapex);
13331 if (acuteorg || acutedest) {
13332 segmentlength =
sqrt((edest[0] - eorg[0]) * (edest[0] - eorg[0]) +
13333 (edest[1] - eorg[1]) * (edest[1] - eorg[1]));
13336 nearestpoweroftwo = 1.0;
13337 while (segmentlength > 3.0 * nearestpoweroftwo) {
13338 nearestpoweroftwo *= 2.0;
13340 while (segmentlength < 1.5 * nearestpoweroftwo) {
13341 nearestpoweroftwo *= 0.5;
13344 split = nearestpoweroftwo / segmentlength;
13346 split = 1.0 - split;
13357 for (i = 0; i < 2 +
m->nextras; i++) {
13358 newvertex[i] = eorg[i] + split * (edest[i] - eorg[i]);
13366 divisor = ((eorg[0] - edest[0]) * (eorg[0] - edest[0]) +
13367 (eorg[1] - edest[1]) * (eorg[1] - edest[1]));
13368 if ((multiplier != 0.0) && (divisor != 0.0)) {
13369 multiplier = multiplier / divisor;
13371 if (multiplier == multiplier) {
13372 newvertex[0] += multiplier * (edest[1] - eorg[1]);
13373 newvertex[1] += multiplier * (eorg[0] - edest[0]);
13380 if (
b->verbose > 1) {
13382 " Splitting subsegment (%.12g, %.12g) (%.12g, %.12g) at (%.12g, %.12g).\n",
13383 eorg[0], eorg[1], edest[0], edest[1],
13384 newvertex[0], newvertex[1]);
13387 if (((newvertex[0] == eorg[0]) && (newvertex[1] == eorg[1])) ||
13388 ((newvertex[0] == edest[0]) && (newvertex[1] == edest[1]))) {
13389 printf(
"Error: Ran out of precision at (%.12g, %.12g).\n",
13390 newvertex[0], newvertex[1]);
13391 printf(
"I attempted to split a segment to a smaller size than\n");
13392 printf(
" can be accommodated by the finite precision of\n");
13393 printf(
" floating point arithmetic.\n");
13401 printf(
"Internal error in splitencsegs():\n");
13402 printf(
" Failure to split a segment.\n");
13405 if (
m->steinerleft > 0) {
13409 dummy = checkseg4encroach(
m,
b, ¤tenc);
13411 dummy = checkseg4encroach(
m,
b, ¤tenc);
13414 badsubsegdealloc(
m, encloop);
13415 encloop = badsubsegtraverse(
m);
13430#ifdef ANSI_DECLARATORS
13433void tallyfaces(
m,
b)
13439 struct otri triangleloop;
13442 printf(
" Making a list of bad triangles.\n");
13445 triangleloop.orient = 0;
13447 while (triangleloop.tri != (
triangle *) NULL) {
13449 testtriangle(
m,
b, &triangleloop);
13466#ifdef ANSI_DECLARATORS
13470void splittriangle(
m,
b, badtri)
13477 struct otri badotri;
13478 vertex borg, bdest, bapex;
13486 org(badotri, borg);
13487 dest(badotri, bdest);
13488 apex(badotri, bapex);
13494 if (
b->verbose > 1) {
13495 printf(
" Splitting this triangle at its circumcenter:\n");
13496 printf(
" (%.12g, %.12g) (%.12g, %.12g) (%.12g, %.12g)\n", borg[0],
13497 borg[1], bdest[0], bdest[1], bapex[0], bapex[1]);
13506 if (((newvertex[0] == borg[0]) && (newvertex[1] == borg[1])) ||
13507 ((newvertex[0] == bdest[0]) && (newvertex[1] == bdest[1])) ||
13508 ((newvertex[0] == bapex[0]) && (newvertex[1] == bapex[1]))) {
13511 "Warning: New vertex (%.12g, %.12g) falls on existing vertex.\n",
13512 newvertex[0], newvertex[1]);
13517 for (i = 2; i < 2 +
m->nextras; i++) {
13519 newvertex[i] = borg[i] + xi * (bdest[i] - borg[i])
13520 + eta * (bapex[i] - borg[i]);
13543 if (
m->steinerleft > 0) {
13550 if (
b->verbose > 1) {
13551 printf(
" Rejecting (%.12g, %.12g).\n", newvertex[0], newvertex[1]);
13562 "Warning: New vertex (%.12g, %.12g) falls on existing vertex.\n",
13563 newvertex[0], newvertex[1]);
13571 printf(
" The new vertex is at the circumcenter of triangle\n");
13572 printf(
" (%.12g, %.12g) (%.12g, %.12g) (%.12g, %.12g)\n",
13573 borg[0], borg[1], bdest[0], bdest[1], bapex[0], bapex[1]);
13575 printf(
"This probably means that I am trying to refine triangles\n");
13576 printf(
" to a smaller size than can be accommodated by the finite\n");
13577 printf(
" precision of floating point arithmetic. (You can be\n");
13578 printf(
" sure of this if I fail to terminate.)\n");
13595#ifdef ANSI_DECLARATORS
13598void enforcequality(
m,
b)
13608 printf(
"Adding Steiner points to enforce quality.\n");
13614 printf(
" Looking for encroached subsegments.\n");
13618 if (
b->verbose && (
m->badsubsegs.items > 0)) {
13619 printf(
" Splitting encroached subsegments.\n");
13622 splitencsegs(
m,
b, 0);
13627 if ((
b->minangle > 0.0) ||
b->vararea ||
b->fixedarea ||
b->usertest) {
13632 for (i = 0; i < 4096; i++) {
13633 m->queuefront[i] = (
struct badtriang *) NULL;
13635 m->firstnonemptyq = -1;
13641 m->checkquality = 1;
13643 printf(
" Splitting bad triangles.\n");
13645 while ((
m->badtriangles.items > 0) && (
m->steinerleft != 0)) {
13647 badtri = dequeuebadtriang(
m);
13648 splittriangle(
m,
b, badtri);
13649 if (
m->badsubsegs.items > 0) {
13651 enqueuebadtriang(
m,
b, badtri);
13654 splitencsegs(
m,
b, 1);
13666 if (!
b->quiet &&
b->conformdel && (
m->badsubsegs.items > 0) &&
13667 (
m->steinerleft == 0)) {
13668 printf(
"\nWarning: I ran out of Steiner points, but the mesh has\n");
13669 if (
m->badsubsegs.items == 1) {
13670 printf(
" one encroached subsegment, and therefore might not be truly\n"
13673 printf(
" %ld encroached subsegments, and therefore might not be truly\n"
13674 ,
m->badsubsegs.items);
13676 printf(
" Delaunay. If the Delaunay property is important to you,\n");
13677 printf(
" try increasing the number of Steiner points (controlled by\n");
13678 printf(
" the -S switch) slightly and try again.\n\n");
13694#ifdef ANSI_DECLARATORS
13703 struct otri triangleloop, trisym;
13704 struct osub checkmark;
13712 printf(
"Adding vertices for second-order triangles.\n");
13719 m->vertices.deaditemstack = (
VOID *) NULL;
13730 for (triangleloop.
orient = 0; triangleloop.
orient < 3;
13731 triangleloop.
orient++) {
13732 sym(triangleloop, trisym);
13733 if ((triangleloop.
tri < trisym.
tri) || (trisym.
tri ==
m->dummytri)) {
13734 org(triangleloop, torg);
13735 dest(triangleloop, tdest);
13739 for (i = 0; i < 2 +
m->nextras; i++) {
13740 newvertex[i] = 0.5 * (torg[i] + tdest[i]);
13747 if (
b->usesegments) {
13748 tspivot(triangleloop, checkmark);
13750 if (checkmark.
ss !=
m->dummysub) {
13755 if (
b->verbose > 1) {
13756 printf(
" Creating (%.12g, %.12g).\n", newvertex[0], newvertex[1]);
13759 triangleloop.
tri[
m->highorderindex + triangleloop.
orient] =
13761 if (trisym.
tri !=
m->dummytri) {
13785#ifdef ANSI_DECLARATORS
13786char *
readline(
char *
string, FILE *infile,
char *infilename)
13788char *
readline(
string, infile, infilename)
13800 if (
result == (
char *) NULL) {
13801 printf(
" Error: Unexpected end of file in %s.\n", infilename);
13829#ifdef ANSI_DECLARATORS
13830char *findfield(
char *
string)
13832char *findfield(
string)
13870#ifdef ANSI_DECLARATORS
13871void readnodes(
struct mesh *
m,
struct behavior *
b,
char *nodefilename,
13872 char *polyfilename, FILE **polyfile)
13874void readnodes(
m,
b, nodefilename, polyfilename, polyfile)
13897 printf(
"Opening %s.\n", polyfilename);
13899 *polyfile = fopen(polyfilename,
"r");
13900 if (*polyfile == (FILE *) NULL) {
13901 printf(
" Error: Cannot access file %s.\n", polyfilename);
13906 stringptr =
readline(inputline, *polyfile, polyfilename);
13907 m->invertices = (
int) strtol(stringptr, &stringptr, 0);
13908 stringptr = findfield(stringptr);
13909 if (*stringptr ==
'\0') {
13912 m->mesh_dim = (
int) strtol(stringptr, &stringptr, 0);
13914 stringptr = findfield(stringptr);
13915 if (*stringptr ==
'\0') {
13918 m->nextras = (
int) strtol(stringptr, &stringptr, 0);
13920 stringptr = findfield(stringptr);
13921 if (*stringptr ==
'\0') {
13924 nodemarkers = (
int) strtol(stringptr, &stringptr, 0);
13926 if (
m->invertices > 0) {
13927 infile = *polyfile;
13928 infilename = polyfilename;
13929 m->readnodefile = 0;
13933 m->readnodefile = 1;
13934 infilename = nodefilename;
13937 m->readnodefile = 1;
13938 infilename = nodefilename;
13939 *polyfile = (FILE *) NULL;
13942 if (
m->readnodefile) {
13945 printf(
"Opening %s.\n", nodefilename);
13947 infile = fopen(nodefilename,
"r");
13948 if (infile == (FILE *) NULL) {
13949 printf(
" Error: Cannot access file %s.\n", nodefilename);
13954 stringptr =
readline(inputline, infile, nodefilename);
13955 m->invertices = (
int) strtol(stringptr, &stringptr, 0);
13956 stringptr = findfield(stringptr);
13957 if (*stringptr ==
'\0') {
13960 m->mesh_dim = (
int) strtol(stringptr, &stringptr, 0);
13962 stringptr = findfield(stringptr);
13963 if (*stringptr ==
'\0') {
13966 m->nextras = (
int) strtol(stringptr, &stringptr, 0);
13968 stringptr = findfield(stringptr);
13969 if (*stringptr ==
'\0') {
13972 nodemarkers = (
int) strtol(stringptr, &stringptr, 0);
13976 if (
m->invertices < 3) {
13977 printf(
"Error: Input must have at least three input vertices.\n");
13980 if (
m->mesh_dim != 2) {
13981 printf(
"Error: Triangle only works with two-dimensional meshes.\n");
13984 if (
m->nextras == 0) {
13991 for (i = 0; i <
m->invertices; i++) {
13993 stringptr =
readline(inputline, infile, infilename);
13995 firstnode = (
int) strtol(stringptr, &stringptr, 0);
13996 if ((firstnode == 0) || (firstnode == 1)) {
13997 b->firstnumber = firstnode;
14000 stringptr = findfield(stringptr);
14001 if (*stringptr ==
'\0') {
14002 printf(
"Error: Vertex %d has no x coordinate.\n",
b->firstnumber + i);
14005 x = (
REAL) strtod(stringptr, &stringptr);
14006 stringptr = findfield(stringptr);
14007 if (*stringptr ==
'\0') {
14008 printf(
"Error: Vertex %d has no y coordinate.\n",
b->firstnumber + i);
14011 y = (
REAL) strtod(stringptr, &stringptr);
14015 for (j = 2; j < 2 +
m->nextras; j++) {
14016 stringptr = findfield(stringptr);
14017 if (*stringptr ==
'\0') {
14018 vertexloop[j] = 0.0;
14020 vertexloop[j] = (
REAL) strtod(stringptr, &stringptr);
14025 stringptr = findfield(stringptr);
14026 if (*stringptr ==
'\0') {
14029 currentmarker = (
int) strtol(stringptr, &stringptr, 0);
14039 m->xmin =
m->xmax =
x;
14040 m->ymin =
m->ymax =
y;
14048 if (
m->readnodefile) {
14054 m->xminextreme = 10 *
m->xmin - 9 *
m->xmax;
14067#ifdef ANSI_DECLARATORS
14069 REAL *pointattriblist,
int *pointmarkerlist,
14070 int numberofpoints,
int numberofpointattribs)
14073 numberofpoints, numberofpointattribs)
14077REAL *pointattriblist;
14078int *pointmarkerlist;
14080int numberofpointattribs;
14090 m->invertices = numberofpoints;
14092 m->nextras = numberofpointattribs;
14093 m->readnodefile = 0;
14094 if (
m->invertices < 3) {
14095 printf(
"Error: Input must have at least three input vertices.\n");
14098 if (
m->nextras == 0) {
14107 for (i = 0; i <
m->invertices; i++) {
14110 x = vertexloop[0] = pointlist[coordindex++];
14111 y = vertexloop[1] = pointlist[coordindex++];
14113 for (j = 0; j < numberofpointattribs; j++) {
14114 vertexloop[2 + j] = pointattriblist[attribindex++];
14116 if (pointmarkerlist != (
int *) NULL) {
14126 m->xmin =
m->xmax =
x;
14127 m->ymin =
m->ymax =
y;
14129 m->xmin = (
x <
m->xmin) ?
x :
m->xmin;
14130 m->xmax = (
x >
m->xmax) ?
x :
m->xmax;
14131 m->ymin = (
y <
m->ymin) ?
y :
m->ymin;
14132 m->ymax = (
y >
m->ymax) ?
y :
m->ymax;
14138 m->xminextreme = 10 *
m->xmin - 9 *
m->xmax;
14152#ifdef ANSI_DECLARATORS
14154 FILE *polyfile,
char *polyfilename,
REAL **hlist,
int *holes,
14155 REAL **rlist,
int *regions)
14157void readholes(
m,
b, polyfile, polyfilename, hlist, holes, rlist, regions)
14177 stringptr =
readline(inputline, polyfile, polyfilename);
14178 *holes = (
int) strtol(stringptr, &stringptr, 0);
14182 for (i = 0; i < 2 * *holes; i += 2) {
14183 stringptr =
readline(inputline, polyfile, polyfilename);
14184 stringptr = findfield(stringptr);
14185 if (*stringptr ==
'\0') {
14186 printf(
"Error: Hole %d has no x coordinate.\n",
14187 b->firstnumber + (i >> 1));
14190 holelist[i] = (
REAL) strtod(stringptr, &stringptr);
14192 stringptr = findfield(stringptr);
14193 if (*stringptr ==
'\0') {
14194 printf(
"Error: Hole %d has no y coordinate.\n",
14195 b->firstnumber + (i >> 1));
14198 holelist[i + 1] = (
REAL) strtod(stringptr, &stringptr);
14202 *hlist = (
REAL *) NULL;
14206 if ((
b->regionattrib ||
b->vararea) && !
b->refine) {
14208 stringptr =
readline(inputline, polyfile, polyfilename);
14209 *regions = (
int) strtol(stringptr, &stringptr, 0);
14210 if (*regions > 0) {
14212 *rlist = regionlist;
14214 for (i = 0; i < *regions; i++) {
14215 stringptr =
readline(inputline, polyfile, polyfilename);
14216 stringptr = findfield(stringptr);
14217 if (*stringptr ==
'\0') {
14218 printf(
"Error: Region %d has no x coordinate.\n",
14219 b->firstnumber + i);
14222 regionlist[
index++] = (
REAL) strtod(stringptr, &stringptr);
14224 stringptr = findfield(stringptr);
14225 if (*stringptr ==
'\0') {
14226 printf(
"Error: Region %d has no y coordinate.\n",
14227 b->firstnumber + i);
14230 regionlist[
index++] = (
REAL) strtod(stringptr, &stringptr);
14232 stringptr = findfield(stringptr);
14233 if (*stringptr ==
'\0') {
14235 "Error: Region %d has no region attribute or area constraint.\n",
14236 b->firstnumber + i);
14239 regionlist[
index++] = (
REAL) strtod(stringptr, &stringptr);
14241 stringptr = findfield(stringptr);
14242 if (*stringptr ==
'\0') {
14245 regionlist[
index] = (
REAL) strtod(stringptr, &stringptr);
14253 *rlist = (
REAL *) NULL;
14271#ifdef ANSI_DECLARATORS
14272void finishfile(FILE *outfile,
int argc,
char **argv)
14274void finishfile(outfile, argc, argv)
14283 fprintf(outfile,
"# Generated by");
14284 for (i = 0; i < argc; i++) {
14285 fprintf(outfile,
" ");
14286 fputs(argv[i], outfile);
14288 fprintf(outfile,
"\n");
14305#ifdef ANSI_DECLARATORS
14307 REAL **pointattriblist,
int **pointmarkerlist)
14309void writenodes(
m,
b, pointlist, pointattriblist, pointmarkerlist)
14313REAL **pointattriblist;
14314int **pointmarkerlist;
14319#ifdef ANSI_DECLARATORS
14321 int argc,
char **argv)
14349 outvertices =
m->vertices.items -
m->undeads;
14351 outvertices =
m->vertices.items;
14356 printf(
"Writing vertices.\n");
14359 if (*pointlist == (
REAL *) NULL) {
14363 if ((
m->nextras > 0) && (*pointattriblist == (
REAL *) NULL)) {
14364 *pointattriblist = (
REAL *)
trimalloc((
int) (outvertices *
m->nextras *
14368 if (!
b->nobound && (*pointmarkerlist == (
int *) NULL)) {
14369 *pointmarkerlist = (
int *)
trimalloc((
int) (outvertices *
sizeof(
int)));
14371 plist = *pointlist;
14372 palist = *pointattriblist;
14373 pmlist = *pointmarkerlist;
14378 printf(
"Writing %s.\n", nodefilename);
14380 outfile = fopen(nodefilename,
"w");
14381 if (outfile == (FILE *) NULL) {
14382 printf(
" Error: Cannot create file %s.\n", nodefilename);
14387 fprintf(outfile,
"%ld %d %d %d\n", outvertices,
m->mesh_dim,
14388 m->nextras, 1 -
b->nobound);
14392 vertexnumber =
b->firstnumber;
14394 while (vertexloop != (
vertex) NULL) {
14398 plist[coordindex++] = vertexloop[0];
14399 plist[coordindex++] = vertexloop[1];
14401 for (i = 0; i <
m->nextras; i++) {
14402 palist[attribindex++] = vertexloop[2 + i];
14406 pmlist[vertexnumber -
b->firstnumber] =
vertexmark(vertexloop);
14410 fprintf(outfile,
"%4d %.17g %.17g", vertexnumber, vertexloop[0],
14412 for (i = 0; i <
m->nextras; i++) {
14414 fprintf(outfile,
" %.17g", vertexloop[i + 2]);
14417 fprintf(outfile,
"\n");
14420 fprintf(outfile,
" %d\n",
vertexmark(vertexloop));
14431 finishfile(outfile, argc, argv);
14445#ifdef ANSI_DECLARATORS
14458 vertexnumber =
b->firstnumber;
14460 while (vertexloop != (
vertex) NULL) {
14477#ifdef ANSI_DECLARATORS
14479 int **trianglelist,
REAL **triangleattriblist)
14485REAL **triangleattriblist;
14490#ifdef ANSI_DECLARATORS
14492 int argc,
char **argv)
14513 struct otri triangleloop;
14515 vertex mid1, mid2, mid3;
14517 long elementnumber;
14523 printf(
"Writing triangles.\n");
14526 if (*trianglelist == (
int *) NULL) {
14527 *trianglelist = (
int *)
trimalloc((
int) (
m->triangles.items *
14528 ((
b->order + 1) * (
b->order + 2) /
14529 2) *
sizeof(
int)));
14532 if ((
m->eextras > 0) && (*triangleattriblist == (
REAL *) NULL)) {
14533 *triangleattriblist = (
REAL *)
trimalloc((
int) (
m->triangles.items *
14537 tlist = *trianglelist;
14538 talist = *triangleattriblist;
14543 printf(
"Writing %s.\n", elefilename);
14545 outfile = fopen(elefilename,
"w");
14546 if (outfile == (FILE *) NULL) {
14547 printf(
" Error: Cannot create file %s.\n", elefilename);
14551 fprintf(outfile,
"%ld %d %d\n",
m->triangles.items,
14552 (
b->order + 1) * (
b->order + 2) / 2,
m->eextras);
14557 triangleloop.
orient = 0;
14559 elementnumber =
b->firstnumber;
14562 org(triangleloop, p1);
14563 dest(triangleloop, p2);
14564 apex(triangleloop, p3);
14565 if (
b->order == 1) {
14572 fprintf(outfile,
"%4ld %4d %4d %4d", elementnumber,
14576 mid1 = (
vertex) triangleloop.
tri[
m->highorderindex + 1];
14577 mid2 = (
vertex) triangleloop.
tri[
m->highorderindex + 2];
14578 mid3 = (
vertex) triangleloop.
tri[
m->highorderindex];
14588 fprintf(outfile,
"%4ld %4d %4d %4d %4d %4d %4d", elementnumber,
14595 for (i = 0; i <
m->eextras; i++) {
14599 for (i = 0; i <
m->eextras; i++) {
14600 fprintf(outfile,
" %.17g",
elemattribute(triangleloop, i));
14602 fprintf(outfile,
"\n");
14612 finishfile(outfile, argc, argv);
14624#ifdef ANSI_DECLARATORS
14626 int **segmentlist,
int **segmentmarkerlist)
14628void writepoly(
m,
b, segmentlist, segmentmarkerlist)
14632int **segmentmarkerlist;
14637#ifdef ANSI_DECLARATORS
14639 REAL *holelist,
int holes,
REAL *regionlist,
int regions,
14640 int argc,
char **argv)
14642void writepoly(
m,
b, polyfilename, holelist, holes, regionlist, regions,
14664 long holenumber, regionnumber;
14666 struct osub subsegloop;
14667 vertex endpoint1, endpoint2;
14672 printf(
"Writing segments.\n");
14675 if (*segmentlist == (
int *) NULL) {
14676 *segmentlist = (
int *)
trimalloc((
int) (
m->subsegs.items * 2 *
14680 if (!
b->nobound && (*segmentmarkerlist == (
int *) NULL)) {
14681 *segmentmarkerlist = (
int *)
trimalloc((
int) (
m->subsegs.items *
14684 slist = *segmentlist;
14685 smlist = *segmentmarkerlist;
14689 printf(
"Writing %s.\n", polyfilename);
14691 outfile = fopen(polyfilename,
"w");
14692 if (outfile == (FILE *) NULL) {
14693 printf(
" Error: Cannot create file %s.\n", polyfilename);
14699 fprintf(outfile,
"%d %d %d %d\n", 0,
m->mesh_dim,
m->nextras,
14702 fprintf(outfile,
"%ld %d\n",
m->subsegs.items, 1 -
b->nobound);
14708 subsegnumber =
b->firstnumber;
14709 while (subsegloop.
ss != (
subseg *) NULL) {
14710 sorg(subsegloop, endpoint1);
14711 sdest(subsegloop, endpoint2);
14718 smlist[subsegnumber -
b->firstnumber] =
mark(subsegloop);
14723 fprintf(outfile,
"%4ld %4d %4d\n", subsegnumber,
14726 fprintf(outfile,
"%4ld %4d %4d %4d\n", subsegnumber,
14737 fprintf(outfile,
"%d\n", holes);
14739 for (holenumber = 0; holenumber < holes; holenumber++) {
14741 fprintf(outfile,
"%4ld %.17g %.17g\n",
b->firstnumber + holenumber,
14742 holelist[2 * holenumber], holelist[2 * holenumber + 1]);
14746 fprintf(outfile,
"%d\n", regions);
14747 for (regionnumber = 0; regionnumber < regions; regionnumber++) {
14749 fprintf(outfile,
"%4ld %.17g %.17g %.17g %.17g\n",
14750 b->firstnumber + regionnumber,
14751 regionlist[4 * regionnumber], regionlist[4 * regionnumber + 1],
14752 regionlist[4 * regionnumber + 2],
14753 regionlist[4 * regionnumber + 3]);
14758 finishfile(outfile, argc, argv);
14770#ifdef ANSI_DECLARATORS
14772 int **edgelist,
int **edgemarkerlist)
14778int **edgemarkerlist;
14783#ifdef ANSI_DECLARATORS
14785 int argc,
char **argv)
14805 struct otri triangleloop, trisym;
14806 struct osub checkmark;
14814 printf(
"Writing edges.\n");
14817 if (*edgelist == (
int *) NULL) {
14818 *edgelist = (
int *)
trimalloc((
int) (
m->edges * 2 *
sizeof(
int)));
14821 if (!
b->nobound && (*edgemarkerlist == (
int *) NULL)) {
14822 *edgemarkerlist = (
int *)
trimalloc((
int) (
m->edges *
sizeof(
int)));
14825 emlist = *edgemarkerlist;
14829 printf(
"Writing %s.\n", edgefilename);
14831 outfile = fopen(edgefilename,
"w");
14832 if (outfile == (FILE *) NULL) {
14833 printf(
" Error: Cannot create file %s.\n", edgefilename);
14837 fprintf(outfile,
"%ld %d\n",
m->edges, 1 -
b->nobound);
14842 edgenumber =
b->firstnumber;
14850 for (triangleloop.
orient = 0; triangleloop.
orient < 3;
14851 triangleloop.
orient++) {
14852 sym(triangleloop, trisym);
14853 if ((triangleloop.
tri < trisym.
tri) || (trisym.
tri ==
m->dummytri)) {
14854 org(triangleloop, p1);
14855 dest(triangleloop, p2);
14863 fprintf(outfile,
"%4ld %d %d\n", edgenumber,
14869 if (
b->usesegments) {
14870 tspivot(triangleloop, checkmark);
14871 if (checkmark.
ss ==
m->dummysub) {
14873 emlist[edgenumber -
b->firstnumber] = 0;
14875 fprintf(outfile,
"%4ld %d %d %d\n", edgenumber,
14880 emlist[edgenumber -
b->firstnumber] =
mark(checkmark);
14882 fprintf(outfile,
"%4ld %d %d %d\n", edgenumber,
14888 emlist[edgenumber -
b->firstnumber] = trisym.
tri ==
m->dummytri;
14890 fprintf(outfile,
"%4ld %d %d %d\n", edgenumber,
14902 finishfile(outfile, argc, argv);
14924#ifdef ANSI_DECLARATORS
14926 REAL **vpointattriblist,
int **vpointmarkerlist,
14927 int **vedgelist,
int **vedgemarkerlist,
REAL **vnormlist)
14929void writevoronoi(
m,
b, vpointlist, vpointattriblist, vpointmarkerlist,
14930 vedgelist, vedgemarkerlist, vnormlist)
14934REAL **vpointattriblist;
14935int **vpointmarkerlist;
14937int **vedgemarkerlist;
14943#ifdef ANSI_DECLARATORS
14945 char *vedgefilename,
int argc,
char **argv)
14947void writevoronoi(
m,
b, vnodefilename, vedgefilename, argc, argv)
14950char *vnodefilename;
14951char *vedgefilename;
14969 struct otri triangleloop, trisym;
14970 vertex torg, tdest, tapex;
14971 REAL circumcenter[2];
14983 printf(
"Writing Voronoi vertices.\n");
14986 if (*vpointlist == (
REAL *) NULL) {
14987 *vpointlist = (
REAL *)
trimalloc((
int) (
m->triangles.items * 2 *
14991 if (*vpointattriblist == (
REAL *) NULL) {
14992 *vpointattriblist = (
REAL *)
trimalloc((
int) (
m->triangles.items *
14993 m->nextras *
sizeof(
REAL)));
14995 *vpointmarkerlist = (
int *) NULL;
14996 plist = *vpointlist;
14997 palist = *vpointattriblist;
15002 printf(
"Writing %s.\n", vnodefilename);
15004 outfile = fopen(vnodefilename,
"w");
15005 if (outfile == (FILE *) NULL) {
15006 printf(
" Error: Cannot create file %s.\n", vnodefilename);
15011 fprintf(outfile,
"%ld %d %d %d\n",
m->triangles.items, 2,
m->nextras, 0);
15016 triangleloop.
orient = 0;
15017 vnodenumber =
b->firstnumber;
15019 org(triangleloop, torg);
15020 dest(triangleloop, tdest);
15021 apex(triangleloop, tapex);
15025 plist[coordindex++] = circumcenter[0];
15026 plist[coordindex++] = circumcenter[1];
15027 for (i = 2; i < 2 +
m->nextras; i++) {
15029 palist[attribindex++] = torg[i] + xi * (tdest[i] - torg[i])
15030 + eta * (tapex[i] - torg[i]);
15034 fprintf(outfile,
"%4ld %.17g %.17g", vnodenumber, circumcenter[0],
15036 for (i = 2; i < 2 +
m->nextras; i++) {
15038 fprintf(outfile,
" %.17g", torg[i] + xi * (tdest[i] - torg[i])
15039 + eta * (tapex[i] - torg[i]));
15041 fprintf(outfile,
"\n");
15044 * (
int *) (triangleloop.
tri + 6) = (
int) vnodenumber;
15050 finishfile(outfile, argc, argv);
15055 printf(
"Writing Voronoi edges.\n");
15058 if (*vedgelist == (
int *) NULL) {
15059 *vedgelist = (
int *)
trimalloc((
int) (
m->edges * 2 *
sizeof(
int)));
15061 *vedgemarkerlist = (
int *) NULL;
15063 if (*vnormlist == (
REAL *) NULL) {
15066 elist = *vedgelist;
15067 normlist = *vnormlist;
15071 printf(
"Writing %s.\n", vedgefilename);
15073 outfile = fopen(vedgefilename,
"w");
15074 if (outfile == (FILE *) NULL) {
15075 printf(
" Error: Cannot create file %s.\n", vedgefilename);
15079 fprintf(outfile,
"%ld %d\n",
m->edges, 0);
15085 vedgenumber =
b->firstnumber;
15094 for (triangleloop.
orient = 0; triangleloop.
orient < 3;
15095 triangleloop.
orient++) {
15096 sym(triangleloop, trisym);
15097 if ((triangleloop.
tri < trisym.
tri) || (trisym.
tri ==
m->dummytri)) {
15099 p1 = * (
int *) (triangleloop.
tri + 6);
15100 if (trisym.
tri ==
m->dummytri) {
15101 org(triangleloop, torg);
15102 dest(triangleloop, tdest);
15105 elist[coordindex] = p1;
15106 normlist[coordindex++] = tdest[1] - torg[1];
15107 elist[coordindex] = -1;
15108 normlist[coordindex++] = torg[0] - tdest[0];
15113 fprintf(outfile,
"%4ld %d %d %.17g %.17g\n", vedgenumber,
15114 p1, -1, tdest[1] - torg[1], torg[0] - tdest[0]);
15118 p2 = * (
int *) (trisym.
tri + 6);
15121 elist[coordindex] = p1;
15122 normlist[coordindex++] = 0.0;
15123 elist[coordindex] = p2;
15124 normlist[coordindex++] = 0.0;
15126 fprintf(outfile,
"%4ld %d %d\n", vedgenumber, p1, p2);
15138 finishfile(outfile, argc, argv);
15144#ifdef ANSI_DECLARATORS
15155#ifdef ANSI_DECLARATORS
15157 int argc,
char **argv)
15162char *neighborfilename;
15176 struct otri triangleloop, trisym;
15177 long elementnumber;
15178 int neighbor1, neighbor2, neighbor3;
15183 printf(
"Writing neighbors.\n");
15186 if (*neighborlist == (
int *) NULL) {
15187 *neighborlist = (
int *)
trimalloc((
int) (
m->triangles.items * 3 *
15190 nlist = *neighborlist;
15194 printf(
"Writing %s.\n", neighborfilename);
15196 outfile = fopen(neighborfilename,
"w");
15197 if (outfile == (FILE *) NULL) {
15198 printf(
" Error: Cannot create file %s.\n", neighborfilename);
15202 fprintf(outfile,
"%ld %d\n",
m->triangles.items, 3);
15207 triangleloop.
orient = 0;
15208 elementnumber =
b->firstnumber;
15210 * (
int *) (triangleloop.
tri + 6) = (
int) elementnumber;
15214 * (
int *) (
m->dummytri + 6) = -1;
15218 elementnumber =
b->firstnumber;
15220 triangleloop.
orient = 1;
15221 sym(triangleloop, trisym);
15222 neighbor1 = * (
int *) (trisym.
tri + 6);
15223 triangleloop.
orient = 2;
15224 sym(triangleloop, trisym);
15225 neighbor2 = * (
int *) (trisym.
tri + 6);
15226 triangleloop.
orient = 0;
15227 sym(triangleloop, trisym);
15228 neighbor3 = * (
int *) (trisym.
tri + 6);
15230 nlist[
index++] = neighbor1;
15231 nlist[
index++] = neighbor2;
15232 nlist[
index++] = neighbor3;
15235 fprintf(outfile,
"%4ld %d %d %d\n", elementnumber,
15236 neighbor1, neighbor2, neighbor3);
15244 finishfile(outfile, argc, argv);
15259#ifdef ANSI_DECLARATORS
15260void writeoff(
struct mesh *
m,
struct behavior *
b,
char *offfilename,
15261 int argc,
char **argv)
15263void writeoff(
m,
b, offfilename, argc, argv)
15273 struct otri triangleloop;
15279 printf(
"Writing %s.\n", offfilename);
15283 outvertices =
m->vertices.items -
m->undeads;
15285 outvertices =
m->vertices.items;
15288 outfile = fopen(offfilename,
"w");
15289 if (outfile == (FILE *) NULL) {
15290 printf(
" Error: Cannot create file %s.\n", offfilename);
15294 fprintf(outfile,
"OFF\n%ld %ld %ld\n", outvertices,
m->triangles.items,
15300 while (vertexloop != (
vertex) NULL) {
15303 fprintf(outfile,
" %.17g %.17g %.17g\n", vertexloop[0], vertexloop[1],
15312 triangleloop.orient = 0;
15313 while (triangleloop.tri != (
triangle *) NULL) {
15314 org(triangleloop, p1);
15315 dest(triangleloop, p2);
15316 apex(triangleloop, p3);
15318 fprintf(outfile,
" 3 %4d %4d %4d\n",
vertexmark(p1) -
b->firstnumber,
15322 finishfile(outfile, argc, argv);
15337#ifdef ANSI_DECLARATORS
15346 struct otri triangleloop;
15348 REAL cossquaretable[8];
15349 REAL ratiotable[16];
15351 REAL edgelength[3];
15355 REAL shortest, longest;
15357 REAL smallestarea, biggestarea;
15358 REAL triminaltitude2;
15362 REAL smallestangle, biggestangle;
15363 REAL radconst, degconst;
15364 int angletable[18];
15365 int aspecttable[16];
15371 printf(
"Mesh quality statistics:\n\n");
15372 radconst =
PI / 18.0;
15373 degconst = 180.0 /
PI;
15374 for (i = 0; i < 8; i++) {
15375 cossquaretable[i] = cos(radconst * (
REAL) (i + 1));
15376 cossquaretable[i] = cossquaretable[i] * cossquaretable[i];
15378 for (i = 0; i < 18; i++) {
15382 ratiotable[0] = 1.5; ratiotable[1] = 2.0;
15383 ratiotable[2] = 2.5; ratiotable[3] = 3.0;
15384 ratiotable[4] = 4.0; ratiotable[5] = 6.0;
15385 ratiotable[6] = 10.0; ratiotable[7] = 15.0;
15386 ratiotable[8] = 25.0; ratiotable[9] = 50.0;
15387 ratiotable[10] = 100.0; ratiotable[11] = 300.0;
15388 ratiotable[12] = 1000.0; ratiotable[13] = 10000.0;
15389 ratiotable[14] = 100000.0; ratiotable[15] = 0.0;
15390 for (i = 0; i < 16; i++) {
15391 aspecttable[i] = 0;
15395 minaltitude =
m->xmax -
m->xmin +
m->ymax -
m->ymin;
15396 minaltitude = minaltitude * minaltitude;
15397 shortest = minaltitude;
15399 smallestarea = minaltitude;
15402 smallestangle = 0.0;
15403 biggestangle = 2.0;
15408 triangleloop.
orient = 0;
15410 org(triangleloop,
p[0]);
15411 dest(triangleloop,
p[1]);
15412 apex(triangleloop,
p[2]);
15415 for (i = 0; i < 3; i++) {
15418 dx[i] =
p[j][0] -
p[k][0];
15419 dy[i] =
p[j][1] -
p[k][1];
15420 edgelength[i] = dx[i] * dx[i] + dy[i] * dy[i];
15421 if (edgelength[i] > trilongest2) {
15422 trilongest2 = edgelength[i];
15424 if (edgelength[i] > longest) {
15425 longest = edgelength[i];
15427 if (edgelength[i] < shortest) {
15428 shortest = edgelength[i];
15433 if (triarea < smallestarea) {
15434 smallestarea = triarea;
15436 if (triarea > biggestarea) {
15437 biggestarea = triarea;
15439 triminaltitude2 = triarea * triarea / trilongest2;
15440 if (triminaltitude2 < minaltitude) {
15441 minaltitude = triminaltitude2;
15443 triaspect2 = trilongest2 / triminaltitude2;
15444 if (triaspect2 > worstaspect) {
15445 worstaspect = triaspect2;
15448 while ((aspectindex < 15) &&
15449 (triaspect2 > ratiotable[aspectindex] * ratiotable[aspectindex])
15453 aspecttable[aspectindex]++;
15455 for (i = 0; i < 3; i++) {
15458 dotproduct = dx[j] * dx[k] + dy[j] * dy[k];
15459 cossquare = dotproduct * dotproduct / (edgelength[j] * edgelength[k]);
15461 for (ii = 7; ii >= 0; ii--) {
15462 if (cossquare > cossquaretable[ii]) {
15466 if (dotproduct <= 0.0) {
15467 angletable[tendegree]++;
15468 if (cossquare > smallestangle) {
15469 smallestangle = cossquare;
15471 if (acutebiggest && (cossquare < biggestangle)) {
15472 biggestangle = cossquare;
15475 angletable[17 - tendegree]++;
15476 if (acutebiggest || (cossquare > biggestangle)) {
15477 biggestangle = cossquare;
15485 shortest = sqrt(shortest);
15486 longest = sqrt(longest);
15487 minaltitude = sqrt(minaltitude);
15488 worstaspect = sqrt(worstaspect);
15489 smallestarea *= 0.5;
15490 biggestarea *= 0.5;
15491 if (smallestangle >= 1.0) {
15492 smallestangle = 0.0;
15494 smallestangle = degconst * acos(sqrt(smallestangle));
15496 if (biggestangle >= 1.0) {
15497 biggestangle = 180.0;
15499 if (acutebiggest) {
15500 biggestangle = degconst * acos(sqrt(biggestangle));
15502 biggestangle = 180.0 - degconst * acos(sqrt(biggestangle));
15506 printf(
" Smallest area: %16.5g | Largest area: %16.5g\n",
15507 smallestarea, biggestarea);
15508 printf(
" Shortest edge: %16.5g | Longest edge: %16.5g\n",
15509 shortest, longest);
15510 printf(
" Shortest altitude: %12.5g | Largest aspect ratio: %8.5g\n\n",
15511 minaltitude, worstaspect);
15513 printf(
" Triangle aspect ratio histogram:\n");
15514 printf(
" 1.1547 - %-6.6g : %8d | %6.6g - %-6.6g : %8d\n",
15515 ratiotable[0], aspecttable[0], ratiotable[7], ratiotable[8],
15517 for (i = 1; i < 7; i++) {
15518 printf(
" %6.6g - %-6.6g : %8d | %6.6g - %-6.6g : %8d\n",
15519 ratiotable[i - 1], ratiotable[i], aspecttable[i],
15520 ratiotable[i + 7], ratiotable[i + 8], aspecttable[i + 8]);
15522 printf(
" %6.6g - %-6.6g : %8d | %6.6g - : %8d\n",
15523 ratiotable[6], ratiotable[7], aspecttable[7], ratiotable[14],
15525 printf(
" (Aspect ratio is longest edge divided by shortest altitude)\n\n");
15527 printf(
" Smallest angle: %15.5g | Largest angle: %15.5g\n\n",
15528 smallestangle, biggestangle);
15530 printf(
" Angle histogram:\n");
15531 for (i = 0; i < 9; i++) {
15532 printf(
" %3d - %3d degrees: %8d | %3d - %3d degrees: %8d\n",
15533 i * 10, i * 10 + 10, angletable[i],
15534 i * 10 + 90, i * 10 + 100, angletable[i + 9]);
15545#ifdef ANSI_DECLARATORS
15554 printf(
"\nStatistics:\n\n");
15555 printf(
" Input vertices: %d\n",
m->invertices);
15557 printf(
" Input triangles: %d\n",
m->inelements);
15560 printf(
" Input segments: %d\n",
m->insegments);
15562 printf(
" Input holes: %d\n",
m->holes);
15566 printf(
"\n Mesh vertices: %ld\n",
m->vertices.items -
m->undeads);
15567 printf(
" Mesh triangles: %ld\n",
m->triangles.items);
15568 printf(
" Mesh edges: %ld\n",
m->edges);
15569 printf(
" Mesh exterior boundary edges: %ld\n",
m->hullsize);
15570 if (
b->poly ||
b->refine) {
15571 printf(
" Mesh interior boundary edges: %ld\n",
15572 m->subsegs.items -
m->hullsize);
15573 printf(
" Mesh subsegments (constrained edges): %ld\n",
15580 printf(
"Memory allocation statistics:\n\n");
15581 printf(
" Maximum number of vertices: %ld\n",
m->vertices.maxitems);
15582 printf(
" Maximum number of triangles: %ld\n",
m->triangles.maxitems);
15583 if (
m->subsegs.maxitems > 0) {
15584 printf(
" Maximum number of subsegments: %ld\n",
m->subsegs.maxitems);
15586 if (
m->viri.maxitems > 0) {
15587 printf(
" Maximum number of viri: %ld\n",
m->viri.maxitems);
15589 if (
m->badsubsegs.maxitems > 0) {
15590 printf(
" Maximum number of encroached subsegments: %ld\n",
15591 m->badsubsegs.maxitems);
15593 if (
m->badtriangles.maxitems > 0) {
15594 printf(
" Maximum number of bad triangles: %ld\n",
15595 m->badtriangles.maxitems);
15597 if (
m->flipstackers.maxitems > 0) {
15598 printf(
" Maximum number of stacked triangle flips: %ld\n",
15599 m->flipstackers.maxitems);
15601 if (
m->splaynodes.maxitems > 0) {
15602 printf(
" Maximum number of splay tree nodes: %ld\n",
15603 m->splaynodes.maxitems);
15605 printf(
" Approximate heap memory use (bytes): %ld\n\n",
15606 m->vertices.maxitems *
m->vertices.itembytes +
15607 m->triangles.maxitems *
m->triangles.itembytes +
15608 m->subsegs.maxitems *
m->subsegs.itembytes +
15609 m->viri.maxitems *
m->viri.itembytes +
15610 m->badsubsegs.maxitems *
m->badsubsegs.itembytes +
15611 m->badtriangles.maxitems *
m->badtriangles.itembytes +
15612 m->flipstackers.maxitems *
m->flipstackers.itembytes +
15613 m->splaynodes.maxitems *
m->splaynodes.itembytes);
15615 printf(
"Algorithmic statistics:\n\n");
15616 if (!
b->weighted) {
15617 printf(
" Number of incircle tests: %ld\n",
m->incirclecount);
15619 printf(
" Number of 3D orientation tests: %ld\n",
m->orient3dcount);
15621 printf(
" Number of 2D orientation tests: %ld\n",
m->counterclockcount);
15622 if (
m->hyperbolacount > 0) {
15623 printf(
" Number of right-of-hyperbola tests: %ld\n",
15624 m->hyperbolacount);
15626 if (
m->circletopcount > 0) {
15627 printf(
" Number of circle top computations: %ld\n",
15628 m->circletopcount);
15630 if (
m->circumcentercount > 0) {
15631 printf(
" Number of triangle circumcenter computations: %ld\n",
15632 m->circumcentercount);
15665#ifdef ANSI_DECLARATORS
15678#ifdef ANSI_DECLARATORS
15679int main(
int argc,
char **argv)
15681int main(argc, argv)
15699 struct timeval tv0,
tv1,
tv2,
tv3, tv4, tv5, tv6;
15700 struct timezone tz;
15704 gettimeofday(&tv0, &tz);
15713 m.steinerleft =
b.steiner;
15720 readnodes(&
m, &
b,
b.innodefilename,
b.inpolyfilename, &polyfile);
15725 gettimeofday(&
tv1, &tz);
15742 m.hullsize = reconstruct(&
m, &
b,
b.inelefilename,
b.areafilename,
15743 b.inpolyfilename, polyfile);
15752 gettimeofday(&
tv2, &tz);
15754 printf(
"Mesh reconstruction");
15756 printf(
"Delaunay");
15758 printf(
" milliseconds: %ld\n", 1000l * (
tv2.tv_sec -
tv1.tv_sec) +
15759 (
tv2.tv_usec -
tv1.tv_usec) / 1000l);
15765 m.infvertex1 = (
vertex) NULL;
15766 m.infvertex2 = (
vertex) NULL;
15767 m.infvertex3 = (
vertex) NULL;
15769 if (
b.usesegments) {
15770 m.checksegments = 1;
15784 gettimeofday(&
tv3, &tz);
15785 if (
b.usesegments && !
b.refine) {
15786 printf(
"Segment milliseconds: %ld\n",
15787 1000l * (
tv3.tv_sec -
tv2.tv_sec) +
15788 (
tv3.tv_usec -
tv2.tv_usec) / 1000l);
15793 if (
b.poly && (
m.triangles.items > 0)) {
15800 readholes(&
m, &
b, polyfile,
b.inpolyfilename, &holearray, &
m.holes,
15801 ®ionarray, &
m.regions);
15817 gettimeofday(&tv4, &tz);
15818 if (
b.poly && !
b.refine) {
15819 printf(
"Hole milliseconds: %ld\n", 1000l * (tv4.tv_sec -
tv3.tv_sec) +
15820 (tv4.tv_usec -
tv3.tv_usec) / 1000l);
15826 if (
b.quality && (
m.triangles.items > 0)) {
15827 enforcequality(&
m, &
b);
15833 gettimeofday(&tv5, &tz);
15836 printf(
"Quality milliseconds: %ld\n",
15837 1000l * (tv5.tv_sec - tv4.tv_sec) +
15838 (tv5.tv_usec - tv4.tv_usec) / 1000l);
15845 m.edges = (3l *
m.triangles.items +
m.hullsize) / 2l;
15856 out->numberofpoints =
m.vertices.items -
m.undeads;
15858 out->numberofpoints =
m.vertices.items;
15860 out->numberofpointattributes =
m.nextras;
15861 out->numberoftriangles =
m.triangles.items;
15862 out->numberofcorners = (
b.order + 1) * (
b.order + 2) / 2;
15863 out->numberoftriangleattributes =
m.eextras;
15864 out->numberofedges =
m.edges;
15865 if (
b.usesegments) {
15866 out->numberofsegments =
m.subsegs.items;
15868 out->numberofsegments =
m.hullsize;
15878 if (
b.nonodewritten || (
b.noiterationnum &&
m.readnodefile)) {
15881 printf(
"NOT writing vertices.\n");
15883 printf(
"NOT writing a .node file.\n");
15890 writenodes(&
m, &
b, &out->pointlist, &out->pointattributelist,
15891 &out->pointmarkerlist);
15896 if (
b.noelewritten) {
15899 printf(
"NOT writing triangles.\n");
15901 printf(
"NOT writing an .ele file.\n");
15906 writeelements(&
m, &
b, &out->trianglelist, &out->triangleattributelist);
15913 if (
b.poly ||
b.convex) {
15915 if (
b.nopolywritten ||
b.noiterationnum) {
15918 printf(
"NOT writing segments.\n");
15920 printf(
"NOT writing a .poly file.\n");
15925 writepoly(&
m, &
b, &out->segmentlist, &out->segmentmarkerlist);
15926 out->numberofholes =
m.holes;
15927 out->numberofregions =
m.regions;
15932 out->holelist = (
REAL *) NULL;
15933 out->regionlist = (
REAL *) NULL;
15936 writepoly(&
m, &
b,
b.outpolyfilename, holearray,
m.holes, regionarray,
15937 m.regions, argc, argv);
15943 if (
m.regions > 0) {
15951 writeoff(&
m, &
b,
b.offfilename, argc, argv);
15956 writeedges(&
m, &
b, &out->edgelist, &out->edgemarkerlist);
15980 gettimeofday(&tv6, &tz);
15981 printf(
"\nOutput milliseconds: %ld\n",
15982 1000l * (tv6.tv_sec - tv5.tv_sec) +
15983 (tv6.tv_usec - tv5.tv_usec) / 1000l);
15984 printf(
"Total running milliseconds: %ld\n",
15985 1000l * (tv6.tv_sec - tv0.tv_sec) +
15986 (tv6.tv_usec - tv0.tv_usec) / 1000l);
15995 checkdelaunay(&
m, &
b);
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t result
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t index
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h length
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t child
Option_t Option_t TPoint TPoint angle
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
struct badtriang * nexttriang
struct flipstacker * prevflip
struct memorypool badtriangles
struct badtriang * queuetail[4096]
struct memorypool flipstackers
struct memorypool subsegs
struct memorypool triangles
struct memorypool splaynodes
struct flipstacker * lastflip
struct memorypool vertices
struct badtriang * queuefront[4096]
struct memorypool badsubsegs
struct splaynode * lchild
struct splaynode * rchild
int numberoftriangleattributes
double * pointattributelist
double * trianglearealist
int numberofpointattributes
double * triangleattributelist
static uint64_t sum(uint64_t i)
#define sdest(osub, vertexptr)
#define sorg(osub, vertexptr)
void highorder(struct mesh *m, struct behavior *b)
#define apex(otri, vertexptr)
void makevertexmap(struct mesh *m, struct behavior *b)
#define dnext(otri1, otri2)
void numbernodes(struct mesh *m, struct behavior *b)
void writepoly(struct mesh *m, struct behavior *b, int **segmentlist, int **segmentmarkerlist)
vertex getvertex(struct mesh *m, struct behavior *b, int number)
#define setelemattribute(otri, attnum, value)
void * trimalloc(int size)
double counterclockwise(struct mesh *m, struct behavior *b, vertex pa, vertex pb, vertex pc)
#define setsorg(osub, vertexptr)
#define segdest(osub, vertexptr)
void poolrestart(struct memorypool *pool)
#define segorg(osub, vertexptr)
void * traverse(struct memorypool *pool)
#define Two_Product(a, b, x, y)
void vertexmedian(vertex *sortarray, int arraysize, int median, int axis)
void * poolalloc(struct memorypool *pool)
#define sdecode(sptr, osub)
enum insertvertexresult insertvertex(struct mesh *m, struct behavior *b, vertex newvertex, struct otri *searchtri, struct osub *splitseg, int segmentflaws, int triflaws)
#define setorg(otri, vertexptr)
enum locateresult preciselocate(struct mesh *m, struct behavior *b, vertex searchpoint, struct otri *searchtri, int stopatsubsegment)
#define elemattribute(otri, attnum)
void makesubseg(struct mesh *m, struct osub *newsubseg)
#define bond(otri1, otri2)
void triangledealloc(struct mesh *m, triangle *dyingtriangle)
#define setsegorg(osub, vertexptr)
void alternateaxes(vertex *sortarray, int arraysize, int axis)
#define lnext(otri1, otri2)
void writeneighbors(struct mesh *m, struct behavior *b, int **neighborlist)
void writevoronoi(struct mesh *m, struct behavior *b, double **vpointlist, double **vpointattriblist, int **vpointmarkerlist, int **vedgelist, int **vedgemarkerlist, double **vnormlist)
void regionplague(struct mesh *m, struct behavior *b, double attribute, double area)
void segmentintersection(struct mesh *m, struct behavior *b, struct otri *splittri, struct osub *splitsubseg, vertex endpoint2)
long divconqdelaunay(struct mesh *m, struct behavior *b)
#define Two_Product_Presplit(a, b, bhi, blo, x, y)
#define onext(otri1, otri2)
void delaunayfixup(struct mesh *m, struct behavior *b, struct otri *fixuptri, int leftside)
void triangledeinit(struct mesh *m, struct behavior *b)
void poolinit(struct memorypool *pool, int bytecount, int itemcount, int firstitemcount, int alignment)
enum locateresult locate(struct mesh *m, struct behavior *b, vertex searchpoint, struct otri *searchtri)
vertex vertextraverse(struct mesh *m)
#define lprev(otri1, otri2)
triangle * triangletraverse(struct mesh *m)
#define tspivot(otri, osub)
double nonregular(struct mesh *m, struct behavior *b, vertex pa, vertex pb, vertex pc, vertex pd)
double estimate(int elen, double *e)
double incircleadapt(vertex pa, vertex pb, vertex pc, vertex pd, double permanent)
void maketriangle(struct mesh *m, struct behavior *b, struct otri *newotri)
void printtriangle(struct mesh *m, struct behavior *b, struct otri *t)
#define Two_Diff_Tail(a, b, x, y)
void poolzero(struct memorypool *pool)
void insertsubseg(struct mesh *m, struct behavior *b, struct otri *tri, int subsegmark)
#define dprev(otri1, otri2)
void formskeleton(struct mesh *m, struct behavior *b, int *segmentlist, int *segmentmarkerlist, int numberofsegments)
#define setvertexmark(vx, value)
double incircle(struct mesh *m, struct behavior *b, vertex pa, vertex pb, vertex pc, vertex pd)
#define SPLAYNODEPERBLOCK
void triangulate(char *triswitches, struct triangulateio *in, struct triangulateio *out, struct triangulateio *vorout)
#define setvertex2tri(vx, value)
void initializevertexpool(struct mesh *m, struct behavior *b)
int triunsuitable(vertex triorg, vertex tridest, vertex triapex, double area)
void findcircumcenter(struct mesh *m, struct behavior *b, vertex torg, vertex tdest, vertex tapex, vertex circumcenter, double *xi, double *eta, int offcenter)
void pooldealloc(struct memorypool *pool, void *dyingitem)
void triangulatepolygon(struct mesh *m, struct behavior *b, struct otri *firstedge, struct otri *lastedge, int edgecount, int doflip, int triflaws)
void mergehulls(struct mesh *m, struct behavior *b, struct otri *farleft, struct otri *innerleft, struct otri *innerright, struct otri *farright, int axis)
void parsecommandline(int argc, char **argv, struct behavior *b)
double counterclockwiseadapt(vertex pa, vertex pb, vertex pc, double detsum)
#define oprev(otri1, otri2)
#define BADSUBSEGPERBLOCK
void subsegdealloc(struct mesh *m, subseg *dyingsubseg)
long removeghosts(struct mesh *m, struct behavior *b, struct otri *startghost)
void vertexdealloc(struct mesh *m, vertex dyingvertex)
void infecthull(struct mesh *m, struct behavior *b)
void writenodes(struct mesh *m, struct behavior *b, double **pointlist, double **pointattriblist, int **pointmarkerlist)
#define setapex(otri, vertexptr)
double orient3dadapt(vertex pa, vertex pb, vertex pc, vertex pd, double aheight, double bheight, double cheight, double dheight, double permanent)
#define dest(otri, vertexptr)
#define otricopy(otri1, otri2)
#define Two_Sum(a, b, x, y)
#define tsbond(otri, osub)
#define setmark(osub, value)
void initializetrisubpools(struct mesh *m, struct behavior *b)
double orient3d(struct mesh *m, struct behavior *b, vertex pa, vertex pb, vertex pc, vertex pd, double aheight, double bheight, double cheight, double dheight)
void carveholes(struct mesh *m, struct behavior *b, double *holelist, int holes, double *regionlist, int regions)
void writeedges(struct mesh *m, struct behavior *b, int **edgelist, int **edgemarkerlist)
#define Two_One_Product(a1, a0, b, x3, x2, x1, x0)
#define Two_Two_Diff(a1, a0, b1, b0, x3, x2, x1, x0)
long delaunay(struct mesh *m, struct behavior *b)
#define setdest(otri, vertexptr)
void markhull(struct mesh *m, struct behavior *b)
#define Split(a, ahi, aLo)
void vertexsort(vertex *sortarray, int arraysize)
#define stpivot(osub, otri)
#define decode(ptr, otri)
#define otriequal(otri1, otri2)
#define org(otri, vertexptr)
void traversalinit(struct memorypool *pool)
void plague(struct mesh *m, struct behavior *b)
#define setsegdest(osub, vertexptr)
#define FLIPSTACKERPERBLOCK
int scale_expansion_zeroelim(int elen, double *e, double b, double *h)
void unflip(struct mesh *m, struct behavior *b, struct otri *flipedge)
#define setvertextype(vx, value)
uintptr_t randomnation(unsigned int choices)
void dummyinit(struct mesh *m, struct behavior *b, int trianglebytes, int subsegbytes)
#define sym(otri1, otri2)
void quality_statistics(struct mesh *m, struct behavior *b)
void trifree(void *memptr)
void flip(struct mesh *m, struct behavior *b, struct otri *flipedge)
void pooldeinit(struct memorypool *pool)
enum finddirectionresult finddirection(struct mesh *m, struct behavior *b, struct otri *searchtri, vertex searchpoint)
void constrainededge(struct mesh *m, struct behavior *b, struct otri *starttri, vertex endpoint2, int newmark)
#define setsdest(osub, vertexptr)
void transfernodes(struct mesh *m, struct behavior *b, double *pointlist, double *pointattriblist, int *pointmarkerlist, int numberofpoints, int numberofpointattribs)
#define spivot(osub1, osub2)
#define ssym(osub1, osub2)
void insertsegment(struct mesh *m, struct behavior *b, vertex endpoint1, vertex endpoint2, int newmark)
#define Two_Two_Sum(a1, a0, b1, b0, x3, x2, x1, x0)
void printsubseg(struct mesh *m, struct behavior *b, struct osub *s)
void divconqrecurse(struct mesh *m, struct behavior *b, vertex *sortarray, int vertices, int axis, struct otri *farleft, struct otri *farright)
#define Fast_Two_Sum(a, b, x, y)
int scoutsegment(struct mesh *m, struct behavior *b, struct otri *searchtri, vertex endpoint2, int newmark)
void writeelements(struct mesh *m, struct behavior *b, int **trianglelist, double **triangleattriblist)
subseg * subsegtraverse(struct mesh *m)
void statistics(struct mesh *m, struct behavior *b)
int fast_expansion_sum_zeroelim(int elen, double *e, int flen, double *f, double *h)
#define sbond(osub1, osub2)
#define setareabound(otri, value)
void triangleinit(struct mesh *m)