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);