284 #define FILENAMESIZE 2048
289 #define INPUTLINESIZE 1024
295 #define TRIPERBLOCK 4092
296 #define SUBSEGPERBLOCK 508
297 #define VERTEXPERBLOCK 4092
298 #define VIRUSPERBLOCK 1020
300 #define BADSUBSEGPERBLOCK 252
302 #define BADTRIPERBLOCK 4092
304 #define FLIPSTACKERPERBLOCK 252
306 #define SPLAYNODEPERBLOCK 508
312 #define INPUTVERTEX 0
313 #define SEGMENTVERTEX 1
315 #define DEADVERTEX -32768
316 #define UNDEADVERTEX -32767
336 #define SAMPLEFACTOR 11
342 #define SAMPLERATE 10
346 #define PI 3.141592653589793238462643383279502884197169399375105820974944592308
350 #define SQUAREROOTTWO 1.4142135623730950488016887242096980785696718753769480732
354 #define ONETHIRD 0.333333333333333333333333333333333333333333333333333333333333
363 #include <WinSock2.h>
366 #if defined(_MSC_EXTENSIONS)
367 #define DELTA_EPOCH_IN_MICROSECS 11644473600000000Ui64
369 #define DELTA_EPOCH_IN_MICROSECS 11644473600000000ULL
380 int gettimeofday(
struct timeval *tv,
struct timezone *tz)
386 unsigned __int64 tmpres = 0;
387 static int tzflag = 0;
391 GetSystemTimeAsFileTime(&ft);
396 tmpres |= ft.dwHighDateTime;
398 tmpres |= ft.dwLowDateTime;
405 tmpres -= DELTA_EPOCH_IN_MICROSECS;
409 tv->tv_sec = (long)(tmpres / 1000000UL);
410 tv->tv_usec = (long)(tmpres % 1000000UL);
422 tz->tz_minuteswest = _timezone / 60;
423 tz->tz_dsttime = _daylight;
428 #include <sys/time.h>
435 #include <fpu_control.h>
438 #include "triangle.h"
452 enum locateresult {INTRIANGLE, ONEDGE, ONVERTEX, OUTSIDE};
460 enum insertvertexresult {SUCCESSFULVERTEX, ENCROACHINGVERTEX, VIOLATINGVERTEX,
468 enum finddirectionresult {WITHIN, LEFTCOLLINEAR, RIGHTCOLLINEAR};
587 typedef REAL **triangle;
604 typedef REAL **subseg;
621 typedef REAL *vertex;
628 vertex subsegorg, subsegdest;
638 vertex triangorg, triangdest, triangapex;
708 VOID **firstblock, **nowblock;
717 long items, maxitems;
718 int unallocateditems;
728 REAL ccwerrboundA, ccwerrboundB, ccwerrboundC;
729 REAL iccerrboundA, iccerrboundB, iccerrboundC;
730 REAL o3derrboundA, o3derrboundB, o3derrboundC;
734 unsigned long randomseed;
760 int nextnonemptyq[4096];
769 REAL xmin, xmax, ymin, ymax;
794 long counterclockcount;
797 long circumcentercount;
802 vertex infvertex1, infvertex2, infvertex3;
807 triangle *dummytribase;
814 subseg *dummysubbase;
819 struct otri recenttri;
862 int poly, refine, quality, vararea, fixedarea, usertest;
863 int regionattrib, convex, weighted, jettison;
865 int edgesout, voronoi, neighbors, geomview;
866 int nobound, nopolywritten, nonodewritten, noelewritten, noiterationnum;
867 int noholes, noexact, conformdel;
868 int incremental, sweepline, dwyer;
876 REAL minangle, goodangle, offconstant;
882 char innodefilename[FILENAMESIZE];
883 char inelefilename[FILENAMESIZE];
884 char inpolyfilename[FILENAMESIZE];
885 char areafilename[FILENAMESIZE];
886 char outnodefilename[FILENAMESIZE];
887 char outelefilename[FILENAMESIZE];
888 char outpolyfilename[FILENAMESIZE];
889 char edgefilename[FILENAMESIZE];
890 char vnodefilename[FILENAMESIZE];
891 char vedgefilename[FILENAMESIZE];
892 char neighborfilename[FILENAMESIZE];
893 char offfilename[FILENAMESIZE];
1011 int plus1mod3[3] = {1, 2, 0};
1012 int minus1mod3[3] = {2, 0, 1};
1021 #define decode(ptr, otri) \
1022 (otri).orient = (int) ((unsigned long) (ptr) & (unsigned long) 3l); \
1023 (otri).tri = (triangle *) \
1024 ((unsigned long) (ptr) ^ (unsigned long) (otri).orient)
1030 #define encode(otri) \
1031 (triangle) ((unsigned long) (otri).tri | (unsigned long) (otri).orient)
1041 #define sym(otri1, otri2) \
1042 ptr = (otri1).tri[(otri1).orient]; \
1045 #define symself(otri) \
1046 ptr = (otri).tri[(otri).orient]; \
1051 #define lnext(otri1, otri2) \
1052 (otri2).tri = (otri1).tri; \
1053 (otri2).orient = plus1mod3[(otri1).orient]
1055 #define lnextself(otri) \
1056 (otri).orient = plus1mod3[(otri).orient]
1060 #define lprev(otri1, otri2) \
1061 (otri2).tri = (otri1).tri; \
1062 (otri2).orient = minus1mod3[(otri1).orient]
1064 #define lprevself(otri) \
1065 (otri).orient = minus1mod3[(otri).orient]
1071 #define onext(otri1, otri2) \
1072 lprev(otri1, otri2); \
1075 #define onextself(otri) \
1083 #define oprev(otri1, otri2) \
1084 sym(otri1, otri2); \
1087 #define oprevself(otri) \
1095 #define dnext(otri1, otri2) \
1096 sym(otri1, otri2); \
1099 #define dnextself(otri) \
1107 #define dprev(otri1, otri2) \
1108 lnext(otri1, otri2); \
1111 #define dprevself(otri) \
1119 #define rnext(otri1, otri2) \
1120 sym(otri1, otri2); \
1124 #define rnextself(otri) \
1133 #define rprev(otri1, otri2) \
1134 sym(otri1, otri2); \
1138 #define rprevself(otri) \
1146 #define org(otri, vertexptr) \
1147 vertexptr = (vertex) (otri).tri[plus1mod3[(otri).orient] + 3]
1149 #define dest(otri, vertexptr) \
1150 vertexptr = (vertex) (otri).tri[minus1mod3[(otri).orient] + 3]
1152 #define apex(otri, vertexptr) \
1153 vertexptr = (vertex) (otri).tri[(otri).orient + 3]
1155 #define setorg(otri, vertexptr) \
1156 (otri).tri[plus1mod3[(otri).orient] + 3] = (triangle) vertexptr
1158 #define setdest(otri, vertexptr) \
1159 (otri).tri[minus1mod3[(otri).orient] + 3] = (triangle) vertexptr
1161 #define setapex(otri, vertexptr) \
1162 (otri).tri[(otri).orient + 3] = (triangle) vertexptr
1166 #define bond(otri1, otri2) \
1167 (otri1).tri[(otri1).orient] = encode(otri2); \
1168 (otri2).tri[(otri2).orient] = encode(otri1)
1175 #define dissolve(otri) \
1176 (otri).tri[(otri).orient] = (triangle) m->dummytri
1180 #define otricopy(otri1, otri2) \
1181 (otri2).tri = (otri1).tri; \
1182 (otri2).orient = (otri1).orient
1186 #define otriequal(otri1, otri2) \
1187 (((otri1).tri == (otri2).tri) && \
1188 ((otri1).orient == (otri2).orient))
1193 #define infect(otri) \
1194 (otri).tri[6] = (triangle) \
1195 ((unsigned long) (otri).tri[6] | (unsigned long) 2l)
1197 #define uninfect(otri) \
1198 (otri).tri[6] = (triangle) \
1199 ((unsigned long) (otri).tri[6] & ~ (unsigned long) 2l)
1203 #define infected(otri) \
1204 (((unsigned long) (otri).tri[6] & (unsigned long) 2l) != 0l)
1208 #define elemattribute(otri, attnum) \
1209 ((REAL *) (otri).tri)[m->elemattribindex + (attnum)]
1211 #define setelemattribute(otri, attnum, value) \
1212 ((REAL *) (otri).tri)[m->elemattribindex + (attnum)] = value
1216 #define areabound(otri) ((REAL *) (otri).tri)[m->areaboundindex]
1218 #define setareabound(otri, value) \
1219 ((REAL *) (otri).tri)[m->areaboundindex] = value
1226 #define deadtri(tria) ((tria)[1] == (triangle) NULL)
1228 #define killtri(tria) \
1229 (tria)[1] = (triangle) NULL; \
1230 (tria)[3] = (triangle) NULL
1241 #define sdecode(sptr, osub) \
1242 (osub).ssorient = (int) ((unsigned long) (sptr) & (unsigned long) 1l); \
1243 (osub).ss = (subseg *) \
1244 ((unsigned long) (sptr) & ~ (unsigned long) 3l)
1250 #define sencode(osub) \
1251 (subseg) ((unsigned long) (osub).ss | (unsigned long) (osub).ssorient)
1255 #define ssym(osub1, osub2) \
1256 (osub2).ss = (osub1).ss; \
1257 (osub2).ssorient = 1 - (osub1).ssorient
1259 #define ssymself(osub) \
1260 (osub).ssorient = 1 - (osub).ssorient
1265 #define spivot(osub1, osub2) \
1266 sptr = (osub1).ss[(osub1).ssorient]; \
1267 sdecode(sptr, osub2)
1269 #define spivotself(osub) \
1270 sptr = (osub).ss[(osub).ssorient]; \
1276 #define snext(osub1, osub2) \
1277 sptr = (osub1).ss[1 - (osub1).ssorient]; \
1278 sdecode(sptr, osub2)
1280 #define snextself(osub) \
1281 sptr = (osub).ss[1 - (osub).ssorient]; \
1287 #define sorg(osub, vertexptr) \
1288 vertexptr = (vertex) (osub).ss[2 + (osub).ssorient]
1290 #define sdest(osub, vertexptr) \
1291 vertexptr = (vertex) (osub).ss[3 - (osub).ssorient]
1293 #define setsorg(osub, vertexptr) \
1294 (osub).ss[2 + (osub).ssorient] = (subseg) vertexptr
1296 #define setsdest(osub, vertexptr) \
1297 (osub).ss[3 - (osub).ssorient] = (subseg) vertexptr
1299 #define segorg(osub, vertexptr) \
1300 vertexptr = (vertex) (osub).ss[4 + (osub).ssorient]
1302 #define segdest(osub, vertexptr) \
1303 vertexptr = (vertex) (osub).ss[5 - (osub).ssorient]
1305 #define setsegorg(osub, vertexptr) \
1306 (osub).ss[4 + (osub).ssorient] = (subseg) vertexptr
1308 #define setsegdest(osub, vertexptr) \
1309 (osub).ss[5 - (osub).ssorient] = (subseg) vertexptr
1315 #define mark(osub) (* (int *) ((osub).ss + 8))
1317 #define setmark(osub, value) \
1318 * (int *) ((osub).ss + 8) = value
1322 #define sbond(osub1, osub2) \
1323 (osub1).ss[(osub1).ssorient] = sencode(osub2); \
1324 (osub2).ss[(osub2).ssorient] = sencode(osub1)
1329 #define sdissolve(osub) \
1330 (osub).ss[(osub).ssorient] = (subseg) m->dummysub
1334 #define subsegcopy(osub1, osub2) \
1335 (osub2).ss = (osub1).ss; \
1336 (osub2).ssorient = (osub1).ssorient
1340 #define subsegequal(osub1, osub2) \
1341 (((osub1).ss == (osub2).ss) && \
1342 ((osub1).ssorient == (osub2).ssorient))
1349 #define deadsubseg(sub) ((sub)[1] == (subseg) NULL)
1351 #define killsubseg(sub) \
1352 (sub)[1] = (subseg) NULL; \
1353 (sub)[2] = (subseg) NULL
1361 #define tspivot(otri, osub) \
1362 sptr = (subseg) (otri).tri[6 + (otri).orient]; \
1368 #define stpivot(osub, otri) \
1369 ptr = (triangle) (osub).ss[6 + (osub).ssorient]; \
1374 #define tsbond(otri, osub) \
1375 (otri).tri[6 + (otri).orient] = (triangle) sencode(osub); \
1376 (osub).ss[6 + (osub).ssorient] = (subseg) encode(otri)
1380 #define tsdissolve(otri) \
1381 (otri).tri[6 + (otri).orient] = (triangle) m->dummysub
1385 #define stdissolve(osub) \
1386 (osub).ss[6 + (osub).ssorient] = (subseg) m->dummytri
1392 #define vertexmark(vx) ((int *) (vx))[m->vertexmarkindex]
1394 #define setvertexmark(vx, value) \
1395 ((int *) (vx))[m->vertexmarkindex] = value
1397 #define vertextype(vx) ((int *) (vx))[m->vertexmarkindex + 1]
1399 #define setvertextype(vx, value) \
1400 ((int *) (vx))[m->vertexmarkindex + 1] = value
1402 #define vertex2tri(vx) ((triangle *) (vx))[m->vertex2triindex]
1404 #define setvertex2tri(vx, value) \
1405 ((triangle *) (vx))[m->vertex2triindex] = value
1438 #ifdef EXTERNAL_TEST
1440 int triunsuitable();
1444 #ifdef ANSI_DECLARATORS
1445 int triunsuitable(vertex triorg, vertex tridest, vertex triapex, REAL area)
1447 int triunsuitable(triorg, tridest, triapex, area)
1455 REAL dxoa, dxda, dxod;
1456 REAL dyoa, dyda, dyod;
1457 REAL oalen, dalen, odlen;
1460 dxoa = triorg[0] - triapex[0];
1461 dyoa = triorg[1] - triapex[1];
1462 dxda = tridest[0] - triapex[0];
1463 dyda = tridest[1] - triapex[1];
1464 dxod = triorg[0] - tridest[0];
1465 dyod = triorg[1] - tridest[1];
1467 oalen = dxoa * dxoa + dyoa * dyoa;
1468 dalen = dxda * dxda + dyda * dyda;
1469 odlen = dxod * dxod + dyod * dyod;
1471 maxlen = (dalen > oalen) ? dalen : oalen;
1472 maxlen = (odlen > maxlen) ? odlen : maxlen;
1474 if (maxlen > 0.05 * (triorg[0] * triorg[0] + triorg[1] * triorg[1]) + 0.02) {
1491 #ifdef ANSI_DECLARATORS
1492 void triexit(
int status)
1494 void triexit(status)
1502 #ifdef ANSI_DECLARATORS
1503 VOID *trimalloc(
int size)
1505 VOID *trimalloc(size)
1512 memptr = (VOID *) malloc((
unsigned int) size);
1513 if (memptr == (VOID *) NULL) {
1514 printf(
"Error: Out of memory.\n");
1520 #ifdef ANSI_DECLARATORS
1521 void trifree(VOID *memptr)
1523 void trifree(memptr)
1551 printf(
"triangle [-pAcjevngBPNEIOXzo_lQVh] input_file\n");
1553 printf(
"triangle [-pAcjevngBPNEIOXzo_iFlCQVh] input_file\n");
1557 printf(
"triangle [-prq__a__uAcDjevngBPNEIOXzo_YS__lQVh] input_file\n");
1559 printf(
"triangle [-prq__a__uAcDjevngBPNEIOXzo_YS__iFlsCQVh] input_file\n");
1563 printf(
" -p Triangulates a Planar Straight Line Graph (.poly file).\n");
1565 printf(
" -r Refines a previously generated mesh.\n");
1567 " -q Quality mesh generation. A minimum angle may be specified.\n");
1568 printf(
" -a Applies a maximum triangle area constraint.\n");
1569 printf(
" -u Applies a user-defined triangle constraint.\n");
1572 " -A Applies attributes to identify triangles in certain regions.\n");
1573 printf(
" -c Encloses the convex hull with segments.\n");
1575 printf(
" -D Conforming Delaunay: all triangles are truly Delaunay.\n");
1581 printf(
" -j Jettison unused vertices from output .node file.\n");
1582 printf(
" -e Generates an edge list.\n");
1583 printf(
" -v Generates a Voronoi diagram.\n");
1584 printf(
" -n Generates a list of triangle neighbors.\n");
1585 printf(
" -g Generates an .off file for Geomview.\n");
1586 printf(
" -B Suppresses output of boundary information.\n");
1587 printf(
" -P Suppresses output of .poly file.\n");
1588 printf(
" -N Suppresses output of .node file.\n");
1589 printf(
" -E Suppresses output of .ele file.\n");
1590 printf(
" -I Suppresses mesh iteration numbers.\n");
1591 printf(
" -O Ignores holes in .poly file.\n");
1592 printf(
" -X Suppresses use of exact arithmetic.\n");
1593 printf(
" -z Numbers all items starting from zero (rather than one).\n");
1594 printf(
" -o2 Generates second-order subparametric elements.\n");
1596 printf(
" -Y Suppresses boundary segment splitting.\n");
1597 printf(
" -S Specifies maximum number of added Steiner points.\n");
1600 printf(
" -i Uses incremental method, rather than divide-and-conquer.\n");
1601 printf(
" -F Uses Fortune's sweepline algorithm, rather than d-and-c.\n");
1603 printf(
" -l Uses vertical cuts only, rather than alternating cuts.\n");
1607 " -s Force segments into mesh by splitting (instead of using CDT).\n");
1609 printf(
" -C Check consistency of final mesh.\n");
1611 printf(
" -Q Quiet: No terminal output except errors.\n");
1612 printf(
" -V Verbose: Detailed information on what I'm doing.\n");
1613 printf(
" -h Help: Detailed instructions for Triangle.\n");
1629 printf(
"Triangle\n");
1631 "A Two-Dimensional Quality Mesh Generator and Delaunay Triangulator.\n");
1632 printf(
"Version 1.6\n\n");
1634 "Copyright 1993, 1995, 1997, 1998, 2002, 2005 Jonathan Richard Shewchuk\n");
1635 printf(
"2360 Woolsey #H / Berkeley, California 94705-1927\n");
1636 printf(
"Bugs/comments to jrs@cs.berkeley.edu\n");
1638 "Created as part of the Quake project (tools for earthquake simulation).\n");
1640 "Supported in part by NSF Grant CMS-9318163 and an NSERC 1967 Scholarship.\n");
1641 printf(
"There is no warranty whatsoever. Use at your own risk.\n");
1643 printf(
"This executable is compiled for single precision arithmetic.\n\n\n");
1645 printf(
"This executable is compiled for double precision arithmetic.\n\n\n");
1648 "Triangle generates exact Delaunay triangulations, constrained Delaunay\n");
1650 "triangulations, conforming Delaunay triangulations, Voronoi diagrams, and\n");
1652 "high-quality triangular meshes. The latter can be generated with no small\n"
1655 "or large angles, and are thus suitable for finite element analysis. If no\n"
1658 "command line switch is specified, your .node input file is read, and the\n");
1660 "Delaunay triangulation is returned in .node and .ele output files. The\n");
1661 printf(
"command syntax is:\n\n");
1662 printf(
"triangle [-prq__a__uAcDjevngBPNEIOXzo_YS__iFlsCQVh] input_file\n\n");
1664 "Underscores indicate that numbers may optionally follow certain switches.\n");
1666 "Do not leave any space between a switch and its numeric parameter.\n");
1668 "input_file must be a file with extension .node, or extension .poly if the\n");
1670 "-p switch is used. If -r is used, you must supply .node and .ele files,\n");
1672 "and possibly a .poly file and an .area file as well. The formats of these\n"
1674 printf(
"files are described below.\n\n");
1675 printf(
"Command Line Switches:\n\n");
1677 " -p Reads a Planar Straight Line Graph (.poly file), which can specify\n"
1680 " vertices, segments, holes, regional attributes, and regional area\n");
1682 " constraints. Generates a constrained Delaunay triangulation (CDT)\n"
1685 " fitting the input; or, if -s, -q, -a, or -u is used, a conforming\n");
1687 " constrained Delaunay triangulation (CCDT). If you want a truly\n");
1689 " Delaunay (not just constrained Delaunay) triangulation, use -D as\n");
1691 " well. When -p is not used, Triangle reads a .node file by default.\n"
1694 " -r Refines a previously generated mesh. The mesh is read from a .node\n"
1697 " file and an .ele file. If -p is also used, a .poly file is read\n");
1699 " and used to constrain segments in the mesh. If -a is also used\n");
1701 " (with no number following), an .area file is read and used to\n");
1703 " impose area constraints on the mesh. Further details on refinement\n"
1705 printf(
" appear below.\n");
1707 " -q Quality mesh generation by Delaunay refinement (a hybrid of Paul\n");
1709 " Chew's and Jim Ruppert's algorithms). Adds vertices to the mesh to\n"
1712 " ensure that all angles are between 20 and 140 degrees. An\n");
1714 " alternative bound on the minimum angle, replacing 20 degrees, may\n");
1716 " be specified after the `q'. The specified angle may include a\n");
1718 " decimal point, but not exponential notation. Note that a bound of\n"
1721 " theta degrees on the smallest angle also implies a bound of\n");
1723 " (180 - 2 theta) on the largest angle. If the minimum angle is 28.6\n"
1726 " degrees or smaller, Triangle is mathematically guaranteed to\n");
1728 " terminate (assuming infinite precision arithmetic--Triangle may\n");
1730 " fail to terminate if you run out of precision). In practice,\n");
1732 " Triangle often succeeds for minimum angles up to 34 degrees. For\n");
1734 " some meshes, however, you might need to reduce the minimum angle to\n"
1737 " avoid problems associated with insufficient floating-point\n");
1738 printf(
" precision.\n");
1740 " -a Imposes a maximum triangle area. If a number follows the `a', no\n");
1742 " triangle is generated whose area is larger than that number. If no\n"
1745 " number is specified, an .area file (if -r is used) or .poly file\n");
1747 " (if -r is not used) specifies a set of maximum area constraints.\n");
1749 " An .area file contains a separate area constraint for each\n");
1751 " triangle, and is useful for refining a finite element mesh based on\n"
1754 " a posteriori error estimates. A .poly file can optionally contain\n"
1757 " an area constraint for each segment-bounded region, thereby\n");
1759 " controlling triangle densities in a first triangulation of a PSLG.\n"
1762 " You can impose both a fixed area constraint and a varying area\n");
1764 " constraint by invoking the -a switch twice, once with and once\n");
1766 " without a number following. Each area specified may include a\n");
1767 printf(
" decimal point.\n");
1769 " -u Imposes a user-defined constraint on triangle size. There are two\n"
1772 " ways to use this feature. One is to edit the triunsuitable()\n");
1774 " procedure in triangle.c to encode any constraint you like, then\n");
1776 " recompile Triangle. The other is to compile triangle.c with the\n");
1778 " EXTERNAL_TEST symbol set (compiler switch -DEXTERNAL_TEST), then\n");
1780 " link Triangle with a separate object file that implements\n");
1782 " triunsuitable(). In either case, the -u switch causes the user-\n");
1783 printf(
" defined test to be applied to every triangle.\n");
1785 " -A Assigns an additional floating-point attribute to each triangle\n");
1787 " that identifies what segment-bounded region each triangle belongs\n");
1789 " to. Attributes are assigned to regions by the .poly file. If a\n");
1791 " region is not explicitly marked by the .poly file, triangles in\n");
1793 " that region are assigned an attribute of zero. The -A switch has\n");
1795 " an effect only when the -p switch is used and the -r switch is not.\n"
1798 " -c Creates segments on the convex hull of the triangulation. If you\n");
1800 " are triangulating a vertex set, this switch causes a .poly file to\n"
1803 " be written, containing all edges of the convex hull. If you are\n");
1805 " triangulating a PSLG, this switch specifies that the whole convex\n");
1807 " hull of the PSLG should be triangulated, regardless of what\n");
1809 " segments the PSLG has. If you do not use this switch when\n");
1811 " triangulating a PSLG, Triangle assumes that you have identified the\n"
1814 " region to be triangulated by surrounding it with segments of the\n");
1816 " input PSLG. Beware: if you are not careful, this switch can cause\n"
1819 " the introduction of an extremely thin angle between a PSLG segment\n"
1822 " and a convex hull segment, which can cause overrefinement (and\n");
1824 " possibly failure if Triangle runs out of precision). If you are\n");
1826 " refining a mesh, the -c switch works differently: it causes a\n");
1828 " .poly file to be written containing the boundary edges of the mesh\n"
1830 printf(
" (useful if no .poly file was read).\n");
1832 " -D Conforming Delaunay triangulation: use this switch if you want to\n"
1835 " ensure that all the triangles in the mesh are Delaunay, and not\n");
1837 " merely constrained Delaunay; or if you want to ensure that all the\n"
1840 " Voronoi vertices lie within the triangulation. (Some finite volume\n"
1843 " methods have this requirement.) This switch invokes Ruppert's\n");
1845 " original algorithm, which splits every subsegment whose diametral\n");
1847 " circle is encroached. It usually increases the number of vertices\n"
1849 printf(
" and triangles.\n");
1851 " -j Jettisons vertices that are not part of the final triangulation\n");
1853 " from the output .node file. By default, Triangle copies all\n");
1855 " vertices in the input .node file to the output .node file, in the\n");
1857 " same order, so their indices do not change. The -j switch prevents\n"
1860 " duplicated input vertices, or vertices `eaten' by holes, from\n");
1862 " appearing in the output .node file. Thus, if two input vertices\n");
1864 " have exactly the same coordinates, only the first appears in the\n");
1866 " output. If any vertices are jettisoned, the vertex numbering in\n");
1868 " the output .node file differs from that of the input .node file.\n");
1870 " -e Outputs (to an .edge file) a list of edges of the triangulation.\n");
1872 " -v Outputs the Voronoi diagram associated with the triangulation.\n");
1874 " Does not attempt to detect degeneracies, so some Voronoi vertices\n");
1876 " may be duplicated. See the discussion of Voronoi diagrams below.\n");
1878 " -n Outputs (to a .neigh file) a list of triangles neighboring each\n");
1879 printf(
" triangle.\n");
1881 " -g Outputs the mesh to an Object File Format (.off) file, suitable for\n"
1883 printf(
" viewing with the Geometry Center's Geomview package.\n");
1885 " -B No boundary markers in the output .node, .poly, and .edge output\n");
1887 " files. See the detailed discussion of boundary markers below.\n");
1889 " -P No output .poly file. Saves disk space, but you lose the ability\n");
1891 " to maintain constraining segments on later refinements of the mesh.\n"
1893 printf(
" -N No output .node file.\n");
1894 printf(
" -E No output .ele file.\n");
1896 " -I No iteration numbers. Suppresses the output of .node and .poly\n");
1898 " files, so your input files won't be overwritten. (If your input is\n"
1901 " a .poly file only, a .node file is written.) Cannot be used with\n");
1903 " the -r switch, because that would overwrite your input .ele file.\n");
1905 " Shouldn't be used with the -q, -a, -u, or -s switch if you are\n");
1907 " using a .node file for input, because no .node file is written, so\n"
1909 printf(
" there is no record of any added Steiner points.\n");
1910 printf(
" -O No holes. Ignores the holes in the .poly file.\n");
1912 " -X No exact arithmetic. Normally, Triangle uses exact floating-point\n"
1915 " arithmetic for certain tests if it thinks the inexact tests are not\n"
1918 " accurate enough. Exact arithmetic ensures the robustness of the\n");
1920 " triangulation algorithms, despite floating-point roundoff error.\n");
1922 " Disabling exact arithmetic with the -X switch causes a small\n");
1924 " improvement in speed and creates the possibility that Triangle will\n"
1926 printf(
" fail to produce a valid mesh. Not recommended.\n");
1928 " -z Numbers all items starting from zero (rather than one). Note that\n"
1931 " this switch is normally overridden by the value used to number the\n"
1934 " first vertex of the input .node or .poly file. However, this\n");
1936 " switch is useful when calling Triangle from another program.\n");
1938 " -o2 Generates second-order subparametric elements with six nodes each.\n"
1941 " -Y No new vertices on the boundary. This switch is useful when the\n");
1943 " mesh boundary must be preserved so that it conforms to some\n");
1945 " adjacent mesh. Be forewarned that you will probably sacrifice much\n"
1948 " of the quality of the mesh; Triangle will try, but the resulting\n");
1950 " mesh may contain poorly shaped triangles. Works well if all the\n");
1952 " boundary vertices are closely spaced. Specify this switch twice\n");
1954 " (`-YY') to prevent all segment splitting, including internal\n");
1955 printf(
" boundaries.\n");
1957 " -S Specifies the maximum number of Steiner points (vertices that are\n");
1959 " not in the input, but are added to meet the constraints on minimum\n"
1962 " angle and maximum area). The default is to allow an unlimited\n");
1964 " number. If you specify this switch with no number after it,\n");
1966 " the limit is set to zero. Triangle always adds vertices at segment\n"
1969 " intersections, even if it needs to use more vertices than the limit\n"
1972 " you set. When Triangle inserts segments by splitting (-s), it\n");
1974 " always adds enough vertices to ensure that all the segments of the\n"
1976 printf(
" PLSG are recovered, ignoring the limit if necessary.\n");
1978 " -i Uses an incremental rather than a divide-and-conquer algorithm to\n");
1980 " construct a Delaunay triangulation. Try it if the divide-and-\n");
1981 printf(
" conquer algorithm fails.\n");
1983 " -F Uses Steven Fortune's sweepline algorithm to construct a Delaunay\n");
1985 " triangulation. Warning: does not use exact arithmetic for all\n");
1986 printf(
" calculations. An exact result is not guaranteed.\n");
1988 " -l Uses only vertical cuts in the divide-and-conquer algorithm. By\n");
1990 " default, Triangle alternates between vertical and horizontal cuts,\n"
1993 " which usually improve the speed except with vertex sets that are\n");
1995 " small or short and wide. This switch is primarily of theoretical\n");
1996 printf(
" interest.\n");
1998 " -s Specifies that segments should be forced into the triangulation by\n"
2001 " recursively splitting them at their midpoints, rather than by\n");
2003 " generating a constrained Delaunay triangulation. Segment splitting\n"
2006 " is true to Ruppert's original algorithm, but can create needlessly\n"
2009 " small triangles. This switch is primarily of theoretical interest.\n"
2012 " -C Check the consistency of the final mesh. Uses exact arithmetic for\n"
2015 " checking, even if the -X switch is used. Useful if you suspect\n");
2016 printf(
" Triangle is buggy.\n");
2018 " -Q Quiet: Suppresses all explanation of what Triangle is doing,\n");
2019 printf(
" unless an error occurs.\n");
2021 " -V Verbose: Gives detailed information about what Triangle is doing.\n"
2024 " Add more `V's for increasing amount of detail. `-V' is most\n");
2026 " useful; itgives information on algorithmic progress and much more\n");
2028 " detailed statistics. `-VV' gives vertex-by-vertex details, and\n");
2030 " prints so much that Triangle runs much more slowly. `-VVVV' gives\n"
2032 printf(
" information only a debugger could love.\n");
2033 printf(
" -h Help: Displays these instructions.\n");
2035 printf(
"Definitions:\n");
2038 " A Delaunay triangulation of a vertex set is a triangulation whose\n");
2040 " vertices are the vertex set, that covers the convex hull of the vertex\n");
2042 " set. A Delaunay triangulation has the property that no vertex lies\n");
2044 " inside the circumscribing circle (circle that passes through all three\n");
2045 printf(
" vertices) of any triangle in the triangulation.\n\n");
2047 " A Voronoi diagram of a vertex set is a subdivision of the plane into\n");
2049 " polygonal cells (some of which may be unbounded, meaning infinitely\n");
2051 " large), where each cell is the set of points in the plane that are closer\n"
2054 " to some input vertex than to any other input vertex. The Voronoi diagram\n"
2056 printf(
" is a geometric dual of the Delaunay triangulation.\n\n");
2058 " A Planar Straight Line Graph (PSLG) is a set of vertices and segments.\n");
2060 " Segments are simply edges, whose endpoints are all vertices in the PSLG.\n"
2063 " Segments may intersect each other only at their endpoints. The file\n");
2064 printf(
" format for PSLGs (.poly files) is described below.\n\n");
2066 " A constrained Delaunay triangulation (CDT) of a PSLG is similar to a\n");
2068 " Delaunay triangulation, but each PSLG segment is present as a single edge\n"
2071 " of the CDT. (A constrained Delaunay triangulation is not truly a\n");
2073 " Delaunay triangulation, because some of its triangles might not be\n");
2075 " Delaunay.) By definition, a CDT does not have any vertices other than\n");
2077 " those specified in the input PSLG. Depending on context, a CDT might\n");
2079 " cover the convex hull of the PSLG, or it might cover only a segment-\n");
2080 printf(
" bounded region (e.g. a polygon).\n\n");
2082 " A conforming Delaunay triangulation of a PSLG is a triangulation in which\n"
2085 " each triangle is truly Delaunay, and each PSLG segment is represented by\n"
2088 " a linear contiguous sequence of edges of the triangulation. New vertices\n"
2091 " (not part of the PSLG) may appear, and each input segment may have been\n");
2093 " subdivided into shorter edges (subsegments) by these additional vertices.\n"
2096 " The new vertices are frequently necessary to maintain the Delaunay\n");
2097 printf(
" property while ensuring that every segment is represented.\n\n");
2099 " A conforming constrained Delaunay triangulation (CCDT) of a PSLG is a\n");
2101 " triangulation of a PSLG whose triangles are constrained Delaunay. New\n");
2102 printf(
" vertices may appear, and input segments may be subdivided into\n");
2104 " subsegments, but not to guarantee that segments are respected; rather, to\n"
2107 " improve the quality of the triangles. The high-quality meshes produced\n");
2109 " by the -q switch are usually CCDTs, but can be made conforming Delaunay\n");
2110 printf(
" with the -D switch.\n\n");
2111 printf(
"File Formats:\n\n");
2113 " All files may contain comments prefixed by the character '#'. Vertices,\n"
2116 " triangles, edges, holes, and maximum area constraints must be numbered\n");
2118 " consecutively, starting from either 1 or 0. Whichever you choose, all\n");
2120 " input files must be consistent; if the vertices are numbered from 1, so\n");
2122 " must be all other objects. Triangle automatically detects your choice\n");
2124 " while reading the .node (or .poly) file. (When calling Triangle from\n");
2126 " another program, use the -z switch if you wish to number objects from\n");
2127 printf(
" zero.) Examples of these file formats are given below.\n\n");
2128 printf(
" .node files:\n");
2130 " First line: <# of vertices> <dimension (must be 2)> <# of attributes>\n"
2133 " <# of boundary markers (0 or 1)>\n"
2136 " Remaining lines: <vertex #> <x> <y> [attributes] [boundary marker]\n");
2139 " The attributes, which are typically floating-point values of physical\n");
2141 " quantities (such as mass or conductivity) associated with the nodes of\n"
2144 " a finite element mesh, are copied unchanged to the output mesh. If -q,\n"
2147 " -a, -u, -D, or -s is selected, each new Steiner point added to the mesh\n"
2149 printf(
" has attributes assigned to it by linear interpolation.\n\n");
2151 " If the fourth entry of the first line is `1', the last column of the\n");
2153 " remainder of the file is assumed to contain boundary markers. Boundary\n"
2156 " markers are used to identify boundary vertices and vertices resting on\n"
2159 " PSLG segments; a complete description appears in a section below. The\n"
2162 " .node file produced by Triangle contains boundary markers in the last\n");
2163 printf(
" column unless they are suppressed by the -B switch.\n\n");
2164 printf(
" .ele files:\n");
2166 " First line: <# of triangles> <nodes per triangle> <# of attributes>\n");
2168 " Remaining lines: <triangle #> <node> <node> <node> ... [attributes]\n");
2171 " Nodes are indices into the corresponding .node file. The first three\n");
2173 " nodes are the corner vertices, and are listed in counterclockwise order\n"
2176 " around each triangle. (The remaining nodes, if any, depend on the type\n"
2178 printf(
" of finite element used.)\n\n");
2180 " The attributes are just like those of .node files. Because there is no\n"
2183 " simple mapping from input to output triangles, Triangle attempts to\n");
2185 " interpolate attributes, and may cause a lot of diffusion of attributes\n"
2188 " among nearby triangles as the triangulation is refined. Attributes do\n"
2190 printf(
" not diffuse across segments, so attributes used to identify\n");
2191 printf(
" segment-bounded regions remain intact.\n\n");
2193 " In .ele files produced by Triangle, each triangular element has three\n");
2195 " nodes (vertices) unless the -o2 switch is used, in which case\n");
2197 " subparametric quadratic elements with six nodes each are generated.\n");
2199 " The first three nodes are the corners in counterclockwise order, and\n");
2201 " the fourth, fifth, and sixth nodes lie on the midpoints of the edges\n");
2203 " opposite the first, second, and third vertices, respectively.\n");
2205 printf(
" .poly files:\n");
2207 " First line: <# of vertices> <dimension (must be 2)> <# of attributes>\n"
2210 " <# of boundary markers (0 or 1)>\n"
2213 " Following lines: <vertex #> <x> <y> [attributes] [boundary marker]\n");
2214 printf(
" One line: <# of segments> <# of boundary markers (0 or 1)>\n");
2216 " Following lines: <segment #> <endpoint> <endpoint> [boundary marker]\n");
2217 printf(
" One line: <# of holes>\n");
2218 printf(
" Following lines: <hole #> <x> <y>\n");
2220 " Optional line: <# of regional attributes and/or area constraints>\n");
2222 " Optional following lines: <region #> <x> <y> <attribute> <max area>\n");
2225 " A .poly file represents a PSLG, as well as some additional information.\n"
2228 " The first section lists all the vertices, and is identical to the\n");
2230 " format of .node files. <# of vertices> may be set to zero to indicate\n"
2233 " that the vertices are listed in a separate .node file; .poly files\n");
2235 " produced by Triangle always have this format. A vertex set represented\n"
2238 " this way has the advantage that it may easily be triangulated with or\n");
2240 " without segments (depending on whether the -p switch is invoked).\n");
2243 " The second section lists the segments. Segments are edges whose\n");
2245 " presence in the triangulation is enforced. (Depending on the choice of\n"
2248 " switches, segment might be subdivided into smaller edges). Each\n");
2250 " segment is specified by listing the indices of its two endpoints. This\n"
2253 " means that you must include its endpoints in the vertex list. Each\n");
2254 printf(
" segment, like each point, may have a boundary marker.\n\n");
2256 " If -q, -a, -u, and -s are not selected, Triangle produces a constrained\n"
2259 " Delaunay triangulation (CDT), in which each segment appears as a single\n"
2262 " edge in the triangulation. If -q, -a, -u, or -s is selected, Triangle\n"
2265 " produces a conforming constrained Delaunay triangulation (CCDT), in\n");
2267 " which segments may be subdivided into smaller edges. If -D is\n");
2269 " selected, Triangle produces a conforming Delaunay triangulation, so\n");
2271 " that every triangle is Delaunay, and not just constrained Delaunay.\n");
2274 " The third section lists holes (and concavities, if -c is selected) in\n");
2276 " the triangulation. Holes are specified by identifying a point inside\n");
2278 " each hole. After the triangulation is formed, Triangle creates holes\n");
2280 " by eating triangles, spreading out from each hole point until its\n");
2282 " progress is blocked by segments in the PSLG. You must be careful to\n");
2284 " enclose each hole in segments, or your whole triangulation might be\n");
2286 " eaten away. If the two triangles abutting a segment are eaten, the\n");
2288 " segment itself is also eaten. Do not place a hole directly on a\n");
2289 printf(
" segment; if you do, Triangle chooses one side of the segment\n");
2290 printf(
" arbitrarily.\n\n");
2292 " The optional fourth section lists regional attributes (to be assigned\n");
2294 " to all triangles in a region) and regional constraints on the maximum\n");
2296 " triangle area. Triangle reads this section only if the -A switch is\n");
2298 " used or the -a switch is used without a number following it, and the -r\n"
2301 " switch is not used. Regional attributes and area constraints are\n");
2303 " propagated in the same manner as holes: you specify a point for each\n");
2305 " attribute and/or constraint, and the attribute and/or constraint\n");
2307 " affects the whole region (bounded by segments) containing the point.\n");
2309 " If two values are written on a line after the x and y coordinate, the\n");
2311 " first such value is assumed to be a regional attribute (but is only\n");
2313 " applied if the -A switch is selected), and the second value is assumed\n"
2316 " to be a regional area constraint (but is only applied if the -a switch\n"
2319 " is selected). You may specify just one value after the coordinates,\n");
2321 " which can serve as both an attribute and an area constraint, depending\n"
2324 " on the choice of switches. If you are using the -A and -a switches\n");
2326 " simultaneously and wish to assign an attribute to some region without\n");
2327 printf(
" imposing an area constraint, use a negative maximum area.\n\n");
2329 " When a triangulation is created from a .poly file, you must either\n");
2331 " enclose the entire region to be triangulated in PSLG segments, or\n");
2333 " use the -c switch, which automatically creates extra segments that\n");
2335 " enclose the convex hull of the PSLG. If you do not use the -c switch,\n"
2338 " Triangle eats all triangles that are not enclosed by segments; if you\n");
2340 " are not careful, your whole triangulation may be eaten away. If you do\n"
2343 " use the -c switch, you can still produce concavities by the appropriate\n"
2346 " placement of holes just inside the boundary of the convex hull.\n");
2349 " An ideal PSLG has no intersecting segments, nor any vertices that lie\n");
2351 " upon segments (except, of course, the endpoints of each segment). You\n"
2354 " aren't required to make your .poly files ideal, but you should be aware\n"
2357 " of what can go wrong. Segment intersections are relatively safe--\n");
2359 " Triangle calculates the intersection points for you and adds them to\n");
2361 " the triangulation--as long as your machine's floating-point precision\n");
2363 " doesn't become a problem. You are tempting the fates if you have three\n"
2366 " segments that cross at the same location, and expect Triangle to figure\n"
2369 " out where the intersection point is. Thanks to floating-point roundoff\n"
2372 " error, Triangle will probably decide that the three segments intersect\n"
2375 " at three different points, and you will find a minuscule triangle in\n");
2377 " your output--unless Triangle tries to refine the tiny triangle, uses\n");
2379 " up the last bit of machine precision, and fails to terminate at all.\n");
2381 " You're better off putting the intersection point in the input files,\n");
2383 " and manually breaking up each segment into two. Similarly, if you\n");
2385 " place a vertex at the middle of a segment, and hope that Triangle will\n"
2388 " break up the segment at that vertex, you might get lucky. On the other\n"
2391 " hand, Triangle might decide that the vertex doesn't lie precisely on\n");
2393 " the segment, and you'll have a needle-sharp triangle in your output--or\n"
2395 printf(
" a lot of tiny triangles if you're generating a quality mesh.\n");
2398 " When Triangle reads a .poly file, it also writes a .poly file, which\n");
2400 " includes all the subsegments--the edges that are parts of input\n");
2402 " segments. If the -c switch is used, the output .poly file also\n");
2404 " includes all of the edges on the convex hull. Hence, the output .poly\n"
2407 " file is useful for finding edges associated with input segments and for\n"
2410 " setting boundary conditions in finite element simulations. Moreover,\n");
2412 " you will need the output .poly file if you plan to refine the output\n");
2414 " mesh, and don't want segments to be missing in later triangulations.\n");
2416 printf(
" .area files:\n");
2417 printf(
" First line: <# of triangles>\n");
2418 printf(
" Following lines: <triangle #> <maximum area>\n");
2421 " An .area file associates with each triangle a maximum area that is used\n"
2424 " for mesh refinement. As with other file formats, every triangle must\n");
2426 " be represented, and the triangles must be numbered consecutively. A\n");
2428 " triangle may be left unconstrained by assigning it a negative maximum\n");
2429 printf(
" area.\n\n");
2430 printf(
" .edge files:\n");
2431 printf(
" First line: <# of edges> <# of boundary markers (0 or 1)>\n");
2433 " Following lines: <edge #> <endpoint> <endpoint> [boundary marker]\n");
2436 " Endpoints are indices into the corresponding .node file. Triangle can\n"
2439 " produce .edge files (use the -e switch), but cannot read them. The\n");
2441 " optional column of boundary markers is suppressed by the -B switch.\n");
2444 " In Voronoi diagrams, one also finds a special kind of edge that is an\n");
2446 " infinite ray with only one endpoint. For these edges, a different\n");
2447 printf(
" format is used:\n\n");
2448 printf(
" <edge #> <endpoint> -1 <direction x> <direction y>\n\n");
2450 " The `direction' is a floating-point vector that indicates the direction\n"
2452 printf(
" of the infinite ray.\n\n");
2453 printf(
" .neigh files:\n");
2455 " First line: <# of triangles> <# of neighbors per triangle (always 3)>\n"
2458 " Following lines: <triangle #> <neighbor> <neighbor> <neighbor>\n");
2461 " Neighbors are indices into the corresponding .ele file. An index of -1\n"
2464 " indicates no neighbor (because the triangle is on an exterior\n");
2466 " boundary). The first neighbor of triangle i is opposite the first\n");
2467 printf(
" corner of triangle i, and so on.\n\n");
2469 " Triangle can produce .neigh files (use the -n switch), but cannot read\n"
2471 printf(
" them.\n\n");
2472 printf(
"Boundary Markers:\n\n");
2474 " Boundary markers are tags used mainly to identify which output vertices\n");
2476 " and edges are associated with which PSLG segment, and to identify which\n");
2478 " vertices and edges occur on a boundary of the triangulation. A common\n");
2480 " use is to determine where boundary conditions should be applied to a\n");
2482 " finite element mesh. You can prevent boundary markers from being written\n"
2484 printf(
" into files produced by Triangle by using the -B switch.\n\n");
2486 " The boundary marker associated with each segment in an output .poly file\n"
2488 printf(
" and each edge in an output .edge file is chosen as follows:\n");
2490 " - If an output edge is part or all of a PSLG segment with a nonzero\n");
2492 " boundary marker, then the edge is assigned the same marker.\n");
2494 " - Otherwise, if the edge lies on a boundary of the triangulation\n");
2496 " (even the boundary of a hole), then the edge is assigned the marker\n");
2497 printf(
" one (1).\n");
2498 printf(
" - Otherwise, the edge is assigned the marker zero (0).\n");
2500 " The boundary marker associated with each vertex in an output .node file\n");
2501 printf(
" is chosen as follows:\n");
2503 " - If a vertex is assigned a nonzero boundary marker in the input file,\n"
2506 " then it is assigned the same marker in the output .node file.\n");
2508 " - Otherwise, if the vertex lies on a PSLG segment (even if it is an\n");
2510 " endpoint of the segment) with a nonzero boundary marker, then the\n");
2512 " vertex is assigned the same marker. If the vertex lies on several\n");
2513 printf(
" such segments, one of the markers is chosen arbitrarily.\n");
2515 " - Otherwise, if the vertex occurs on a boundary of the triangulation,\n");
2516 printf(
" then the vertex is assigned the marker one (1).\n");
2517 printf(
" - Otherwise, the vertex is assigned the marker zero (0).\n");
2520 " If you want Triangle to determine for you which vertices and edges are on\n"
2523 " the boundary, assign them the boundary marker zero (or use no markers at\n"
2526 " all) in your input files. In the output files, all boundary vertices,\n");
2527 printf(
" edges, and segments will be assigned the value one.\n\n");
2528 printf(
"Triangulation Iteration Numbers:\n\n");
2530 " Because Triangle can read and refine its own triangulations, input\n");
2532 " and output files have iteration numbers. For instance, Triangle might\n");
2534 " read the files mesh.3.node, mesh.3.ele, and mesh.3.poly, refine the\n");
2536 " triangulation, and output the files mesh.4.node, mesh.4.ele, and\n");
2537 printf(
" mesh.4.poly. Files with no iteration number are treated as if\n");
2539 " their iteration number is zero; hence, Triangle might read the file\n");
2541 " points.node, triangulate it, and produce the files points.1.node and\n");
2542 printf(
" points.1.ele.\n\n");
2544 " Iteration numbers allow you to create a sequence of successively finer\n");
2546 " meshes suitable for multigrid methods. They also allow you to produce a\n"
2549 " sequence of meshes using error estimate-driven mesh refinement.\n");
2552 " If you're not using refinement or quality meshing, and you don't like\n");
2554 " iteration numbers, use the -I switch to disable them. This switch also\n");
2556 " disables output of .node and .poly files to prevent your input files from\n"
2559 " being overwritten. (If the input is a .poly file that contains its own\n");
2561 " points, a .node file is written. This can be quite convenient for\n");
2562 printf(
" computing CDTs or quality meshes.)\n\n");
2563 printf(
"Examples of How to Use Triangle:\n\n");
2565 " `triangle dots' reads vertices from dots.node, and writes their Delaunay\n"
2568 " triangulation to dots.1.node and dots.1.ele. (dots.1.node is identical\n");
2570 " to dots.node.) `triangle -I dots' writes the triangulation to dots.ele\n");
2572 " instead. (No additional .node file is needed, so none is written.)\n");
2575 " `triangle -pe object.1' reads a PSLG from object.1.poly (and possibly\n");
2577 " object.1.node, if the vertices are omitted from object.1.poly) and writes\n"
2580 " its constrained Delaunay triangulation to object.2.node and object.2.ele.\n"
2583 " The segments are copied to object.2.poly, and all edges are written to\n");
2584 printf(
" object.2.edge.\n\n");
2586 " `triangle -pq31.5a.1 object' reads a PSLG from object.poly (and possibly\n"
2589 " object.node), generates a mesh whose angles are all between 31.5 and 117\n"
2592 " degrees and whose triangles all have areas of 0.1 or less, and writes the\n"
2595 " mesh to object.1.node and object.1.ele. Each segment may be broken up\n");
2596 printf(
" into multiple subsegments; these are written to object.1.poly.\n");
2599 " Here is a sample file `box.poly' describing a square with a square hole:\n"
2603 " # A box with eight vertices in 2D, no attributes, one boundary marker.\n"
2605 printf(
" 8 2 0 1\n");
2606 printf(
" # Outer box has these vertices:\n");
2607 printf(
" 1 0 0 0\n");
2608 printf(
" 2 0 3 0\n");
2609 printf(
" 3 3 0 0\n");
2610 printf(
" 4 3 3 33 # A special marker for this vertex.\n");
2611 printf(
" # Inner square has these vertices:\n");
2612 printf(
" 5 1 1 0\n");
2613 printf(
" 6 1 2 0\n");
2614 printf(
" 7 2 1 0\n");
2615 printf(
" 8 2 2 0\n");
2616 printf(
" # Five segments with boundary markers.\n");
2618 printf(
" 1 1 2 5 # Left side of outer box.\n");
2619 printf(
" # Square hole has these segments:\n");
2620 printf(
" 2 5 7 0\n");
2621 printf(
" 3 7 8 0\n");
2622 printf(
" 4 8 6 10\n");
2623 printf(
" 5 6 5 0\n");
2624 printf(
" # One hole in the middle of the inner square.\n");
2626 printf(
" 1 1.5 1.5\n");
2629 " Note that some segments are missing from the outer square, so you must\n");
2631 " use the `-c' switch. After `triangle -pqc box.poly', here is the output\n"
2634 " file `box.1.node', with twelve vertices. The last four vertices were\n");
2636 " added to meet the angle constraint. Vertices 1, 2, and 9 have markers\n");
2638 " from segment 1. Vertices 6 and 8 have markers from segment 4. All the\n");
2640 " other vertices but 4 have been marked to indicate that they lie on a\n");
2641 printf(
" boundary.\n\n");
2642 printf(
" 12 2 0 1\n");
2643 printf(
" 1 0 0 5\n");
2644 printf(
" 2 0 3 5\n");
2645 printf(
" 3 3 0 1\n");
2646 printf(
" 4 3 3 33\n");
2647 printf(
" 5 1 1 1\n");
2648 printf(
" 6 1 2 10\n");
2649 printf(
" 7 2 1 1\n");
2650 printf(
" 8 2 2 10\n");
2651 printf(
" 9 0 1.5 5\n");
2652 printf(
" 10 1.5 0 1\n");
2653 printf(
" 11 3 1.5 1\n");
2654 printf(
" 12 1.5 3 1\n");
2655 printf(
" # Generated by triangle -pqc box.poly\n");
2657 printf(
" Here is the output file `box.1.ele', with twelve triangles.\n");
2659 printf(
" 12 3 0\n");
2660 printf(
" 1 5 6 9\n");
2661 printf(
" 2 10 3 7\n");
2662 printf(
" 3 6 8 12\n");
2663 printf(
" 4 9 1 5\n");
2664 printf(
" 5 6 2 9\n");
2665 printf(
" 6 7 3 11\n");
2666 printf(
" 7 11 4 8\n");
2667 printf(
" 8 7 5 10\n");
2668 printf(
" 9 12 2 6\n");
2669 printf(
" 10 8 7 11\n");
2670 printf(
" 11 5 1 10\n");
2671 printf(
" 12 8 4 12\n");
2672 printf(
" # Generated by triangle -pqc box.poly\n\n");
2674 " Here is the output file `box.1.poly'. Note that segments have been added\n"
2677 " to represent the convex hull, and some segments have been subdivided by\n");
2679 " newly added vertices. Note also that <# of vertices> is set to zero to\n");
2680 printf(
" indicate that the vertices should be read from the .node file.\n");
2682 printf(
" 0 2 0 1\n");
2684 printf(
" 1 1 9 5\n");
2685 printf(
" 2 5 7 1\n");
2686 printf(
" 3 8 7 1\n");
2687 printf(
" 4 6 8 10\n");
2688 printf(
" 5 5 6 1\n");
2689 printf(
" 6 3 10 1\n");
2690 printf(
" 7 4 11 1\n");
2691 printf(
" 8 2 12 1\n");
2692 printf(
" 9 9 2 5\n");
2693 printf(
" 10 10 1 1\n");
2694 printf(
" 11 11 3 1\n");
2695 printf(
" 12 12 4 1\n");
2697 printf(
" 1 1.5 1.5\n");
2698 printf(
" # Generated by triangle -pqc box.poly\n");
2700 printf(
"Refinement and Area Constraints:\n");
2703 " The -r switch causes a mesh (.node and .ele files) to be read and\n");
2705 " refined. If the -p switch is also used, a .poly file is read and used to\n"
2708 " specify edges that are constrained and cannot be eliminated (although\n");
2710 " they can be subdivided into smaller edges) by the refinement process.\n");
2713 " When you refine a mesh, you generally want to impose tighter constraints.\n"
2716 " One way to accomplish this is to use -q with a larger angle, or -a\n");
2718 " followed by a smaller area than you used to generate the mesh you are\n");
2720 " refining. Another way to do this is to create an .area file, which\n");
2722 " specifies a maximum area for each triangle, and use the -a switch\n");
2724 " (without a number following). Each triangle's area constraint is applied\n"
2727 " to that triangle. Area constraints tend to diffuse as the mesh is\n");
2729 " refined, so if there are large variations in area constraint between\n");
2731 " adjacent triangles, you may not get the results you want. In that case,\n"
2734 " consider instead using the -u switch and writing a C procedure that\n");
2735 printf(
" determines which triangles are too large.\n\n");
2737 " If you are refining a mesh composed of linear (three-node) elements, the\n"
2740 " output mesh contains all the nodes present in the input mesh, in the same\n"
2743 " order, with new nodes added at the end of the .node file. However, the\n");
2745 " refinement is not hierarchical: there is no guarantee that each output\n");
2747 " element is contained in a single input element. Often, an output element\n"
2750 " can overlap two or three input elements, and some input edges are not\n");
2752 " present in the output mesh. Hence, a sequence of refined meshes forms a\n"
2755 " hierarchy of nodes, but not a hierarchy of elements. If you refine a\n");
2757 " mesh of higher-order elements, the hierarchical property applies only to\n"
2760 " the nodes at the corners of an element; the midpoint nodes on each edge\n");
2761 printf(
" are discarded before the mesh is refined.\n\n");
2763 " Maximum area constraints in .poly files operate differently from those in\n"
2766 " .area files. A maximum area in a .poly file applies to the whole\n");
2768 " (segment-bounded) region in which a point falls, whereas a maximum area\n");
2770 " in an .area file applies to only one triangle. Area constraints in .poly\n"
2773 " files are used only when a mesh is first generated, whereas area\n");
2775 " constraints in .area files are used only to refine an existing mesh, and\n"
2778 " are typically based on a posteriori error estimates resulting from a\n");
2779 printf(
" finite element simulation on that mesh.\n\n");
2781 " `triangle -rq25 object.1' reads object.1.node and object.1.ele, then\n");
2783 " refines the triangulation to enforce a 25 degree minimum angle, and then\n"
2786 " writes the refined triangulation to object.2.node and object.2.ele.\n");
2789 " `triangle -rpaa6.2 z.3' reads z.3.node, z.3.ele, z.3.poly, and z.3.area.\n"
2792 " After reconstructing the mesh and its subsegments, Triangle refines the\n");
2794 " mesh so that no triangle has area greater than 6.2, and furthermore the\n");
2796 " triangles satisfy the maximum area constraints in z.3.area. No angle\n");
2798 " bound is imposed at all. The output is written to z.4.node, z.4.ele, and\n"
2800 printf(
" z.4.poly.\n\n");
2802 " The sequence `triangle -qa1 x', `triangle -rqa.3 x.1', `triangle -rqa.1\n");
2804 " x.2' creates a sequence of successively finer meshes x.1, x.2, and x.3,\n");
2805 printf(
" suitable for multigrid.\n\n");
2806 printf(
"Convex Hulls and Mesh Boundaries:\n\n");
2808 " If the input is a vertex set (not a PSLG), Triangle produces its convex\n");
2810 " hull as a by-product in the output .poly file if you use the -c switch.\n");
2812 " There are faster algorithms for finding a two-dimensional convex hull\n");
2813 printf(
" than triangulation, of course, but this one comes for free.\n\n");
2815 " If the input is an unconstrained mesh (you are using the -r switch but\n");
2817 " not the -p switch), Triangle produces a list of its boundary edges\n");
2819 " (including hole boundaries) as a by-product when you use the -c switch.\n");
2821 " If you also use the -p switch, the output .poly file contains all the\n");
2822 printf(
" segments from the input .poly file as well.\n\n");
2823 printf(
"Voronoi Diagrams:\n\n");
2825 " The -v switch produces a Voronoi diagram, in files suffixed .v.node and\n");
2827 " .v.edge. For example, `triangle -v points' reads points.node, produces\n");
2829 " its Delaunay triangulation in points.1.node and points.1.ele, and\n");
2831 " produces its Voronoi diagram in points.1.v.node and points.1.v.edge. The\n"
2834 " .v.node file contains a list of all Voronoi vertices, and the .v.edge\n");
2836 " file contains a list of all Voronoi edges, some of which may be infinite\n"
2839 " rays. (The choice of filenames makes it easy to run the set of Voronoi\n");
2840 printf(
" vertices through Triangle, if so desired.)\n\n");
2842 " This implementation does not use exact arithmetic to compute the Voronoi\n"
2845 " vertices, and does not check whether neighboring vertices are identical.\n"
2848 " Be forewarned that if the Delaunay triangulation is degenerate or\n");
2850 " near-degenerate, the Voronoi diagram may have duplicate vertices or\n");
2851 printf(
" crossing edges.\n\n");
2853 " The result is a valid Voronoi diagram only if Triangle's output is a true\n"
2856 " Delaunay triangulation. The Voronoi output is usually meaningless (and\n");
2858 " may contain crossing edges and other pathology) if the output is a CDT or\n"
2861 " CCDT, or if it has holes or concavities. If the triangulated domain is\n");
2863 " convex and has no holes, you can use -D switch to force Triangle to\n");
2865 " construct a conforming Delaunay triangulation instead of a CCDT, so the\n");
2866 printf(
" Voronoi diagram will be valid.\n\n");
2867 printf(
"Mesh Topology:\n\n");
2869 " You may wish to know which triangles are adjacent to a certain Delaunay\n");
2871 " edge in an .edge file, which Voronoi cells are adjacent to a certain\n");
2873 " Voronoi edge in a .v.edge file, or which Voronoi cells are adjacent to\n");
2875 " each other. All of this information can be found by cross-referencing\n");
2877 " output files with the recollection that the Delaunay triangulation and\n");
2878 printf(
" the Voronoi diagram are planar duals.\n\n");
2880 " Specifically, edge i of an .edge file is the dual of Voronoi edge i of\n");
2882 " the corresponding .v.edge file, and is rotated 90 degrees counterclock-\n");
2884 " wise from the Voronoi edge. Triangle j of an .ele file is the dual of\n");
2886 " vertex j of the corresponding .v.node file. Voronoi cell k is the dual\n");
2887 printf(
" of vertex k of the corresponding .node file.\n\n");
2889 " Hence, to find the triangles adjacent to a Delaunay edge, look at the\n");
2891 " vertices of the corresponding Voronoi edge. If the endpoints of a\n");
2893 " Voronoi edge are Voronoi vertices 2 and 6 respectively, then triangles 2\n"
2896 " and 6 adjoin the left and right sides of the corresponding Delaunay edge,\n"
2899 " respectively. To find the Voronoi cells adjacent to a Voronoi edge, look\n"
2902 " at the endpoints of the corresponding Delaunay edge. If the endpoints of\n"
2905 " a Delaunay edge are input vertices 7 and 12, then Voronoi cells 7 and 12\n"
2908 " adjoin the right and left sides of the corresponding Voronoi edge,\n");
2910 " respectively. To find which Voronoi cells are adjacent to each other,\n");
2911 printf(
" just read the list of Delaunay edges.\n\n");
2913 " Triangle does not write a list of the edges adjoining each Voronoi cell,\n"
2916 " but you can reconstructed it straightforwardly. For instance, to find\n");
2918 " all the edges of Voronoi cell 1, search the output .edge file for every\n");
2920 " edge that has input vertex 1 as an endpoint. The corresponding dual\n");
2922 " edges in the output .v.edge file form the boundary of Voronoi cell 1.\n");
2925 " For each Voronoi vertex, the .neigh file gives a list of the three\n");
2927 " Voronoi vertices attached to it. You might find this more convenient\n");
2928 printf(
" than the .v.edge file.\n\n");
2929 printf(
"Quadratic Elements:\n\n");
2931 " Triangle generates meshes with subparametric quadratic elements if the\n");
2933 " -o2 switch is specified. Quadratic elements have six nodes per element,\n"
2936 " rather than three. `Subparametric' means that the edges of the triangles\n"
2939 " are always straight, so that subparametric quadratic elements are\n");
2941 " geometrically identical to linear elements, even though they can be used\n"
2944 " with quadratic interpolating functions. The three extra nodes of an\n");
2946 " element fall at the midpoints of the three edges, with the fourth, fifth,\n"
2949 " and sixth nodes appearing opposite the first, second, and third corners\n");
2950 printf(
" respectively.\n\n");
2951 printf(
"Domains with Small Angles:\n\n");
2953 " If two input segments adjoin each other at a small angle, clearly the -q\n"
2956 " switch cannot remove the small angle. Moreover, Triangle may have no\n");
2958 " choice but to generate additional triangles whose smallest angles are\n");
2960 " smaller than the specified bound. However, these triangles only appear\n");
2962 " between input segments separated by small angles. Moreover, if you\n");
2964 " request a minimum angle of theta degrees, Triangle will generally produce\n"
2967 " no angle larger than 180 - 2 theta, even if it is forced to compromise on\n"
2969 printf(
" the minimum angle.\n\n");
2970 printf(
"Statistics:\n\n");
2972 " After generating a mesh, Triangle prints a count of entities in the\n");
2974 " output mesh, including the number of vertices, triangles, edges, exterior\n"
2977 " boundary edges (i.e. subsegments on the boundary of the triangulation,\n");
2979 " including hole boundaries), interior boundary edges (i.e. subsegments of\n"
2982 " input segments not on the boundary), and total subsegments. If you've\n");
2984 " forgotten the statistics for an existing mesh, run Triangle on that mesh\n"
2987 " with the -rNEP switches to read the mesh and print the statistics without\n"
2990 " writing any files. Use -rpNEP if you've got a .poly file for the mesh.\n");
2993 " The -V switch produces extended statistics, including a rough estimate\n");
2995 " of memory use, the number of calls to geometric predicates, and\n");
2997 " histograms of the angles and the aspect ratios of the triangles in the\n");
2998 printf(
" mesh.\n\n");
2999 printf(
"Exact Arithmetic:\n\n");
3001 " Triangle uses adaptive exact arithmetic to perform what computational\n");
3003 " geometers call the `orientation' and `incircle' tests. If the floating-\n"
3006 " point arithmetic of your machine conforms to the IEEE 754 standard (as\n");
3008 " most workstations do), and does not use extended precision internal\n");
3010 " floating-point registers, then your output is guaranteed to be an\n");
3012 " absolutely true Delaunay or constrained Delaunay triangulation, roundoff\n"
3015 " error notwithstanding. The word `adaptive' implies that these arithmetic\n"
3018 " routines compute the result only to the precision necessary to guarantee\n"
3021 " correctness, so they are usually nearly as fast as their approximate\n");
3022 printf(
" counterparts.\n\n");
3024 " May CPUs, including Intel x86 processors, have extended precision\n");
3026 " floating-point registers. These must be reconfigured so their precision\n"
3029 " is reduced to memory precision. Triangle does this if it is compiled\n");
3030 printf(
" correctly. See the makefile for details.\n\n");
3032 " The exact tests can be disabled with the -X switch. On most inputs, this\n"
3035 " switch reduces the computation time by about eight percent--it's not\n");
3037 " worth the risk. There are rare difficult inputs (having many collinear\n");
3039 " and cocircular vertices), however, for which the difference in speed\n");
3041 " could be a factor of two. Be forewarned that these are precisely the\n");
3043 " inputs most likely to cause errors if you use the -X switch. Hence, the\n"
3045 printf(
" -X switch is not recommended.\n\n");
3047 " Unfortunately, the exact tests don't solve every numerical problem.\n");
3049 " Exact arithmetic is not used to compute the positions of new vertices,\n");
3051 " because the bit complexity of vertex coordinates would grow without\n");
3053 " bound. Hence, segment intersections aren't computed exactly; in very\n");
3055 " unusual cases, roundoff error in computing an intersection point might\n");
3057 " actually lead to an inverted triangle and an invalid triangulation.\n");
3059 " (This is one reason to specify your own intersection points in your .poly\n"
3062 " files.) Similarly, exact arithmetic is not used to compute the vertices\n"
3064 printf(
" of the Voronoi diagram.\n\n");
3066 " Another pair of problems not solved by the exact arithmetic routines is\n");
3068 " underflow and overflow. If Triangle is compiled for double precision\n");
3070 " arithmetic, I believe that Triangle's geometric predicates work correctly\n"
3073 " if the exponent of every input coordinate falls in the range [-148, 201].\n"
3076 " Underflow can silently prevent the orientation and incircle tests from\n");
3078 " being performed exactly, while overflow typically causes a floating\n");
3079 printf(
" exception.\n\n");
3080 printf(
"Calling Triangle from Another Program:\n\n");
3081 printf(
" Read the file triangle.h for details.\n\n");
3082 printf(
"Troubleshooting:\n\n");
3083 printf(
" Please read this section before mailing me bugs.\n\n");
3084 printf(
" `My output mesh has no triangles!'\n\n");
3086 " If you're using a PSLG, you've probably failed to specify a proper set\n"
3089 " of bounding segments, or forgotten to use the -c switch. Or you may\n");
3091 " have placed a hole badly, thereby eating all your triangles. To test\n");
3092 printf(
" these possibilities, try again with the -c and -O switches.\n");
3094 " Alternatively, all your input vertices may be collinear, in which case\n"
3096 printf(
" you can hardly expect to triangulate them.\n\n");
3097 printf(
" `Triangle doesn't terminate, or just crashes.'\n\n");
3099 " Bad things can happen when triangles get so small that the distance\n");
3101 " between their vertices isn't much larger than the precision of your\n");
3103 " machine's arithmetic. If you've compiled Triangle for single-precision\n"
3106 " arithmetic, you might do better by recompiling it for double-precision.\n"
3109 " Then again, you might just have to settle for more lenient constraints\n"
3112 " on the minimum angle and the maximum area than you had planned.\n");
3115 " You can minimize precision problems by ensuring that the origin lies\n");
3117 " inside your vertex set, or even inside the densest part of your\n");
3119 " mesh. If you're triangulating an object whose x-coordinates all fall\n");
3121 " between 6247133 and 6247134, you're not leaving much floating-point\n");
3122 printf(
" precision for Triangle to work with.\n\n");
3124 " Precision problems can occur covertly if the input PSLG contains two\n");
3126 " segments that meet (or intersect) at an extremely small angle, or if\n");
3128 " such an angle is introduced by the -c switch. If you don't realize\n");
3130 " that a tiny angle is being formed, you might never discover why\n");
3132 " Triangle is crashing. To check for this possibility, use the -S switch\n"
3135 " (with an appropriate limit on the number of Steiner points, found by\n");
3137 " trial-and-error) to stop Triangle early, and view the output .poly file\n"
3140 " with Show Me (described below). Look carefully for regions where dense\n"
3143 " clusters of vertices are forming and for small angles between segments.\n"
3146 " Zoom in closely, as such segments might look like a single segment from\n"
3148 printf(
" a distance.\n\n");
3150 " If some of the input values are too large, Triangle may suffer a\n");
3152 " floating exception due to overflow when attempting to perform an\n");
3154 " orientation or incircle test. (Read the section on exact arithmetic\n");
3156 " above.) Again, I recommend compiling Triangle for double (rather\n");
3157 printf(
" than single) precision arithmetic.\n\n");
3159 " Unexpected problems can arise if you use quality meshing (-q, -a, or\n");
3161 " -u) with an input that is not segment-bounded--that is, if your input\n");
3163 " is a vertex set, or you're using the -c switch. If the convex hull of\n"
3166 " your input vertices has collinear vertices on its boundary, an input\n");
3168 " vertex that you think lies on the convex hull might actually lie just\n");
3170 " inside the convex hull. If so, the vertex and the nearby convex hull\n");
3172 " edge form an extremely thin triangle. When Triangle tries to refine\n");
3174 " the mesh to enforce angle and area constraints, Triangle might generate\n"
3177 " extremely tiny triangles, or it might fail because of insufficient\n");
3178 printf(
" floating-point precision.\n\n");
3180 " `The numbering of the output vertices doesn't match the input vertices.'\n"
3184 " You may have had duplicate input vertices, or you may have eaten some\n");
3186 " of your input vertices with a hole, or by placing them outside the area\n"
3189 " enclosed by segments. In any case, you can solve the problem by not\n");
3190 printf(
" using the -j switch.\n\n");
3192 " `Triangle executes without incident, but when I look at the resulting\n");
3194 " mesh, it has overlapping triangles or other geometric inconsistencies.'\n");
3197 " If you select the -X switch, Triangle occasionally makes mistakes due\n");
3199 " to floating-point roundoff error. Although these errors are rare,\n");
3201 " don't use the -X switch. If you still have problems, please report the\n"
3203 printf(
" bug.\n\n");
3205 " `Triangle executes without incident, but when I look at the resulting\n");
3206 printf(
" Voronoi diagram, it has overlapping edges or other geometric\n");
3207 printf(
" inconsistencies.'\n");
3210 " If your input is a PSLG (-p), you can only expect a meaningful Voronoi\n"
3213 " diagram if the domain you are triangulating is convex and free of\n");
3215 " holes, and you use the -D switch to construct a conforming Delaunay\n");
3216 printf(
" triangulation (instead of a CDT or CCDT).\n\n");
3218 " Strange things can happen if you've taken liberties with your PSLG. Do\n");
3220 " you have a vertex lying in the middle of a segment? Triangle sometimes\n");
3222 " copes poorly with that sort of thing. Do you want to lay out a collinear\n"
3225 " row of evenly spaced, segment-connected vertices? Have you simply\n");
3227 " defined one long segment connecting the leftmost vertex to the rightmost\n"
3230 " vertex, and a bunch of vertices lying along it? This method occasionally\n"
3233 " works, especially with horizontal and vertical lines, but often it\n");
3235 " doesn't, and you'll have to connect each adjacent pair of vertices with a\n"
3237 printf(
" separate segment. If you don't like it, tough.\n\n");
3239 " Furthermore, if you have segments that intersect other than at their\n");
3241 " endpoints, try not to let the intersections fall extremely close to PSLG\n"
3243 printf(
" vertices or each other.\n\n");
3245 " If you have problems refining a triangulation not produced by Triangle:\n");
3247 " Are you sure the triangulation is geometrically valid? Is it formatted\n");
3249 " correctly for Triangle? Are the triangles all listed so the first three\n"
3252 " vertices are their corners in counterclockwise order? Are all of the\n");
3254 " triangles constrained Delaunay? Triangle's Delaunay refinement algorithm\n"
3256 printf(
" assumes that it starts with a CDT.\n\n");
3257 printf(
"Show Me:\n\n");
3259 " Triangle comes with a separate program named `Show Me', whose primary\n");
3261 " purpose is to draw meshes on your screen or in PostScript. Its secondary\n"
3264 " purpose is to check the validity of your input files, and do so more\n");
3266 " thoroughly than Triangle does. Unlike Triangle, Show Me requires that\n");
3268 " you have the X Windows system. Sorry, Microsoft Windows users.\n");
3270 printf(
"Triangle on the Web:\n");
3272 printf(
" To see an illustrated version of these instructions, check out\n");
3274 printf(
" http://www.cs.cmu.edu/~quake/triangle.html\n");
3276 printf(
"A Brief Plea:\n");
3279 " If you use Triangle, and especially if you use it to accomplish real\n");
3281 " work, I would like very much to hear from you. A short letter or email\n");
3283 " (to jrs@cs.berkeley.edu) describing how you use Triangle will mean a lot\n"
3286 " to me. The more people I know are using this program, the more easily I\n"
3289 " can justify spending time on improvements, which in turn will benefit\n");
3291 " you. Also, I can put you on a list to receive email whenever a new\n");
3292 printf(
" version of Triangle is available.\n\n");
3294 " If you use a mesh generated by Triangle in a publication, please include\n"
3297 " an acknowledgment as well. And please spell Triangle with a capital `T'!\n"
3300 " If you want to include a citation, use `Jonathan Richard Shewchuk,\n");
3302 " ``Triangle: Engineering a 2D Quality Mesh Generator and Delaunay\n");
3304 " Triangulator,'' in Applied Computational Geometry: Towards Geometric\n");
3306 " Engineering (Ming C. Lin and Dinesh Manocha, editors), volume 1148 of\n");
3308 " Lecture Notes in Computer Science, pages 203-222, Springer-Verlag,\n");
3310 " Berlin, May 1996. (From the First ACM Workshop on Applied Computational\n"
3312 printf(
" Geometry.)'\n\n");
3313 printf(
"Research credit:\n\n");
3315 " Of course, I can take credit for only a fraction of the ideas that made\n");
3317 " this mesh generator possible. Triangle owes its existence to the efforts\n"
3320 " of many fine computational geometers and other researchers, including\n");
3322 " Marshall Bern, L. Paul Chew, Kenneth L. Clarkson, Boris Delaunay, Rex A.\n"
3325 " Dwyer, David Eppstein, Steven Fortune, Leonidas J. Guibas, Donald E.\n");
3327 " Knuth, Charles L. Lawson, Der-Tsai Lee, Gary L. Miller, Ernst P. Mucke,\n");
3329 " Steven E. Pav, Douglas M. Priest, Jim Ruppert, Isaac Saias, Bruce J.\n");
3331 " Schachter, Micha Sharir, Peter W. Shor, Daniel D. Sleator, Jorge Stolfi,\n"
3333 printf(
" Robert E. Tarjan, Alper Ungor, Christopher J. Van Wyk, Noel J.\n");
3335 " Walkington, and Binhai Zhu. See the comments at the beginning of the\n");
3336 printf(
" source code for references.\n\n");
3348 void internalerror()
3350 printf(
" Please report this bug to jrs@cs.berkeley.edu\n");
3351 printf(
" Include the message above, your input data set, and the exact\n");
3352 printf(
" command line you used to run Triangle.\n");
3363 #ifdef ANSI_DECLARATORS
3364 void parsecommandline(
int argc,
char **argv,
struct behavior *b)
3366 void parsecommandline(argc, argv, b)
3374 #define STARTINDEX 0
3376 #define STARTINDEX 1
3381 char workstring[FILENAMESIZE];
3383 b->poly = b->refine = b->quality = 0;
3384 b->vararea = b->fixedarea = b->usertest = 0;
3385 b->regionattrib = b->convex = b->weighted = b->jettison = 0;
3387 b->edgesout = b->voronoi = b->neighbors = b->geomview = 0;
3388 b->nobound = b->nopolywritten = b->nonodewritten = b->noelewritten = 0;
3389 b->noiterationnum = 0;
3390 b->noholes = b->noexact = 0;
3391 b->incremental = b->sweepline = 0;
3401 b->quiet = b->verbose = 0;
3403 b->innodefilename[0] =
'\0';
3406 for (i = STARTINDEX; i < argc; i++) {
3408 if (argv[i][0] ==
'-') {
3410 for (j = STARTINDEX; argv[i][j] !=
'\0'; j++) {
3411 if (argv[i][j] ==
'p') {
3415 if (argv[i][j] ==
'r') {
3418 if (argv[i][j] ==
'q') {
3420 if (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
3421 (argv[i][j + 1] ==
'.')) {
3423 while (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
3424 (argv[i][j + 1] ==
'.')) {
3426 workstring[k] = argv[i][j];
3429 workstring[k] =
'\0';
3430 b->minangle = (REAL) strtod(workstring, (
char **) NULL);
3435 if (argv[i][j] ==
'a') {
3437 if (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
3438 (argv[i][j + 1] ==
'.')) {
3441 while (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
3442 (argv[i][j + 1] ==
'.')) {
3444 workstring[k] = argv[i][j];
3447 workstring[k] =
'\0';
3448 b->maxarea = (REAL) strtod(workstring, (
char **) NULL);
3449 if (b->maxarea <= 0.0) {
3450 printf(
"Error: Maximum area must be greater than zero.\n");
3457 if (argv[i][j] ==
'u') {
3462 if (argv[i][j] ==
'A') {
3463 b->regionattrib = 1;
3465 if (argv[i][j] ==
'c') {
3468 if (argv[i][j] ==
'w') {
3471 if (argv[i][j] ==
'W') {
3474 if (argv[i][j] ==
'j') {
3477 if (argv[i][j] ==
'z') {
3480 if (argv[i][j] ==
'e') {
3483 if (argv[i][j] ==
'v') {
3486 if (argv[i][j] ==
'n') {
3489 if (argv[i][j] ==
'g') {
3492 if (argv[i][j] ==
'B') {
3495 if (argv[i][j] ==
'P') {
3496 b->nopolywritten = 1;
3498 if (argv[i][j] ==
'N') {
3499 b->nonodewritten = 1;
3501 if (argv[i][j] ==
'E') {
3502 b->noelewritten = 1;
3505 if (argv[i][j] ==
'I') {
3506 b->noiterationnum = 1;
3509 if (argv[i][j] ==
'O') {
3512 if (argv[i][j] ==
'X') {
3515 if (argv[i][j] ==
'o') {
3516 if (argv[i][j + 1] ==
'2') {
3522 if (argv[i][j] ==
'Y') {
3525 if (argv[i][j] ==
'S') {
3527 while ((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) {
3529 b->steiner = b->steiner * 10 + (int) (argv[i][j] -
'0');
3534 if (argv[i][j] ==
'i') {
3537 if (argv[i][j] ==
'F') {
3541 if (argv[i][j] ==
'l') {
3546 if (argv[i][j] ==
's') {
3549 if ((argv[i][j] ==
'D') || (argv[i][j] ==
'L')) {
3554 if (argv[i][j] ==
'C') {
3558 if (argv[i][j] ==
'Q') {
3561 if (argv[i][j] ==
'V') {
3565 if ((argv[i][j] ==
'h') || (argv[i][j] ==
'H') ||
3566 (argv[i][j] ==
'?')) {
3573 strncpy(b->innodefilename, argv[i], FILENAMESIZE - 1);
3574 b->innodefilename[FILENAMESIZE - 1] =
'\0';
3579 if (b->innodefilename[0] ==
'\0') {
3582 if (!strcmp(&b->innodefilename[strlen(b->innodefilename) - 5],
".node")) {
3583 b->innodefilename[strlen(b->innodefilename) - 5] =
'\0';
3585 if (!strcmp(&b->innodefilename[strlen(b->innodefilename) - 5],
".poly")) {
3586 b->innodefilename[strlen(b->innodefilename) - 5] =
'\0';
3590 if (!strcmp(&b->innodefilename[strlen(b->innodefilename) - 4],
".ele")) {
3591 b->innodefilename[strlen(b->innodefilename) - 4] =
'\0';
3594 if (!strcmp(&b->innodefilename[strlen(b->innodefilename) - 5],
".area")) {
3595 b->innodefilename[strlen(b->innodefilename) - 5] =
'\0';
3602 b->usesegments = b->poly || b->refine || b->quality || b->convex;
3603 b->goodangle = cos(b->minangle * PI / 180.0);
3604 if (b->goodangle == 1.0) {
3605 b->offconstant = 0.0;
3607 b->offconstant = 0.475 * sqrt((1.0 + b->goodangle) / (1.0 - b->goodangle));
3609 b->goodangle *= b->goodangle;
3610 if (b->refine && b->noiterationnum) {
3612 "Error: You cannot use the -I switch when refining a triangulation.\n");
3617 if (!b->refine && !b->poly) {
3622 if (b->refine || !b->poly) {
3623 b->regionattrib = 0;
3627 if (b->weighted && (b->poly || b->quality)) {
3630 printf(
"Warning: weighted triangulations (-w, -W) are incompatible\n");
3631 printf(
" with PSLGs (-p) and meshing (-q, -a, -u). Weights ignored.\n"
3635 if (b->jettison && b->nonodewritten && !b->quiet) {
3636 printf(
"Warning: -j and -N switches are somewhat incompatible.\n");
3637 printf(
" If any vertices are jettisoned, you will need the output\n");
3638 printf(
" .node file to reconstruct the new node indices.");
3642 strcpy(b->inpolyfilename, b->innodefilename);
3643 strcpy(b->inelefilename, b->innodefilename);
3644 strcpy(b->areafilename, b->innodefilename);
3646 strcpy(workstring, b->innodefilename);
3648 while (workstring[j] !=
'\0') {
3649 if ((workstring[j] ==
'.') && (workstring[j + 1] !=
'\0')) {
3655 if (increment > 0) {
3658 if ((workstring[j] >=
'0') && (workstring[j] <=
'9')) {
3659 meshnumber = meshnumber * 10 + (int) (workstring[j] -
'0');
3664 }
while (workstring[j] !=
'\0');
3666 if (b->noiterationnum) {
3667 strcpy(b->outnodefilename, b->innodefilename);
3668 strcpy(b->outelefilename, b->innodefilename);
3669 strcpy(b->edgefilename, b->innodefilename);
3670 strcpy(b->vnodefilename, b->innodefilename);
3671 strcpy(b->vedgefilename, b->innodefilename);
3672 strcpy(b->neighborfilename, b->innodefilename);
3673 strcpy(b->offfilename, b->innodefilename);
3674 strcat(b->outnodefilename,
".node");
3675 strcat(b->outelefilename,
".ele");
3676 strcat(b->edgefilename,
".edge");
3677 strcat(b->vnodefilename,
".v.node");
3678 strcat(b->vedgefilename,
".v.edge");
3679 strcat(b->neighborfilename,
".neigh");
3680 strcat(b->offfilename,
".off");
3681 }
else if (increment == 0) {
3682 strcpy(b->outnodefilename, b->innodefilename);
3683 strcpy(b->outpolyfilename, b->innodefilename);
3684 strcpy(b->outelefilename, b->innodefilename);
3685 strcpy(b->edgefilename, b->innodefilename);
3686 strcpy(b->vnodefilename, b->innodefilename);
3687 strcpy(b->vedgefilename, b->innodefilename);
3688 strcpy(b->neighborfilename, b->innodefilename);
3689 strcpy(b->offfilename, b->innodefilename);
3690 strcat(b->outnodefilename,
".1.node");
3691 strcat(b->outpolyfilename,
".1.poly");
3692 strcat(b->outelefilename,
".1.ele");
3693 strcat(b->edgefilename,
".1.edge");
3694 strcat(b->vnodefilename,
".1.v.node");
3695 strcat(b->vedgefilename,
".1.v.edge");
3696 strcat(b->neighborfilename,
".1.neigh");
3697 strcat(b->offfilename,
".1.off");
3699 workstring[increment] =
'%';
3700 workstring[increment + 1] =
'd';
3701 workstring[increment + 2] =
'\0';
3702 sprintf(b->outnodefilename, workstring, meshnumber + 1);
3703 strcpy(b->outpolyfilename, b->outnodefilename);
3704 strcpy(b->outelefilename, b->outnodefilename);
3705 strcpy(b->edgefilename, b->outnodefilename);
3706 strcpy(b->vnodefilename, b->outnodefilename);
3707 strcpy(b->vedgefilename, b->outnodefilename);
3708 strcpy(b->neighborfilename, b->outnodefilename);
3709 strcpy(b->offfilename, b->outnodefilename);
3710 strcat(b->outnodefilename,
".node");
3711 strcat(b->outpolyfilename,
".poly");
3712 strcat(b->outelefilename,
".ele");
3713 strcat(b->edgefilename,
".edge");
3714 strcat(b->vnodefilename,
".v.node");
3715 strcat(b->vedgefilename,
".v.edge");
3716 strcat(b->neighborfilename,
".neigh");
3717 strcat(b->offfilename,
".off");
3719 strcat(b->innodefilename,
".node");
3720 strcat(b->inpolyfilename,
".poly");
3721 strcat(b->inelefilename,
".ele");
3722 strcat(b->areafilename,
".area");
3745 #ifdef ANSI_DECLARATORS
3748 void printtriangle(m, b, t)
3755 struct otri printtri;
3756 struct osub printsh;
3759 printf(
"triangle x%lx with orientation %d:\n", (
unsigned long) t->tri,
3761 decode(t->tri[0], printtri);
3762 if (printtri.tri == m->dummytri) {
3763 printf(
" [0] = Outer space\n");
3765 printf(
" [0] = x%lx %d\n", (
unsigned long) printtri.tri,
3768 decode(t->tri[1], printtri);
3769 if (printtri.tri == m->dummytri) {
3770 printf(
" [1] = Outer space\n");
3772 printf(
" [1] = x%lx %d\n", (
unsigned long) printtri.tri,
3775 decode(t->tri[2], printtri);
3776 if (printtri.tri == m->dummytri) {
3777 printf(
" [2] = Outer space\n");
3779 printf(
" [2] = x%lx %d\n", (
unsigned long) printtri.tri,
3783 org(*t, printvertex);
3784 if (printvertex == (vertex) NULL)
3785 printf(
" Origin[%d] = NULL\n", (t->orient + 1) % 3 + 3);
3787 printf(
" Origin[%d] = x%lx (%.12g, %.12g)\n",
3788 (t->orient + 1) % 3 + 3, (
unsigned long) printvertex,
3789 printvertex[0], printvertex[1]);
3790 dest(*t, printvertex);
3791 if (printvertex == (vertex) NULL)
3792 printf(
" Dest [%d] = NULL\n", (t->orient + 2) % 3 + 3);
3794 printf(
" Dest [%d] = x%lx (%.12g, %.12g)\n",
3795 (t->orient + 2) % 3 + 3, (
unsigned long) printvertex,
3796 printvertex[0], printvertex[1]);
3797 apex(*t, printvertex);
3798 if (printvertex == (vertex) NULL)
3799 printf(
" Apex [%d] = NULL\n", t->orient + 3);
3801 printf(
" Apex [%d] = x%lx (%.12g, %.12g)\n",
3802 t->orient + 3, (
unsigned long) printvertex,
3803 printvertex[0], printvertex[1]);
3805 if (b->usesegments) {
3806 sdecode(t->tri[6], printsh);
3807 if (printsh.ss != m->dummysub) {
3808 printf(
" [6] = x%lx %d\n", (
unsigned long) printsh.ss,
3811 sdecode(t->tri[7], printsh);
3812 if (printsh.ss != m->dummysub) {
3813 printf(
" [7] = x%lx %d\n", (
unsigned long) printsh.ss,
3816 sdecode(t->tri[8], printsh);
3817 if (printsh.ss != m->dummysub) {
3818 printf(
" [8] = x%lx %d\n", (
unsigned long) printsh.ss,
3824 printf(
" Area constraint: %.4g\n", areabound(*t));
3839 #ifdef ANSI_DECLARATORS
3842 void printsubseg(m, b, s)
3849 struct osub printsh;
3850 struct otri printtri;
3853 printf(
"subsegment x%lx with orientation %d and mark %d:\n",
3854 (
unsigned long) s->ss, s->ssorient, mark(*s));
3855 sdecode(s->ss[0], printsh);
3856 if (printsh.ss == m->dummysub) {
3857 printf(
" [0] = No subsegment\n");
3859 printf(
" [0] = x%lx %d\n", (
unsigned long) printsh.ss,
3862 sdecode(s->ss[1], printsh);
3863 if (printsh.ss == m->dummysub) {
3864 printf(
" [1] = No subsegment\n");
3866 printf(
" [1] = x%lx %d\n", (
unsigned long) printsh.ss,
3870 sorg(*s, printvertex);
3871 if (printvertex == (vertex) NULL)
3872 printf(
" Origin[%d] = NULL\n", 2 + s->ssorient);
3874 printf(
" Origin[%d] = x%lx (%.12g, %.12g)\n",
3875 2 + s->ssorient, (
unsigned long) printvertex,
3876 printvertex[0], printvertex[1]);
3877 sdest(*s, printvertex);
3878 if (printvertex == (vertex) NULL)
3879 printf(
" Dest [%d] = NULL\n", 3 - s->ssorient);
3881 printf(
" Dest [%d] = x%lx (%.12g, %.12g)\n",
3882 3 - s->ssorient, (
unsigned long) printvertex,
3883 printvertex[0], printvertex[1]);
3885 decode(s->ss[6], printtri);
3886 if (printtri.tri == m->dummytri) {
3887 printf(
" [6] = Outer space\n");
3889 printf(
" [6] = x%lx %d\n", (
unsigned long) printtri.tri,
3892 decode(s->ss[7], printtri);
3893 if (printtri.tri == m->dummytri) {
3894 printf(
" [7] = Outer space\n");
3896 printf(
" [7] = x%lx %d\n", (
unsigned long) printtri.tri,
3900 segorg(*s, printvertex);
3901 if (printvertex == (vertex) NULL)
3902 printf(
" Segment origin[%d] = NULL\n", 4 + s->ssorient);
3904 printf(
" Segment origin[%d] = x%lx (%.12g, %.12g)\n",
3905 4 + s->ssorient, (
unsigned long) printvertex,
3906 printvertex[0], printvertex[1]);
3907 segdest(*s, printvertex);
3908 if (printvertex == (vertex) NULL)
3909 printf(
" Segment dest [%d] = NULL\n", 5 - s->ssorient);
3911 printf(
" Segment dest [%d] = x%lx (%.12g, %.12g)\n",
3912 5 - s->ssorient, (
unsigned long) printvertex,
3913 printvertex[0], printvertex[1]);
3933 #ifdef ANSI_DECLARATORS
3941 pool->firstblock = (VOID **) NULL;
3942 pool->nowblock = (VOID **) NULL;
3943 pool->nextitem = (VOID *) NULL;
3944 pool->deaditemstack = (VOID *) NULL;
3945 pool->pathblock = (VOID **) NULL;
3946 pool->pathitem = (VOID *) NULL;
3947 pool->alignbytes = 0;
3948 pool->itembytes = 0;
3949 pool->itemsperblock = 0;
3950 pool->itemsfirstblock = 0;
3953 pool->unallocateditems = 0;
3954 pool->pathitemsleft = 0;
3967 #ifdef ANSI_DECLARATORS
3970 void poolrestart(pool)
3975 unsigned long alignptr;
3981 pool->nowblock = pool->firstblock;
3983 alignptr = (
unsigned long) (pool->nowblock + 1);
3985 pool->nextitem = (VOID *)
3986 (alignptr + (
unsigned long) pool->alignbytes -
3987 (alignptr % (
unsigned long) pool->alignbytes));
3989 pool->unallocateditems = pool->itemsfirstblock;
3991 pool->deaditemstack = (VOID *) NULL;
4013 #ifdef ANSI_DECLARATORS
4014 void poolinit(
struct memorypool* pool,
int bytecount,
int itemcount,
4015 int firstitemcount,
unsigned alignment)
4017 void poolinit(pool, bytecount, itemcount, firstitemcount, alignment)
4030 if (alignment >
sizeof(VOID *)) {
4031 pool->alignbytes = alignment;
4033 pool->alignbytes =
sizeof(VOID *);
4035 pool->itembytes = ((bytecount - 1) / pool->alignbytes + 1) *
4037 pool->itemsperblock = itemcount;
4038 if (firstitemcount == 0) {
4039 pool->itemsfirstblock = itemcount;
4041 pool->itemsfirstblock = firstitemcount;
4047 pool->firstblock = (VOID **)
4048 trimalloc(pool->itemsfirstblock * pool->itembytes + (
int)
sizeof(VOID *) +
4051 *(pool->firstblock) = (VOID *) NULL;
4061 #ifdef ANSI_DECLARATORS
4064 void pooldeinit(pool)
4069 while (pool->firstblock != (VOID **) NULL) {
4070 pool->nowblock = (VOID **) *(pool->firstblock);
4071 trifree((VOID *) pool->firstblock);
4072 pool->firstblock = pool->nowblock;
4082 #ifdef ANSI_DECLARATORS
4085 VOID* poolalloc(pool)
4092 unsigned long alignptr;
4096 if (pool->deaditemstack != (VOID *) NULL) {
4097 newitem = pool->deaditemstack;
4098 pool->deaditemstack = * (VOID **) pool->deaditemstack;
4101 if (pool->unallocateditems == 0) {
4103 if (*(pool->nowblock) == (VOID *) NULL) {
4105 newblock = (VOID **) trimalloc(pool->itemsperblock * pool->itembytes +
4106 (
int)
sizeof(VOID *) +
4108 *(pool->nowblock) = (VOID *) newblock;
4110 *newblock = (VOID *) NULL;
4114 pool->nowblock = (VOID **) *(pool->nowblock);
4117 alignptr = (
unsigned long) (pool->nowblock + 1);
4119 pool->nextitem = (VOID *)
4120 (alignptr + (
unsigned long) pool->alignbytes -
4121 (alignptr % (
unsigned long) pool->alignbytes));
4123 pool->unallocateditems = pool->itemsperblock;
4127 newitem = pool->nextitem;
4129 pool->nextitem = (VOID *) ((
char *) pool->nextitem + pool->itembytes);
4130 pool->unallocateditems--;
4145 #ifdef ANSI_DECLARATORS
4146 void pooldealloc(
struct memorypool* pool, VOID *dyingitem)
4148 void pooldealloc(pool, dyingitem)
4155 *((VOID **) dyingitem) = pool->deaditemstack;
4156 pool->deaditemstack = dyingitem;
4168 #ifdef ANSI_DECLARATORS
4171 void traversalinit(pool)
4176 unsigned long alignptr;
4179 pool->pathblock = pool->firstblock;
4181 alignptr = (
unsigned long) (pool->pathblock + 1);
4183 pool->pathitem = (VOID *)
4184 (alignptr + (
unsigned long) pool->alignbytes -
4185 (alignptr % (
unsigned long) pool->alignbytes));
4187 pool->pathitemsleft = pool->itemsfirstblock;
4204 #ifdef ANSI_DECLARATORS
4207 VOID *traverse(pool)
4213 unsigned long alignptr;
4216 if (pool->pathitem == pool->nextitem) {
4217 return (VOID *) NULL;
4221 if (pool->pathitemsleft == 0) {
4223 pool->pathblock = (VOID **) *(pool->pathblock);
4225 alignptr = (
unsigned long) (pool->pathblock + 1);
4227 pool->pathitem = (VOID *)
4228 (alignptr + (
unsigned long) pool->alignbytes -
4229 (alignptr % (
unsigned long) pool->alignbytes));
4231 pool->pathitemsleft = pool->itemsperblock;
4234 newitem = pool->pathitem;
4236 pool->pathitem = (VOID *) ((
char *) pool->pathitem + pool->itembytes);
4237 pool->pathitemsleft--;
4269 #ifdef ANSI_DECLARATORS
4270 void dummyinit(
struct mesh *m,
struct behavior *b,
int trianglebytes,
4273 void dummyinit(m, b, trianglebytes, subsegbytes)
4281 unsigned long alignptr;
4284 m->dummytribase = (triangle *) trimalloc(trianglebytes +
4285 m->triangles.alignbytes);
4287 alignptr = (
unsigned long) m->dummytribase;
4288 m->dummytri = (triangle *)
4289 (alignptr + (
unsigned long) m->triangles.alignbytes -
4290 (alignptr % (
unsigned long) m->triangles.alignbytes));
4295 m->dummytri[0] = (triangle) m->dummytri;
4296 m->dummytri[1] = (triangle) m->dummytri;
4297 m->dummytri[2] = (triangle) m->dummytri;
4299 m->dummytri[3] = (triangle) NULL;
4300 m->dummytri[4] = (triangle) NULL;
4301 m->dummytri[5] = (triangle) NULL;
4303 if (b->usesegments) {
4307 m->dummysubbase = (subseg *) trimalloc(subsegbytes +
4308 m->subsegs.alignbytes);
4310 alignptr = (
unsigned long) m->dummysubbase;
4311 m->dummysub = (subseg *)
4312 (alignptr + (
unsigned long) m->subsegs.alignbytes -
4313 (alignptr % (
unsigned long) m->subsegs.alignbytes));
4318 m->dummysub[0] = (subseg) m->dummysub;
4319 m->dummysub[1] = (subseg) m->dummysub;
4321 m->dummysub[2] = (subseg) NULL;
4322 m->dummysub[3] = (subseg) NULL;
4323 m->dummysub[4] = (subseg) NULL;
4324 m->dummysub[5] = (subseg) NULL;
4326 m->dummysub[6] = (subseg) m->dummytri;
4327 m->dummysub[7] = (subseg) m->dummytri;
4329 * (
int *) (m->dummysub + 8) = 0;
4333 m->dummytri[6] = (triangle) m->dummysub;
4334 m->dummytri[7] = (triangle) m->dummysub;
4335 m->dummytri[8] = (triangle) m->dummysub;
4349 #ifdef ANSI_DECLARATORS
4350 void initializevertexpool(
struct mesh *m,
struct behavior *b)
4352 void initializevertexpool(m, b)
4363 m->vertexmarkindex = ((m->mesh_dim + m->nextras) *
sizeof(REAL) +
4366 vertexsize = (m->vertexmarkindex + 2) *
sizeof(
int);
4370 m->vertex2triindex = (vertexsize +
sizeof(triangle) - 1) /
4372 vertexsize = (m->vertex2triindex + 1) *
sizeof(triangle);
4376 poolinit(&m->vertices, vertexsize, VERTEXPERBLOCK,
4377 m->invertices > VERTEXPERBLOCK ? m->invertices : VERTEXPERBLOCK,
4392 #ifdef ANSI_DECLARATORS
4393 void initializetrisubpools(
struct mesh *m,
struct behavior *b)
4395 void initializetrisubpools(m, b)
4407 m->highorderindex = 6 + (b->usesegments * 3);
4409 trisize = ((b->order + 1) * (b->order + 2) / 2 + (m->highorderindex - 3)) *
4413 m->elemattribindex = (trisize +
sizeof(REAL) - 1) /
sizeof(REAL);
4417 m->areaboundindex = m->elemattribindex + m->eextras + b->regionattrib;
4421 trisize = (m->areaboundindex + 1) *
sizeof(REAL);
4422 }
else if (m->eextras + b->regionattrib > 0) {
4423 trisize = m->areaboundindex *
sizeof(REAL);
4429 if ((b->voronoi || b->neighbors) &&
4430 (trisize < 6 *
sizeof(triangle) +
sizeof(
int))) {
4431 trisize = 6 *
sizeof(triangle) +
sizeof(
int);
4435 poolinit(&m->triangles, trisize, TRIPERBLOCK,
4436 (2 * m->invertices - 2) > TRIPERBLOCK ? (2 * m->invertices - 2) :
4439 if (b->usesegments) {
4442 poolinit(&m->subsegs, 8 *
sizeof(triangle) +
sizeof(
int),
4443 SUBSEGPERBLOCK, SUBSEGPERBLOCK, 4);
4446 dummyinit(m, b, m->triangles.itembytes, m->subsegs.itembytes);
4449 dummyinit(m, b, m->triangles.itembytes, 0);
4459 #ifdef ANSI_DECLARATORS
4460 void triangledealloc(
struct mesh *m, triangle *dyingtriangle)
4462 void triangledealloc(m, dyingtriangle)
4464 triangle *dyingtriangle;
4470 killtri(dyingtriangle);
4471 pooldealloc(&m->triangles, (VOID *) dyingtriangle);
4480 #ifdef ANSI_DECLARATORS
4481 triangle *triangletraverse(
struct mesh *m)
4483 triangle *triangletraverse(m)
4488 triangle *newtriangle;
4491 newtriangle = (triangle *) traverse(&m->triangles);
4492 if (newtriangle == (triangle *) NULL) {
4493 return (triangle *) NULL;
4495 }
while (deadtri(newtriangle));
4505 #ifdef ANSI_DECLARATORS
4506 void subsegdealloc(
struct mesh *m, subseg *dyingsubseg)
4508 void subsegdealloc(m, dyingsubseg)
4510 subseg *dyingsubseg;
4516 killsubseg(dyingsubseg);
4517 pooldealloc(&m->subsegs, (VOID *) dyingsubseg);
4526 #ifdef ANSI_DECLARATORS
4527 subseg *subsegtraverse(
struct mesh *m)
4529 subseg *subsegtraverse(m)
4537 newsubseg = (subseg *) traverse(&m->subsegs);
4538 if (newsubseg == (subseg *) NULL) {
4539 return (subseg *) NULL;
4541 }
while (deadsubseg(newsubseg));
4551 #ifdef ANSI_DECLARATORS
4552 void vertexdealloc(
struct mesh *m, vertex dyingvertex)
4554 void vertexdealloc(m, dyingvertex)
4562 setvertextype(dyingvertex, DEADVERTEX);
4563 pooldealloc(&m->vertices, (VOID *) dyingvertex);
4572 #ifdef ANSI_DECLARATORS
4573 vertex vertextraverse(
struct mesh *m)
4575 vertex vertextraverse(m)
4583 newvertex = (vertex) traverse(&m->vertices);
4584 if (newvertex == (vertex) NULL) {
4585 return (vertex) NULL;
4587 }
while (vertextype(newvertex) == DEADVERTEX);
4600 #ifdef ANSI_DECLARATORS
4601 void badsubsegdealloc(
struct mesh *m,
struct badsubseg *dyingseg)
4603 void badsubsegdealloc(m, dyingseg)
4611 dyingseg->subsegorg = (vertex) NULL;
4612 pooldealloc(&m->badsubsegs, (VOID *) dyingseg);
4625 #ifdef ANSI_DECLARATORS
4636 newseg = (
struct badsubseg *) traverse(&m->badsubsegs);
4637 if (newseg == (
struct badsubseg *) NULL) {
4640 }
while (newseg->subsegorg == (vertex) NULL);
4658 #ifdef ANSI_DECLARATORS
4659 vertex getvertex(
struct mesh *m,
struct behavior *b,
int number)
4661 vertex getvertex(m, b, number)
4670 unsigned long alignptr;
4673 getblock = m->vertices.firstblock;
4674 current = b->firstnumber;
4677 if (current + m->vertices.itemsfirstblock <= number) {
4678 getblock = (VOID **) *getblock;
4679 current += m->vertices.itemsfirstblock;
4680 while (current + m->vertices.itemsperblock <= number) {
4681 getblock = (VOID **) *getblock;
4682 current += m->vertices.itemsperblock;
4687 alignptr = (
unsigned long) (getblock + 1);
4688 foundvertex = (
char *) (alignptr + (
unsigned long) m->vertices.alignbytes -
4689 (alignptr % (
unsigned long) m->vertices.alignbytes));
4690 return (vertex) (foundvertex + m->vertices.itembytes * (number - current));
4699 #ifdef ANSI_DECLARATORS
4700 void triangledeinit(
struct mesh *m,
struct behavior *b)
4702 void triangledeinit(m, b)
4708 pooldeinit(&m->triangles);
4709 trifree((VOID *) m->dummytribase);
4710 if (b->usesegments) {
4711 pooldeinit(&m->subsegs);
4712 trifree((VOID *) m->dummysubbase);
4714 pooldeinit(&m->vertices);
4717 pooldeinit(&m->badsubsegs);
4718 if ((b->minangle > 0.0) || b->vararea || b->fixedarea || b->usertest) {
4719 pooldeinit(&m->badtriangles);
4720 pooldeinit(&m->flipstackers);
4740 #ifdef ANSI_DECLARATORS
4741 void maketriangle(
struct mesh *m,
struct behavior *b,
struct otri *newotri)
4743 void maketriangle(m, b, newotri)
4746 struct
otri *newotri;
4752 newotri->tri = (triangle *) poolalloc(&m->triangles);
4754 newotri->tri[0] = (triangle) m->dummytri;
4755 newotri->tri[1] = (triangle) m->dummytri;
4756 newotri->tri[2] = (triangle) m->dummytri;
4758 newotri->tri[3] = (triangle) NULL;
4759 newotri->tri[4] = (triangle) NULL;
4760 newotri->tri[5] = (triangle) NULL;
4761 if (b->usesegments) {
4764 newotri->tri[6] = (triangle) m->dummysub;
4765 newotri->tri[7] = (triangle) m->dummysub;
4766 newotri->tri[8] = (triangle) m->dummysub;
4768 for (i = 0; i < m->eextras; i++) {
4769 setelemattribute(*newotri, i, 0.0);
4772 setareabound(*newotri, -1.0);
4775 newotri->orient = 0;
4784 #ifdef ANSI_DECLARATORS
4785 void makesubseg(
struct mesh *m,
struct osub *newsubseg)
4787 void makesubseg(m, newsubseg)
4789 struct
osub *newsubseg;
4793 newsubseg->ss = (subseg *) poolalloc(&m->subsegs);
4796 newsubseg->ss[0] = (subseg) m->dummysub;
4797 newsubseg->ss[1] = (subseg) m->dummysub;
4799 newsubseg->ss[2] = (subseg) NULL;
4800 newsubseg->ss[3] = (subseg) NULL;
4801 newsubseg->ss[4] = (subseg) NULL;
4802 newsubseg->ss[5] = (subseg) NULL;
4804 newsubseg->ss[6] = (subseg) m->dummytri;
4805 newsubseg->ss[7] = (subseg) m->dummytri;
4807 setmark(*newsubseg, 0);
4809 newsubseg->ssorient = 0;
4833 #define Absolute(a) ((a) >= 0.0 ? (a) : -(a))
4849 #define Fast_Two_Sum_Tail(a, b, x, y) \
4853 #define Fast_Two_Sum(a, b, x, y) \
4854 x = (REAL) (a + b); \
4855 Fast_Two_Sum_Tail(a, b, x, y)
4857 #define Two_Sum_Tail(a, b, x, y) \
4858 bvirt = (REAL) (x - a); \
4859 avirt = x - bvirt; \
4860 bround = b - bvirt; \
4861 around = a - avirt; \
4864 #define Two_Sum(a, b, x, y) \
4865 x = (REAL) (a + b); \
4866 Two_Sum_Tail(a, b, x, y)
4868 #define Two_Diff_Tail(a, b, x, y) \
4869 bvirt = (REAL) (a - x); \
4870 avirt = x + bvirt; \
4871 bround = bvirt - b; \
4872 around = a - avirt; \
4875 #define Two_Diff(a, b, x, y) \
4876 x = (REAL) (a - b); \
4877 Two_Diff_Tail(a, b, x, y)
4879 #define Split(a, ahi, alo) \
4880 c = (REAL) (splitter * a); \
4881 abig = (REAL) (c - a); \
4885 #define Two_Product_Tail(a, b, x, y) \
4886 Split(a, ahi, alo); \
4887 Split(b, bhi, blo); \
4888 err1 = x - (ahi * bhi); \
4889 err2 = err1 - (alo * bhi); \
4890 err3 = err2 - (ahi * blo); \
4891 y = (alo * blo) - err3
4893 #define Two_Product(a, b, x, y) \
4894 x = (REAL) (a * b); \
4895 Two_Product_Tail(a, b, x, y)
4900 #define Two_Product_Presplit(a, b, bhi, blo, x, y) \
4901 x = (REAL) (a * b); \
4902 Split(a, ahi, alo); \
4903 err1 = x - (ahi * bhi); \
4904 err2 = err1 - (alo * bhi); \
4905 err3 = err2 - (ahi * blo); \
4906 y = (alo * blo) - err3
4910 #define Square_Tail(a, x, y) \
4911 Split(a, ahi, alo); \
4912 err1 = x - (ahi * ahi); \
4913 err3 = err1 - ((ahi + ahi) * alo); \
4914 y = (alo * alo) - err3
4916 #define Square(a, x, y) \
4917 x = (REAL) (a * a); \
4918 Square_Tail(a, x, y)
4923 #define Two_One_Sum(a1, a0, b, x2, x1, x0) \
4924 Two_Sum(a0, b , _i, x0); \
4925 Two_Sum(a1, _i, x2, x1)
4927 #define Two_One_Diff(a1, a0, b, x2, x1, x0) \
4928 Two_Diff(a0, b , _i, x0); \
4929 Two_Sum( a1, _i, x2, x1)
4931 #define Two_Two_Sum(a1, a0, b1, b0, x3, x2, x1, x0) \
4932 Two_One_Sum(a1, a0, b0, _j, _0, x0); \
4933 Two_One_Sum(_j, _0, b1, x3, x2, x1)
4935 #define Two_Two_Diff(a1, a0, b1, b0, x3, x2, x1, x0) \
4936 Two_One_Diff(a1, a0, b0, _j, _0, x0); \
4937 Two_One_Diff(_j, _0, b1, x3, x2, x1)
4941 #define Two_One_Product(a1, a0, b, x3, x2, x1, x0) \
4942 Split(b, bhi, blo); \
4943 Two_Product_Presplit(a0, b, bhi, blo, _i, x0); \
4944 Two_Product_Presplit(a1, b, bhi, blo, _j, _0); \
4945 Two_Sum(_i, _0, _k, x1); \
4946 Fast_Two_Sum(_j, _k, x3, x2)
4970 REAL check, lastcheck;
4978 _control87(_PC_24, _MCW_PC);
4980 _control87(_PC_53, _MCW_PC);
5009 every_other = !every_other;
5010 check = 1.0 + epsilon;
5011 }
while ((check != 1.0) && (check != lastcheck));
5014 resulterrbound = (3.0 + 8.0 * epsilon) * epsilon;
5015 ccwerrboundA = (3.0 + 16.0 * epsilon) * epsilon;
5016 ccwerrboundB = (2.0 + 12.0 * epsilon) * epsilon;
5017 ccwerrboundC = (9.0 + 64.0 * epsilon) * epsilon * epsilon;
5018 iccerrboundA = (10.0 + 96.0 * epsilon) * epsilon;
5019 iccerrboundB = (4.0 + 48.0 * epsilon) * epsilon;
5020 iccerrboundC = (44.0 + 576.0 * epsilon) * epsilon * epsilon;
5021 o3derrboundA = (7.0 + 56.0 * epsilon) * epsilon;
5022 o3derrboundB = (3.0 + 28.0 * epsilon) * epsilon;
5023 o3derrboundC = (26.0 + 288.0 * epsilon) * epsilon * epsilon;
5040 #ifdef ANSI_DECLARATORS
5041 int fast_expansion_sum_zeroelim(
int elen, REAL *e,
int flen, REAL *f, REAL *h)
5043 int fast_expansion_sum_zeroelim(elen, e, flen, f, h)
5056 REAL avirt, bround, around;
5057 int eindex, findex, hindex;
5062 eindex = findex = 0;
5063 if ((fnow > enow) == (fnow > -enow)) {
5071 if ((eindex < elen) && (findex < flen)) {
5072 if ((fnow > enow) == (fnow > -enow)) {
5073 Fast_Two_Sum(enow, Q, Qnew, hh);
5076 Fast_Two_Sum(fnow, Q, Qnew, hh);
5083 while ((eindex < elen) && (findex < flen)) {
5084 if ((fnow > enow) == (fnow > -enow)) {
5085 Two_Sum(Q, enow, Qnew, hh);
5088 Two_Sum(Q, fnow, Qnew, hh);
5097 while (eindex < elen) {
5098 Two_Sum(Q, enow, Qnew, hh);
5105 while (findex < flen) {
5106 Two_Sum(Q, fnow, Qnew, hh);
5113 if ((Q != 0.0) || (hindex == 0)) {
5134 #ifdef ANSI_DECLARATORS
5135 int scale_expansion_zeroelim(
int elen, REAL *e, REAL b, REAL *h)
5137 int scale_expansion_zeroelim(elen, e, b, h)
5145 INEXACT REAL Q, sum;
5147 INEXACT REAL product1;
5152 REAL avirt, bround, around;
5155 REAL ahi, alo, bhi, blo;
5156 REAL err1, err2, err3;
5159 Two_Product_Presplit(e[0], b, bhi, blo, Q, hh);
5164 for (eindex = 1; eindex < elen; eindex++) {
5166 Two_Product_Presplit(enow, b, bhi, blo, product1, product0);
5167 Two_Sum(Q, product0, sum, hh);
5171 Fast_Two_Sum(product1, sum, Q, hh);
5176 if ((Q != 0.0) || (hindex == 0)) {
5190 #ifdef ANSI_DECLARATORS
5191 REAL estimate(
int elen, REAL *e)
5193 REAL estimate(elen, e)
5203 for (eindex = 1; eindex < elen; eindex++) {
5229 #ifdef ANSI_DECLARATORS
5230 REAL counterclockwiseadapt(vertex pa, vertex pb, vertex pc, REAL detsum)
5232 REAL counterclockwiseadapt(pa, pb, pc, detsum)
5240 INEXACT REAL acx, acy, bcx, bcy;
5241 REAL acxtail, acytail, bcxtail, bcytail;
5242 INEXACT REAL detleft, detright;
5243 REAL detlefttail, detrighttail;
5245 REAL B[4], C1[8], C2[12], D[16];
5247 int C1length, C2length, Dlength;
5250 INEXACT REAL s1, t1;
5254 REAL avirt, bround, around;
5257 REAL ahi, alo, bhi, blo;
5258 REAL err1, err2, err3;
5259 INEXACT REAL _i, _j;
5262 acx = (REAL) (pa[0] - pc[0]);
5263 bcx = (REAL) (pb[0] - pc[0]);
5264 acy = (REAL) (pa[1] - pc[1]);
5265 bcy = (REAL) (pb[1] - pc[1]);
5267 Two_Product(acx, bcy, detleft, detlefttail);
5268 Two_Product(acy, bcx, detright, detrighttail);
5270 Two_Two_Diff(detleft, detlefttail, detright, detrighttail,
5271 B3, B[2], B[1], B[0]);
5274 det = estimate(4, B);
5275 errbound = ccwerrboundB * detsum;
5276 if ((det >= errbound) || (-det >= errbound)) {
5280 Two_Diff_Tail(pa[0], pc[0], acx, acxtail);
5281 Two_Diff_Tail(pb[0], pc[0], bcx, bcxtail);
5282 Two_Diff_Tail(pa[1], pc[1], acy, acytail);
5283 Two_Diff_Tail(pb[1], pc[1], bcy, bcytail);
5285 if ((acxtail == 0.0) && (acytail == 0.0)
5286 && (bcxtail == 0.0) && (bcytail == 0.0)) {
5290 errbound = ccwerrboundC * detsum + resulterrbound * Absolute(det);
5291 det += (acx * bcytail + bcy * acxtail)
5292 - (acy * bcxtail + bcx * acytail);
5293 if ((det >= errbound) || (-det >= errbound)) {
5297 Two_Product(acxtail, bcy, s1, s0);
5298 Two_Product(acytail, bcx, t1, t0);
5299 Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]);
5301 C1length = fast_expansion_sum_zeroelim(4, B, 4, u, C1);
5303 Two_Product(acx, bcytail, s1, s0);
5304 Two_Product(acy, bcxtail, t1, t0);
5305 Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]);
5307 C2length = fast_expansion_sum_zeroelim(C1length, C1, 4, u, C2);
5309 Two_Product(acxtail, bcytail, s1, s0);
5310 Two_Product(acytail, bcxtail, t1, t0);
5311 Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]);
5313 Dlength = fast_expansion_sum_zeroelim(C2length, C2, 4, u, D);
5315 return(D[Dlength - 1]);
5318 #ifdef ANSI_DECLARATORS
5319 REAL counterclockwise(
struct mesh *m,
struct behavior *b,
5320 vertex pa, vertex pb, vertex pc)
5322 REAL counterclockwise(m, b, pa, pb, pc)
5331 REAL detleft, detright, det;
5332 REAL detsum, errbound;
5334 m->counterclockcount++;
5336 detleft = (pa[0] - pc[0]) * (pb[1] - pc[1]);
5337 detright = (pa[1] - pc[1]) * (pb[0] - pc[0]);
5338 det = detleft - detright;
5344 if (detleft > 0.0) {
5345 if (detright <= 0.0) {
5348 detsum = detleft + detright;
5350 }
else if (detleft < 0.0) {
5351 if (detright >= 0.0) {
5354 detsum = -detleft - detright;
5360 errbound = ccwerrboundA * detsum;
5361 if ((det >= errbound) || (-det >= errbound)) {
5365 return counterclockwiseadapt(pa, pb, pc, detsum);
5387 #ifdef ANSI_DECLARATORS
5388 REAL incircleadapt(vertex pa, vertex pb, vertex pc, vertex pd, REAL permanent)
5390 REAL incircleadapt(pa, pb, pc, pd, permanent)
5399 INEXACT REAL adx, bdx, cdx, ady, bdy, cdy;
5402 INEXACT REAL bdxcdy1, cdxbdy1, cdxady1, adxcdy1, adxbdy1, bdxady1;
5403 REAL bdxcdy0, cdxbdy0, cdxady0, adxcdy0, adxbdy0, bdxady0;
5404 REAL bc[4], ca[4], ab[4];
5405 INEXACT REAL bc3, ca3, ab3;
5406 REAL axbc[8], axxbc[16], aybc[8], ayybc[16], adet[32];
5407 int axbclen, axxbclen, aybclen, ayybclen, alen;
5408 REAL bxca[8], bxxca[16], byca[8], byyca[16], bdet[32];
5409 int bxcalen, bxxcalen, bycalen, byycalen, blen;
5410 REAL cxab[8], cxxab[16], cyab[8], cyyab[16], cdet[32];
5411 int cxablen, cxxablen, cyablen, cyyablen, clen;
5414 REAL fin1[1152], fin2[1152];
5415 REAL *finnow, *finother, *finswap;
5418 REAL adxtail, bdxtail, cdxtail, adytail, bdytail, cdytail;
5419 INEXACT REAL adxadx1, adyady1, bdxbdx1, bdybdy1, cdxcdx1, cdycdy1;
5420 REAL adxadx0, adyady0, bdxbdx0, bdybdy0, cdxcdx0, cdycdy0;
5421 REAL aa[4], bb[4], cc[4];
5422 INEXACT REAL aa3, bb3, cc3;
5423 INEXACT REAL ti1, tj1;
5426 INEXACT REAL u3, v3;
5427 REAL temp8[8], temp16a[16], temp16b[16], temp16c[16];
5428 REAL temp32a[32], temp32b[32], temp48[48], temp64[64];
5429 int temp8len, temp16alen, temp16blen, temp16clen;
5430 int temp32alen, temp32blen, temp48len, temp64len;
5431 REAL axtbb[8], axtcc[8], aytbb[8], aytcc[8];
5432 int axtbblen, axtcclen, aytbblen, aytcclen;
5433 REAL bxtaa[8], bxtcc[8], bytaa[8], bytcc[8];
5434 int bxtaalen, bxtcclen, bytaalen, bytcclen;
5435 REAL cxtaa[8], cxtbb[8], cytaa[8], cytbb[8];
5436 int cxtaalen, cxtbblen, cytaalen, cytbblen;
5437 REAL axtbc[8], aytbc[8], bxtca[8], bytca[8], cxtab[8], cytab[8];
5438 int axtbclen=0, aytbclen=0, bxtcalen=0, bytcalen=0, cxtablen=0, cytablen=0;
5439 REAL axtbct[16], aytbct[16], bxtcat[16], bytcat[16], cxtabt[16], cytabt[16];
5440 int axtbctlen, aytbctlen, bxtcatlen, bytcatlen, cxtabtlen, cytabtlen;
5441 REAL axtbctt[8], aytbctt[8], bxtcatt[8];
5442 REAL bytcatt[8], cxtabtt[8], cytabtt[8];
5443 int axtbcttlen, aytbcttlen, bxtcattlen, bytcattlen, cxtabttlen, cytabttlen;
5444 REAL abt[8], bct[8], cat[8];
5445 int abtlen, bctlen, catlen;
5446 REAL abtt[4], bctt[4], catt[4];
5447 int abttlen, bcttlen, cattlen;
5448 INEXACT REAL abtt3, bctt3, catt3;
5452 REAL avirt, bround, around;
5455 REAL ahi, alo, bhi, blo;
5456 REAL err1, err2, err3;
5457 INEXACT REAL _i, _j;
5460 adx = (REAL) (pa[0] - pd[0]);
5461 bdx = (REAL) (pb[0] - pd[0]);
5462 cdx = (REAL) (pc[0] - pd[0]);
5463 ady = (REAL) (pa[1] - pd[1]);
5464 bdy = (REAL) (pb[1] - pd[1]);
5465 cdy = (REAL) (pc[1] - pd[1]);
5467 Two_Product(bdx, cdy, bdxcdy1, bdxcdy0);
5468 Two_Product(cdx, bdy, cdxbdy1, cdxbdy0);
5469 Two_Two_Diff(bdxcdy1, bdxcdy0, cdxbdy1, cdxbdy0, bc3, bc[2], bc[1], bc[0]);
5471 axbclen = scale_expansion_zeroelim(4, bc, adx, axbc);
5472 axxbclen = scale_expansion_zeroelim(axbclen, axbc, adx, axxbc);
5473 aybclen = scale_expansion_zeroelim(4, bc, ady, aybc);
5474 ayybclen = scale_expansion_zeroelim(aybclen, aybc, ady, ayybc);
5475 alen = fast_expansion_sum_zeroelim(axxbclen, axxbc, ayybclen, ayybc, adet);
5477 Two_Product(cdx, ady, cdxady1, cdxady0);
5478 Two_Product(adx, cdy, adxcdy1, adxcdy0);
5479 Two_Two_Diff(cdxady1, cdxady0, adxcdy1, adxcdy0, ca3, ca[2], ca[1], ca[0]);
5481 bxcalen = scale_expansion_zeroelim(4, ca, bdx, bxca);
5482 bxxcalen = scale_expansion_zeroelim(bxcalen, bxca, bdx, bxxca);
5483 bycalen = scale_expansion_zeroelim(4, ca, bdy, byca);
5484 byycalen = scale_expansion_zeroelim(bycalen, byca, bdy, byyca);
5485 blen = fast_expansion_sum_zeroelim(bxxcalen, bxxca, byycalen, byyca, bdet);
5487 Two_Product(adx, bdy, adxbdy1, adxbdy0);
5488 Two_Product(bdx, ady, bdxady1, bdxady0);
5489 Two_Two_Diff(adxbdy1, adxbdy0, bdxady1, bdxady0, ab3, ab[2], ab[1], ab[0]);
5491 cxablen = scale_expansion_zeroelim(4, ab, cdx, cxab);
5492 cxxablen = scale_expansion_zeroelim(cxablen, cxab, cdx, cxxab);
5493 cyablen = scale_expansion_zeroelim(4, ab, cdy, cyab);
5494 cyyablen = scale_expansion_zeroelim(cyablen, cyab, cdy, cyyab);
5495 clen = fast_expansion_sum_zeroelim(cxxablen, cxxab, cyyablen, cyyab, cdet);
5497 ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
5498 finlength = fast_expansion_sum_zeroelim(ablen, abdet, clen, cdet, fin1);
5500 det = estimate(finlength, fin1);
5501 errbound = iccerrboundB * permanent;
5502 if ((det >= errbound) || (-det >= errbound)) {
5506 Two_Diff_Tail(pa[0], pd[0], adx, adxtail);
5507 Two_Diff_Tail(pa[1], pd[1], ady, adytail);
5508 Two_Diff_Tail(pb[0], pd[0], bdx, bdxtail);
5509 Two_Diff_Tail(pb[1], pd[1], bdy, bdytail);
5510 Two_Diff_Tail(pc[0], pd[0], cdx, cdxtail);
5511 Two_Diff_Tail(pc[1], pd[1], cdy, cdytail);
5512 if ((adxtail == 0.0) && (bdxtail == 0.0) && (cdxtail == 0.0)
5513 && (adytail == 0.0) && (bdytail == 0.0) && (cdytail == 0.0)) {
5517 errbound = iccerrboundC * permanent + resulterrbound * Absolute(det);
5518 det += ((adx * adx + ady * ady) * ((bdx * cdytail + cdy * bdxtail)
5519 - (bdy * cdxtail + cdx * bdytail))
5520 + 2.0 * (adx * adxtail + ady * adytail) * (bdx * cdy - bdy * cdx))
5521 + ((bdx * bdx + bdy * bdy) * ((cdx * adytail + ady * cdxtail)
5522 - (cdy * adxtail + adx * cdytail))
5523 + 2.0 * (bdx * bdxtail + bdy * bdytail) * (cdx * ady - cdy * adx))
5524 + ((cdx * cdx + cdy * cdy) * ((adx * bdytail + bdy * adxtail)
5525 - (ady * bdxtail + bdx * adytail))
5526 + 2.0 * (cdx * cdxtail + cdy * cdytail) * (adx * bdy - ady * bdx));
5527 if ((det >= errbound) || (-det >= errbound)) {
5534 if ((bdxtail != 0.0) || (bdytail != 0.0)
5535 || (cdxtail != 0.0) || (cdytail != 0.0)) {
5536 Square(adx, adxadx1, adxadx0);
5537 Square(ady, adyady1, adyady0);
5538 Two_Two_Sum(adxadx1, adxadx0, adyady1, adyady0, aa3, aa[2], aa[1], aa[0]);
5541 if ((cdxtail != 0.0) || (cdytail != 0.0)
5542 || (adxtail != 0.0) || (adytail != 0.0)) {
5543 Square(bdx, bdxbdx1, bdxbdx0);
5544 Square(bdy, bdybdy1, bdybdy0);
5545 Two_Two_Sum(bdxbdx1, bdxbdx0, bdybdy1, bdybdy0, bb3, bb[2], bb[1], bb[0]);
5548 if ((adxtail != 0.0) || (adytail != 0.0)
5549 || (bdxtail != 0.0) || (bdytail != 0.0)) {
5550 Square(cdx, cdxcdx1, cdxcdx0);
5551 Square(cdy, cdycdy1, cdycdy0);
5552 Two_Two_Sum(cdxcdx1, cdxcdx0, cdycdy1, cdycdy0, cc3, cc[2], cc[1], cc[0]);
5556 if (adxtail != 0.0) {
5557 axtbclen = scale_expansion_zeroelim(4, bc, adxtail, axtbc);
5558 temp16alen = scale_expansion_zeroelim(axtbclen, axtbc, 2.0 * adx,
5561 axtcclen = scale_expansion_zeroelim(4, cc, adxtail, axtcc);
5562 temp16blen = scale_expansion_zeroelim(axtcclen, axtcc, bdy, temp16b);
5564 axtbblen = scale_expansion_zeroelim(4, bb, adxtail, axtbb);
5565 temp16clen = scale_expansion_zeroelim(axtbblen, axtbb, -cdy, temp16c);
5567 temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
5568 temp16blen, temp16b, temp32a);
5569 temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
5570 temp32alen, temp32a, temp48);
5571 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
5573 finswap = finnow; finnow = finother; finother = finswap;
5575 if (adytail != 0.0) {
5576 aytbclen = scale_expansion_zeroelim(4, bc, adytail, aytbc);
5577 temp16alen = scale_expansion_zeroelim(aytbclen, aytbc, 2.0 * ady,
5580 aytbblen = scale_expansion_zeroelim(4, bb, adytail, aytbb);
5581 temp16blen = scale_expansion_zeroelim(aytbblen, aytbb, cdx, temp16b);
5583 aytcclen = scale_expansion_zeroelim(4, cc, adytail, aytcc);
5584 temp16clen = scale_expansion_zeroelim(aytcclen, aytcc, -bdx, temp16c);
5586 temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
5587 temp16blen, temp16b, temp32a);
5588 temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
5589 temp32alen, temp32a, temp48);
5590 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
5592 finswap = finnow; finnow = finother; finother = finswap;
5594 if (bdxtail != 0.0) {
5595 bxtcalen = scale_expansion_zeroelim(4, ca, bdxtail, bxtca);
5596 temp16alen = scale_expansion_zeroelim(bxtcalen, bxtca, 2.0 * bdx,
5599 bxtaalen = scale_expansion_zeroelim(4, aa, bdxtail, bxtaa);
5600 temp16blen = scale_expansion_zeroelim(bxtaalen, bxtaa, cdy, temp16b);
5602 bxtcclen = scale_expansion_zeroelim(4, cc, bdxtail, bxtcc);
5603 temp16clen = scale_expansion_zeroelim(bxtcclen, bxtcc, -ady, temp16c);
5605 temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
5606 temp16blen, temp16b, temp32a);
5607 temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
5608 temp32alen, temp32a, temp48);
5609 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
5611 finswap = finnow; finnow = finother; finother = finswap;
5613 if (bdytail != 0.0) {
5614 bytcalen = scale_expansion_zeroelim(4, ca, bdytail, bytca);
5615 temp16alen = scale_expansion_zeroelim(bytcalen, bytca, 2.0 * bdy,
5618 bytcclen = scale_expansion_zeroelim(4, cc, bdytail, bytcc);
5619 temp16blen = scale_expansion_zeroelim(bytcclen, bytcc, adx, temp16b);
5621 bytaalen = scale_expansion_zeroelim(4, aa, bdytail, bytaa);
5622 temp16clen = scale_expansion_zeroelim(bytaalen, bytaa, -cdx, temp16c);
5624 temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
5625 temp16blen, temp16b, temp32a);
5626 temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
5627 temp32alen, temp32a, temp48);
5628 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
5630 finswap = finnow; finnow = finother; finother = finswap;
5632 if (cdxtail != 0.0) {
5633 cxtablen = scale_expansion_zeroelim(4, ab, cdxtail, cxtab);
5634 temp16alen = scale_expansion_zeroelim(cxtablen, cxtab, 2.0 * cdx,
5637 cxtbblen = scale_expansion_zeroelim(4, bb, cdxtail, cxtbb);
5638 temp16blen = scale_expansion_zeroelim(cxtbblen, cxtbb, ady, temp16b);
5640 cxtaalen = scale_expansion_zeroelim(4, aa, cdxtail, cxtaa);
5641 temp16clen = scale_expansion_zeroelim(cxtaalen, cxtaa, -bdy, temp16c);
5643 temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
5644 temp16blen, temp16b, temp32a);
5645 temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
5646 temp32alen, temp32a, temp48);
5647 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
5649 finswap = finnow; finnow = finother; finother = finswap;
5651 if (cdytail != 0.0) {
5652 cytablen = scale_expansion_zeroelim(4, ab, cdytail, cytab);
5653 temp16alen = scale_expansion_zeroelim(cytablen, cytab, 2.0 * cdy,
5656 cytaalen = scale_expansion_zeroelim(4, aa, cdytail, cytaa);
5657 temp16blen = scale_expansion_zeroelim(cytaalen, cytaa, bdx, temp16b);
5659 cytbblen = scale_expansion_zeroelim(4, bb, cdytail, cytbb);
5660 temp16clen = scale_expansion_zeroelim(cytbblen, cytbb, -adx, temp16c);
5662 temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
5663 temp16blen, temp16b, temp32a);
5664 temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c,
5665 temp32alen, temp32a, temp48);
5666 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
5668 finswap = finnow; finnow = finother; finother = finswap;
5671 if ((adxtail != 0.0) || (adytail != 0.0)) {
5672 if ((bdxtail != 0.0) || (bdytail != 0.0)
5673 || (cdxtail != 0.0) || (cdytail != 0.0)) {
5674 Two_Product(bdxtail, cdy, ti1, ti0);
5675 Two_Product(bdx, cdytail, tj1, tj0);
5676 Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]);
5679 Two_Product(cdxtail, negate, ti1, ti0);
5681 Two_Product(cdx, negate, tj1, tj0);
5682 Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]);
5684 bctlen = fast_expansion_sum_zeroelim(4, u, 4, v, bct);
5686 Two_Product(bdxtail, cdytail, ti1, ti0);
5687 Two_Product(cdxtail, bdytail, tj1, tj0);
5688 Two_Two_Diff(ti1, ti0, tj1, tj0, bctt3, bctt[2], bctt[1], bctt[0]);
5698 if (adxtail != 0.0) {
5699 temp16alen = scale_expansion_zeroelim(axtbclen, axtbc, adxtail, temp16a);
5700 axtbctlen = scale_expansion_zeroelim(bctlen, bct, adxtail, axtbct);
5701 temp32alen = scale_expansion_zeroelim(axtbctlen, axtbct, 2.0 * adx,
5703 temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
5704 temp32alen, temp32a, temp48);
5705 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
5707 finswap = finnow; finnow = finother; finother = finswap;
5708 if (bdytail != 0.0) {
5709 temp8len = scale_expansion_zeroelim(4, cc, adxtail, temp8);
5710 temp16alen = scale_expansion_zeroelim(temp8len, temp8, bdytail,
5712 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
5714 finswap = finnow; finnow = finother; finother = finswap;
5716 if (cdytail != 0.0) {
5717 temp8len = scale_expansion_zeroelim(4, bb, -adxtail, temp8);
5718 temp16alen = scale_expansion_zeroelim(temp8len, temp8, cdytail,
5720 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
5722 finswap = finnow; finnow = finother; finother = finswap;
5725 temp32alen = scale_expansion_zeroelim(axtbctlen, axtbct, adxtail,
5727 axtbcttlen = scale_expansion_zeroelim(bcttlen, bctt, adxtail, axtbctt);
5728 temp16alen = scale_expansion_zeroelim(axtbcttlen, axtbctt, 2.0 * adx,
5730 temp16blen = scale_expansion_zeroelim(axtbcttlen, axtbctt, adxtail,
5732 temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
5733 temp16blen, temp16b, temp32b);
5734 temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
5735 temp32blen, temp32b, temp64);
5736 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
5738 finswap = finnow; finnow = finother; finother = finswap;
5740 if (adytail != 0.0) {
5741 temp16alen = scale_expansion_zeroelim(aytbclen, aytbc, adytail, temp16a);
5742 aytbctlen = scale_expansion_zeroelim(bctlen, bct, adytail, aytbct);
5743 temp32alen = scale_expansion_zeroelim(aytbctlen, aytbct, 2.0 * ady,
5745 temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
5746 temp32alen, temp32a, temp48);
5747 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
5749 finswap = finnow; finnow = finother; finother = finswap;
5752 temp32alen = scale_expansion_zeroelim(aytbctlen, aytbct, adytail,
5754 aytbcttlen = scale_expansion_zeroelim(bcttlen, bctt, adytail, aytbctt);
5755 temp16alen = scale_expansion_zeroelim(aytbcttlen, aytbctt, 2.0 * ady,
5757 temp16blen = scale_expansion_zeroelim(aytbcttlen, aytbctt, adytail,
5759 temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
5760 temp16blen, temp16b, temp32b);
5761 temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
5762 temp32blen, temp32b, temp64);
5763 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
5765 finswap = finnow; finnow = finother; finother = finswap;
5768 if ((bdxtail != 0.0) || (bdytail != 0.0)) {
5769 if ((cdxtail != 0.0) || (cdytail != 0.0)
5770 || (adxtail != 0.0) || (adytail != 0.0)) {
5771 Two_Product(cdxtail, ady, ti1, ti0);
5772 Two_Product(cdx, adytail, tj1, tj0);
5773 Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]);
5776 Two_Product(adxtail, negate, ti1, ti0);
5778 Two_Product(adx, negate, tj1, tj0);
5779 Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]);
5781 catlen = fast_expansion_sum_zeroelim(4, u, 4, v, cat);
5783 Two_Product(cdxtail, adytail, ti1, ti0);
5784 Two_Product(adxtail, cdytail, tj1, tj0);
5785 Two_Two_Diff(ti1, ti0, tj1, tj0, catt3, catt[2], catt[1], catt[0]);
5795 if (bdxtail != 0.0) {
5796 temp16alen = scale_expansion_zeroelim(bxtcalen, bxtca, bdxtail, temp16a);
5797 bxtcatlen = scale_expansion_zeroelim(catlen, cat, bdxtail, bxtcat);
5798 temp32alen = scale_expansion_zeroelim(bxtcatlen, bxtcat, 2.0 * bdx,
5800 temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
5801 temp32alen, temp32a, temp48);
5802 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
5804 finswap = finnow; finnow = finother; finother = finswap;
5805 if (cdytail != 0.0) {
5806 temp8len = scale_expansion_zeroelim(4, aa, bdxtail, temp8);
5807 temp16alen = scale_expansion_zeroelim(temp8len, temp8, cdytail,
5809 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
5811 finswap = finnow; finnow = finother; finother = finswap;
5813 if (adytail != 0.0) {
5814 temp8len = scale_expansion_zeroelim(4, cc, -bdxtail, temp8);
5815 temp16alen = scale_expansion_zeroelim(temp8len, temp8, adytail,
5817 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
5819 finswap = finnow; finnow = finother; finother = finswap;
5822 temp32alen = scale_expansion_zeroelim(bxtcatlen, bxtcat, bdxtail,
5824 bxtcattlen = scale_expansion_zeroelim(cattlen, catt, bdxtail, bxtcatt);
5825 temp16alen = scale_expansion_zeroelim(bxtcattlen, bxtcatt, 2.0 * bdx,
5827 temp16blen = scale_expansion_zeroelim(bxtcattlen, bxtcatt, bdxtail,
5829 temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
5830 temp16blen, temp16b, temp32b);
5831 temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
5832 temp32blen, temp32b, temp64);
5833 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
5835 finswap = finnow; finnow = finother; finother = finswap;
5837 if (bdytail != 0.0) {
5838 temp16alen = scale_expansion_zeroelim(bytcalen, bytca, bdytail, temp16a);
5839 bytcatlen = scale_expansion_zeroelim(catlen, cat, bdytail, bytcat);
5840 temp32alen = scale_expansion_zeroelim(bytcatlen, bytcat, 2.0 * bdy,
5842 temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
5843 temp32alen, temp32a, temp48);
5844 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
5846 finswap = finnow; finnow = finother; finother = finswap;
5849 temp32alen = scale_expansion_zeroelim(bytcatlen, bytcat, bdytail,
5851 bytcattlen = scale_expansion_zeroelim(cattlen, catt, bdytail, bytcatt);
5852 temp16alen = scale_expansion_zeroelim(bytcattlen, bytcatt, 2.0 * bdy,
5854 temp16blen = scale_expansion_zeroelim(bytcattlen, bytcatt, bdytail,
5856 temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
5857 temp16blen, temp16b, temp32b);
5858 temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
5859 temp32blen, temp32b, temp64);
5860 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
5862 finswap = finnow; finnow = finother; finother = finswap;
5865 if ((cdxtail != 0.0) || (cdytail != 0.0)) {
5866 if ((adxtail != 0.0) || (adytail != 0.0)
5867 || (bdxtail != 0.0) || (bdytail != 0.0)) {
5868 Two_Product(adxtail, bdy, ti1, ti0);
5869 Two_Product(adx, bdytail, tj1, tj0);
5870 Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]);
5873 Two_Product(bdxtail, negate, ti1, ti0);
5875 Two_Product(bdx, negate, tj1, tj0);
5876 Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]);
5878 abtlen = fast_expansion_sum_zeroelim(4, u, 4, v, abt);
5880 Two_Product(adxtail, bdytail, ti1, ti0);
5881 Two_Product(bdxtail, adytail, tj1, tj0);
5882 Two_Two_Diff(ti1, ti0, tj1, tj0, abtt3, abtt[2], abtt[1], abtt[0]);
5892 if (cdxtail != 0.0) {
5893 temp16alen = scale_expansion_zeroelim(cxtablen, cxtab, cdxtail, temp16a);
5894 cxtabtlen = scale_expansion_zeroelim(abtlen, abt, cdxtail, cxtabt);
5895 temp32alen = scale_expansion_zeroelim(cxtabtlen, cxtabt, 2.0 * cdx,
5897 temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
5898 temp32alen, temp32a, temp48);
5899 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
5901 finswap = finnow; finnow = finother; finother = finswap;
5902 if (adytail != 0.0) {
5903 temp8len = scale_expansion_zeroelim(4, bb, cdxtail, temp8);
5904 temp16alen = scale_expansion_zeroelim(temp8len, temp8, adytail,
5906 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
5908 finswap = finnow; finnow = finother; finother = finswap;
5910 if (bdytail != 0.0) {
5911 temp8len = scale_expansion_zeroelim(4, aa, -cdxtail, temp8);
5912 temp16alen = scale_expansion_zeroelim(temp8len, temp8, bdytail,
5914 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
5916 finswap = finnow; finnow = finother; finother = finswap;
5919 temp32alen = scale_expansion_zeroelim(cxtabtlen, cxtabt, cdxtail,
5921 cxtabttlen = scale_expansion_zeroelim(abttlen, abtt, cdxtail, cxtabtt);
5922 temp16alen = scale_expansion_zeroelim(cxtabttlen, cxtabtt, 2.0 * cdx,
5924 temp16blen = scale_expansion_zeroelim(cxtabttlen, cxtabtt, cdxtail,
5926 temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
5927 temp16blen, temp16b, temp32b);
5928 temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
5929 temp32blen, temp32b, temp64);
5930 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
5932 finswap = finnow; finnow = finother; finother = finswap;
5934 if (cdytail != 0.0) {
5935 temp16alen = scale_expansion_zeroelim(cytablen, cytab, cdytail, temp16a);
5936 cytabtlen = scale_expansion_zeroelim(abtlen, abt, cdytail, cytabt);
5937 temp32alen = scale_expansion_zeroelim(cytabtlen, cytabt, 2.0 * cdy,
5939 temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
5940 temp32alen, temp32a, temp48);
5941 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
5943 finswap = finnow; finnow = finother; finother = finswap;
5946 temp32alen = scale_expansion_zeroelim(cytabtlen, cytabt, cdytail,
5948 cytabttlen = scale_expansion_zeroelim(abttlen, abtt, cdytail, cytabtt);
5949 temp16alen = scale_expansion_zeroelim(cytabttlen, cytabtt, 2.0 * cdy,
5951 temp16blen = scale_expansion_zeroelim(cytabttlen, cytabtt, cdytail,
5953 temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a,
5954 temp16blen, temp16b, temp32b);
5955 temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a,
5956 temp32blen, temp32b, temp64);
5957 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len,
5959 finswap = finnow; finnow = finother; finother = finswap;
5963 return finnow[finlength - 1];
5966 #ifdef ANSI_DECLARATORS
5968 vertex pa, vertex pb, vertex pc, vertex pd)
5970 REAL incircle(m, b, pa, pb, pc, pd)
5980 REAL adx, bdx, cdx, ady, bdy, cdy;
5981 REAL bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady;
5982 REAL alift, blift, clift;
5984 REAL permanent, errbound;
5988 adx = pa[0] - pd[0];
5989 bdx = pb[0] - pd[0];
5990 cdx = pc[0] - pd[0];
5991 ady = pa[1] - pd[1];
5992 bdy = pb[1] - pd[1];
5993 cdy = pc[1] - pd[1];
5997 alift = adx * adx + ady * ady;
6001 blift = bdx * bdx + bdy * bdy;
6005 clift = cdx * cdx + cdy * cdy;
6007 det = alift * (bdxcdy - cdxbdy)
6008 + blift * (cdxady - adxcdy)
6009 + clift * (adxbdy - bdxady);
6015 permanent = (Absolute(bdxcdy) + Absolute(cdxbdy)) * alift
6016 + (Absolute(cdxady) + Absolute(adxcdy)) * blift
6017 + (Absolute(adxbdy) + Absolute(bdxady)) * clift;
6018 errbound = iccerrboundA * permanent;
6019 if ((det > errbound) || (-det > errbound)) {
6023 return incircleadapt(pa, pb, pc, pd, permanent);
6048 #ifdef ANSI_DECLARATORS
6049 REAL orient3dadapt(vertex pa, vertex pb, vertex pc, vertex pd,
6050 REAL aheight, REAL bheight, REAL cheight, REAL dheight,
6053 REAL orient3dadapt(pa, pb, pc, pd,
6054 aheight, bheight, cheight, dheight, permanent)
6067 INEXACT REAL adx, bdx, cdx, ady, bdy, cdy, adheight, bdheight, cdheight;
6070 INEXACT REAL bdxcdy1, cdxbdy1, cdxady1, adxcdy1, adxbdy1, bdxady1;
6071 REAL bdxcdy0, cdxbdy0, cdxady0, adxcdy0, adxbdy0, bdxady0;
6072 REAL bc[4], ca[4], ab[4];
6073 INEXACT REAL bc3, ca3, ab3;
6074 REAL adet[8], bdet[8], cdet[8];
6075 int alen, blen, clen;
6078 REAL *finnow, *finother, *finswap;
6079 REAL fin1[192], fin2[192];
6082 REAL adxtail, bdxtail, cdxtail;
6083 REAL adytail, bdytail, cdytail;
6084 REAL adheighttail, bdheighttail, cdheighttail;
6085 INEXACT REAL at_blarge, at_clarge;
6086 INEXACT REAL bt_clarge, bt_alarge;
6087 INEXACT REAL ct_alarge, ct_blarge;
6088 REAL at_b[4], at_c[4], bt_c[4], bt_a[4], ct_a[4], ct_b[4];
6089 int at_blen, at_clen, bt_clen, bt_alen, ct_alen, ct_blen;
6090 INEXACT REAL bdxt_cdy1, cdxt_bdy1, cdxt_ady1;
6091 INEXACT REAL adxt_cdy1, adxt_bdy1, bdxt_ady1;
6092 REAL bdxt_cdy0, cdxt_bdy0, cdxt_ady0;
6093 REAL adxt_cdy0, adxt_bdy0, bdxt_ady0;
6094 INEXACT REAL bdyt_cdx1, cdyt_bdx1, cdyt_adx1;
6095 INEXACT REAL adyt_cdx1, adyt_bdx1, bdyt_adx1;
6096 REAL bdyt_cdx0, cdyt_bdx0, cdyt_adx0;
6097 REAL adyt_cdx0, adyt_bdx0, bdyt_adx0;
6098 REAL bct[8], cat[8], abt[8];
6099 int bctlen, catlen, abtlen;
6100 INEXACT REAL bdxt_cdyt1, cdxt_bdyt1, cdxt_adyt1;
6101 INEXACT REAL adxt_cdyt1, adxt_bdyt1, bdxt_adyt1;
6102 REAL bdxt_cdyt0, cdxt_bdyt0, cdxt_adyt0;
6103 REAL adxt_cdyt0, adxt_bdyt0, bdxt_adyt0;
6104 REAL u[4], v[12], w[16];
6106 int vlength, wlength;
6110 REAL avirt, bround, around;
6113 REAL ahi, alo, bhi, blo;
6114 REAL err1, err2, err3;
6115 INEXACT REAL _i, _j, _k;
6118 adx = (REAL) (pa[0] - pd[0]);
6119 bdx = (REAL) (pb[0] - pd[0]);
6120 cdx = (REAL) (pc[0] - pd[0]);
6121 ady = (REAL) (pa[1] - pd[1]);
6122 bdy = (REAL) (pb[1] - pd[1]);
6123 cdy = (REAL) (pc[1] - pd[1]);
6124 adheight = (REAL) (aheight - dheight);
6125 bdheight = (REAL) (bheight - dheight);
6126 cdheight = (REAL) (cheight - dheight);
6128 Two_Product(bdx, cdy, bdxcdy1, bdxcdy0);
6129 Two_Product(cdx, bdy, cdxbdy1, cdxbdy0);
6130 Two_Two_Diff(bdxcdy1, bdxcdy0, cdxbdy1, cdxbdy0, bc3, bc[2], bc[1], bc[0]);
6132 alen = scale_expansion_zeroelim(4, bc, adheight, adet);
6134 Two_Product(cdx, ady, cdxady1, cdxady0);
6135 Two_Product(adx, cdy, adxcdy1, adxcdy0);
6136 Two_Two_Diff(cdxady1, cdxady0, adxcdy1, adxcdy0, ca3, ca[2], ca[1], ca[0]);
6138 blen = scale_expansion_zeroelim(4, ca, bdheight, bdet);
6140 Two_Product(adx, bdy, adxbdy1, adxbdy0);
6141 Two_Product(bdx, ady, bdxady1, bdxady0);
6142 Two_Two_Diff(adxbdy1, adxbdy0, bdxady1, bdxady0, ab3, ab[2], ab[1], ab[0]);
6144 clen = scale_expansion_zeroelim(4, ab, cdheight, cdet);
6146 ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
6147 finlength = fast_expansion_sum_zeroelim(ablen, abdet, clen, cdet, fin1);
6149 det = estimate(finlength, fin1);
6150 errbound = o3derrboundB * permanent;
6151 if ((det >= errbound) || (-det >= errbound)) {
6155 Two_Diff_Tail(pa[0], pd[0], adx, adxtail);
6156 Two_Diff_Tail(pb[0], pd[0], bdx, bdxtail);
6157 Two_Diff_Tail(pc[0], pd[0], cdx, cdxtail);
6158 Two_Diff_Tail(pa[1], pd[1], ady, adytail);
6159 Two_Diff_Tail(pb[1], pd[1], bdy, bdytail);
6160 Two_Diff_Tail(pc[1], pd[1], cdy, cdytail);
6161 Two_Diff_Tail(aheight, dheight, adheight, adheighttail);
6162 Two_Diff_Tail(bheight, dheight, bdheight, bdheighttail);
6163 Two_Diff_Tail(cheight, dheight, cdheight, cdheighttail);
6165 if ((adxtail == 0.0) && (bdxtail == 0.0) && (cdxtail == 0.0) &&
6166 (adytail == 0.0) && (bdytail == 0.0) && (cdytail == 0.0) &&
6167 (adheighttail == 0.0) &&
6168 (bdheighttail == 0.0) &&
6169 (cdheighttail == 0.0)) {
6173 errbound = o3derrboundC * permanent + resulterrbound * Absolute(det);
6174 det += (adheight * ((bdx * cdytail + cdy * bdxtail) -
6175 (bdy * cdxtail + cdx * bdytail)) +
6176 adheighttail * (bdx * cdy - bdy * cdx)) +
6177 (bdheight * ((cdx * adytail + ady * cdxtail) -
6178 (cdy * adxtail + adx * cdytail)) +
6179 bdheighttail * (cdx * ady - cdy * adx)) +
6180 (cdheight * ((adx * bdytail + bdy * adxtail) -
6181 (ady * bdxtail + bdx * adytail)) +
6182 cdheighttail * (adx * bdy - ady * bdx));
6183 if ((det >= errbound) || (-det >= errbound)) {
6190 if (adxtail == 0.0) {
6191 if (adytail == 0.0) {
6198 Two_Product(negate, bdx, at_blarge, at_b[0]);
6199 at_b[1] = at_blarge;
6201 Two_Product(adytail, cdx, at_clarge, at_c[0]);
6202 at_c[1] = at_clarge;
6206 if (adytail == 0.0) {
6207 Two_Product(adxtail, bdy, at_blarge, at_b[0]);
6208 at_b[1] = at_blarge;
6211 Two_Product(negate, cdy, at_clarge, at_c[0]);
6212 at_c[1] = at_clarge;
6215 Two_Product(adxtail, bdy, adxt_bdy1, adxt_bdy0);
6216 Two_Product(adytail, bdx, adyt_bdx1, adyt_bdx0);
6217 Two_Two_Diff(adxt_bdy1, adxt_bdy0, adyt_bdx1, adyt_bdx0,
6218 at_blarge, at_b[2], at_b[1], at_b[0]);
6219 at_b[3] = at_blarge;
6221 Two_Product(adytail, cdx, adyt_cdx1, adyt_cdx0);
6222 Two_Product(adxtail, cdy, adxt_cdy1, adxt_cdy0);
6223 Two_Two_Diff(adyt_cdx1, adyt_cdx0, adxt_cdy1, adxt_cdy0,
6224 at_clarge, at_c[2], at_c[1], at_c[0]);
6225 at_c[3] = at_clarge;
6229 if (bdxtail == 0.0) {
6230 if (bdytail == 0.0) {
6237 Two_Product(negate, cdx, bt_clarge, bt_c[0]);
6238 bt_c[1] = bt_clarge;
6240 Two_Product(bdytail, adx, bt_alarge, bt_a[0]);
6241 bt_a[1] = bt_alarge;
6245 if (bdytail == 0.0) {
6246 Two_Product(bdxtail, cdy, bt_clarge, bt_c[0]);
6247 bt_c[1] = bt_clarge;
6250 Two_Product(negate, ady, bt_alarge, bt_a[0]);
6251 bt_a[1] = bt_alarge;
6254 Two_Product(bdxtail, cdy, bdxt_cdy1, bdxt_cdy0);
6255 Two_Product(bdytail, cdx, bdyt_cdx1, bdyt_cdx0);
6256 Two_Two_Diff(bdxt_cdy1, bdxt_cdy0, bdyt_cdx1, bdyt_cdx0,
6257 bt_clarge, bt_c[2], bt_c[1], bt_c[0]);
6258 bt_c[3] = bt_clarge;
6260 Two_Product(bdytail, adx, bdyt_adx1, bdyt_adx0);
6261 Two_Product(bdxtail, ady, bdxt_ady1, bdxt_ady0);
6262 Two_Two_Diff(bdyt_adx1, bdyt_adx0, bdxt_ady1, bdxt_ady0,
6263 bt_alarge, bt_a[2], bt_a[1], bt_a[0]);
6264 bt_a[3] = bt_alarge;
6268 if (cdxtail == 0.0) {
6269 if (cdytail == 0.0) {
6276 Two_Product(negate, adx, ct_alarge, ct_a[0]);
6277 ct_a[1] = ct_alarge;
6279 Two_Product(cdytail, bdx, ct_blarge, ct_b[0]);
6280 ct_b[1] = ct_blarge;
6284 if (cdytail == 0.0) {
6285 Two_Product(cdxtail, ady, ct_alarge, ct_a[0]);
6286 ct_a[1] = ct_alarge;
6289 Two_Product(negate, bdy, ct_blarge, ct_b[0]);
6290 ct_b[1] = ct_blarge;
6293 Two_Product(cdxtail, ady, cdxt_ady1, cdxt_ady0);
6294 Two_Product(cdytail, adx, cdyt_adx1, cdyt_adx0);
6295 Two_Two_Diff(cdxt_ady1, cdxt_ady0, cdyt_adx1, cdyt_adx0,
6296 ct_alarge, ct_a[2], ct_a[1], ct_a[0]);
6297 ct_a[3] = ct_alarge;
6299 Two_Product(cdytail, bdx, cdyt_bdx1, cdyt_bdx0);
6300 Two_Product(cdxtail, bdy, cdxt_bdy1, cdxt_bdy0);
6301 Two_Two_Diff(cdyt_bdx1, cdyt_bdx0, cdxt_bdy1, cdxt_bdy0,
6302 ct_blarge, ct_b[2], ct_b[1], ct_b[0]);
6303 ct_b[3] = ct_blarge;
6308 bctlen = fast_expansion_sum_zeroelim(bt_clen, bt_c, ct_blen, ct_b, bct);
6309 wlength = scale_expansion_zeroelim(bctlen, bct, adheight, w);
6310 finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
6312 finswap = finnow; finnow = finother; finother = finswap;
6314 catlen = fast_expansion_sum_zeroelim(ct_alen, ct_a, at_clen, at_c, cat);
6315 wlength = scale_expansion_zeroelim(catlen, cat, bdheight, w);
6316 finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
6318 finswap = finnow; finnow = finother; finother = finswap;
6320 abtlen = fast_expansion_sum_zeroelim(at_blen, at_b, bt_alen, bt_a, abt);
6321 wlength = scale_expansion_zeroelim(abtlen, abt, cdheight, w);
6322 finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
6324 finswap = finnow; finnow = finother; finother = finswap;
6326 if (adheighttail != 0.0) {
6327 vlength = scale_expansion_zeroelim(4, bc, adheighttail, v);
6328 finlength = fast_expansion_sum_zeroelim(finlength, finnow, vlength, v,
6330 finswap = finnow; finnow = finother; finother = finswap;
6332 if (bdheighttail != 0.0) {
6333 vlength = scale_expansion_zeroelim(4, ca, bdheighttail, v);
6334 finlength = fast_expansion_sum_zeroelim(finlength, finnow, vlength, v,
6336 finswap = finnow; finnow = finother; finother = finswap;
6338 if (cdheighttail != 0.0) {
6339 vlength = scale_expansion_zeroelim(4, ab, cdheighttail, v);
6340 finlength = fast_expansion_sum_zeroelim(finlength, finnow, vlength, v,
6342 finswap = finnow; finnow = finother; finother = finswap;
6345 if (adxtail != 0.0) {
6346 if (bdytail != 0.0) {
6347 Two_Product(adxtail, bdytail, adxt_bdyt1, adxt_bdyt0);
6348 Two_One_Product(adxt_bdyt1, adxt_bdyt0, cdheight, u3, u[2], u[1], u[0]);
6350 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
6352 finswap = finnow; finnow = finother; finother = finswap;
6353 if (cdheighttail != 0.0) {
6354 Two_One_Product(adxt_bdyt1, adxt_bdyt0, cdheighttail,
6355 u3, u[2], u[1], u[0]);
6357 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
6359 finswap = finnow; finnow = finother; finother = finswap;
6362 if (cdytail != 0.0) {
6364 Two_Product(negate, cdytail, adxt_cdyt1, adxt_cdyt0);
6365 Two_One_Product(adxt_cdyt1, adxt_cdyt0, bdheight, u3, u[2], u[1], u[0]);
6367 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
6369 finswap = finnow; finnow = finother; finother = finswap;
6370 if (bdheighttail != 0.0) {
6371 Two_One_Product(adxt_cdyt1, adxt_cdyt0, bdheighttail,
6372 u3, u[2], u[1], u[0]);
6374 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
6376 finswap = finnow; finnow = finother; finother = finswap;
6380 if (bdxtail != 0.0) {
6381 if (cdytail != 0.0) {
6382 Two_Product(bdxtail, cdytail, bdxt_cdyt1, bdxt_cdyt0);
6383 Two_One_Product(bdxt_cdyt1, bdxt_cdyt0, adheight, u3, u[2], u[1], u[0]);
6385 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
6387 finswap = finnow; finnow = finother; finother = finswap;
6388 if (adheighttail != 0.0) {
6389 Two_One_Product(bdxt_cdyt1, bdxt_cdyt0, adheighttail,
6390 u3, u[2], u[1], u[0]);
6392 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
6394 finswap = finnow; finnow = finother; finother = finswap;
6397 if (adytail != 0.0) {
6399 Two_Product(negate, adytail, bdxt_adyt1, bdxt_adyt0);
6400 Two_One_Product(bdxt_adyt1, bdxt_adyt0, cdheight, u3, u[2], u[1], u[0]);
6402 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
6404 finswap = finnow; finnow = finother; finother = finswap;
6405 if (cdheighttail != 0.0) {
6406 Two_One_Product(bdxt_adyt1, bdxt_adyt0, cdheighttail,
6407 u3, u[2], u[1], u[0]);
6409 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
6411 finswap = finnow; finnow = finother; finother = finswap;
6415 if (cdxtail != 0.0) {
6416 if (adytail != 0.0) {
6417 Two_Product(cdxtail, adytail, cdxt_adyt1, cdxt_adyt0);
6418 Two_One_Product(cdxt_adyt1, cdxt_adyt0, bdheight, u3, u[2], u[1], u[0]);
6420 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
6422 finswap = finnow; finnow = finother; finother = finswap;
6423 if (bdheighttail != 0.0) {
6424 Two_One_Product(cdxt_adyt1, cdxt_adyt0, bdheighttail,
6425 u3, u[2], u[1], u[0]);
6427 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
6429 finswap = finnow; finnow = finother; finother = finswap;
6432 if (bdytail != 0.0) {
6434 Two_Product(negate, bdytail, cdxt_bdyt1, cdxt_bdyt0);
6435 Two_One_Product(cdxt_bdyt1, cdxt_bdyt0, adheight, u3, u[2], u[1], u[0]);
6437 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
6439 finswap = finnow; finnow = finother; finother = finswap;
6440 if (adheighttail != 0.0) {
6441 Two_One_Product(cdxt_bdyt1, cdxt_bdyt0, adheighttail,
6442 u3, u[2], u[1], u[0]);
6444 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
6446 finswap = finnow; finnow = finother; finother = finswap;
6451 if (adheighttail != 0.0) {
6452 wlength = scale_expansion_zeroelim(bctlen, bct, adheighttail, w);
6453 finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
6455 finswap = finnow; finnow = finother; finother = finswap;
6457 if (bdheighttail != 0.0) {
6458 wlength = scale_expansion_zeroelim(catlen, cat, bdheighttail, w);
6459 finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
6461 finswap = finnow; finnow = finother; finother = finswap;
6463 if (cdheighttail != 0.0) {
6464 wlength = scale_expansion_zeroelim(abtlen, abt, cdheighttail, w);
6465 finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
6467 finswap = finnow; finnow = finother; finother = finswap;
6470 return finnow[finlength - 1];
6473 #ifdef ANSI_DECLARATORS
6475 vertex pa, vertex pb, vertex pc, vertex pd,
6476 REAL aheight, REAL bheight, REAL cheight, REAL dheight)
6478 REAL orient3d(m, b, pa, pb, pc, pd, aheight, bheight, cheight, dheight)
6492 REAL adx, bdx, cdx, ady, bdy, cdy, adheight, bdheight, cdheight;
6493 REAL bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady;
6495 REAL permanent, errbound;
6499 adx = pa[0] - pd[0];
6500 bdx = pb[0] - pd[0];
6501 cdx = pc[0] - pd[0];
6502 ady = pa[1] - pd[1];
6503 bdy = pb[1] - pd[1];
6504 cdy = pc[1] - pd[1];
6505 adheight = aheight - dheight;
6506 bdheight = bheight - dheight;
6507 cdheight = cheight - dheight;
6518 det = adheight * (bdxcdy - cdxbdy)
6519 + bdheight * (cdxady - adxcdy)
6520 + cdheight * (adxbdy - bdxady);
6526 permanent = (Absolute(bdxcdy) + Absolute(cdxbdy)) * Absolute(adheight)
6527 + (Absolute(cdxady) + Absolute(adxcdy)) * Absolute(bdheight)
6528 + (Absolute(adxbdy) + Absolute(bdxady)) * Absolute(cdheight);
6529 errbound = o3derrboundA * permanent;
6530 if ((det > errbound) || (-det > errbound)) {
6534 return orient3dadapt(pa, pb, pc, pd, aheight, bheight, cheight, dheight,
6556 #ifdef ANSI_DECLARATORS
6558 vertex pa, vertex pb, vertex pc, vertex pd)
6560 REAL nonregular(m, b, pa, pb, pc, pd)
6570 if (b->weighted == 0) {
6571 return incircle(m, b, pa, pb, pc, pd);
6572 }
else if (b->weighted == 1) {
6573 return orient3d(m, b, pa, pb, pc, pd,
6574 pa[0] * pa[0] + pa[1] * pa[1] - pa[2],
6575 pb[0] * pb[0] + pb[1] * pb[1] - pb[2],
6576 pc[0] * pc[0] + pc[1] * pc[1] - pc[2],
6577 pd[0] * pd[0] + pd[1] * pd[1] - pd[2]);
6579 return orient3d(m, b, pa, pb, pc, pd, pa[2], pb[2], pc[2], pd[2]);
6597 #ifdef ANSI_DECLARATORS
6598 void findcircumcenter(
struct mesh *m,
struct behavior *b,
6599 vertex torg, vertex tdest, vertex tapex,
6600 vertex circumcenter, REAL *xi, REAL *eta,
int offcenter)
6602 void findcircumcenter(m, b, torg, tdest, tapex, circumcenter, xi, eta,
6609 vertex circumcenter;
6616 REAL xdo, ydo, xao, yao;
6617 REAL dodist, aodist, dadist;
6619 REAL dx, dy, dxoff, dyoff;
6621 m->circumcentercount++;
6624 xdo = tdest[0] - torg[0];
6625 ydo = tdest[1] - torg[1];
6626 xao = tapex[0] - torg[0];
6627 yao = tapex[1] - torg[1];
6628 dodist = xdo * xdo + ydo * ydo;
6629 aodist = xao * xao + yao * yao;
6630 dadist = (tdest[0] - tapex[0]) * (tdest[0] - tapex[0]) +
6631 (tdest[1] - tapex[1]) * (tdest[1] - tapex[1]);
6633 denominator = 0.5 / (xdo * yao - xao * ydo);
6638 denominator = 0.5 / counterclockwise(m, b, tdest, tapex, torg);
6640 m->counterclockcount--;
6642 dx = (yao * dodist - ydo * aodist) * denominator;
6643 dy = (xdo * aodist - xao * dodist) * denominator;
6650 if ((dodist < aodist) && (dodist < dadist)) {
6651 if (offcenter && (b->offconstant > 0.0)) {
6653 dxoff = 0.5 * xdo - b->offconstant * ydo;
6654 dyoff = 0.5 * ydo + b->offconstant * xdo;
6657 if (dxoff * dxoff + dyoff * dyoff < dx * dx + dy * dy) {
6662 }
else if (aodist < dadist) {
6663 if (offcenter && (b->offconstant > 0.0)) {
6664 dxoff = 0.5 * xao + b->offconstant * yao;
6665 dyoff = 0.5 * yao - b->offconstant * xao;
6668 if (dxoff * dxoff + dyoff * dyoff < dx * dx + dy * dy) {
6674 if (offcenter && (b->offconstant > 0.0)) {
6675 dxoff = 0.5 * (tapex[0] - tdest[0]) -
6676 b->offconstant * (tapex[1] - tdest[1]);
6677 dyoff = 0.5 * (tapex[1] - tdest[1]) +
6678 b->offconstant * (tapex[0] - tdest[0]);
6681 if (dxoff * dxoff + dyoff * dyoff <
6682 (dx - xdo) * (dx - xdo) + (dy - ydo) * (dy - ydo)) {
6689 circumcenter[0] = torg[0] + dx;
6690 circumcenter[1] = torg[1] + dy;
6697 *xi = (yao * dx - xao * dy) * (2.0 * denominator);
6698 *eta = (xdo * dy - ydo * dx) * (2.0 * denominator);
6711 #ifdef ANSI_DECLARATORS
6712 void triangleinit(
struct mesh *m)
6714 void triangleinit(m)
6719 poolzero(&m->vertices);
6720 poolzero(&m->triangles);
6721 poolzero(&m->subsegs);
6723 poolzero(&m->badsubsegs);
6724 poolzero(&m->badtriangles);
6725 poolzero(&m->flipstackers);
6726 poolzero(&m->splaynodes);
6728 m->recenttri.tri = (triangle *) NULL;
6731 m->checksegments = 0;
6732 m->checkquality = 0;
6733 m->incirclecount = m->counterclockcount = m->orient3dcount = 0;
6734 m->hyperbolacount = m->circletopcount = m->circumcentercount = 0;
6750 #ifdef ANSI_DECLARATORS
6751 unsigned long randomnation(
unsigned int choices)
6753 unsigned long randomnation(choices)
6758 randomseed = (randomseed * 1366l + 150889l) % 714025l;
6759 return randomseed / (714025l / choices + 1);
6774 #ifdef ANSI_DECLARATORS
6777 void checkmesh(m, b)
6783 struct otri triangleloop;
6784 struct otri oppotri, oppooppotri;
6785 vertex triorg, tridest, triapex;
6786 vertex oppoorg, oppodest;
6792 saveexact = b->noexact;
6795 printf(
" Checking consistency of mesh...\n");
6799 traversalinit(&m->triangles);
6800 triangleloop.tri = triangletraverse(m);
6801 while (triangleloop.tri != (triangle *) NULL) {
6803 for (triangleloop.orient = 0; triangleloop.orient < 3;
6804 triangleloop.orient++) {
6805 org(triangleloop, triorg);
6806 dest(triangleloop, tridest);
6807 if (triangleloop.orient == 0) {
6809 apex(triangleloop, triapex);
6810 if (counterclockwise(m, b, triorg, tridest, triapex) <= 0.0) {
6811 printf(
" !! !! Inverted ");
6812 printtriangle(m, b, &triangleloop);
6817 sym(triangleloop, oppotri);
6818 if (oppotri.tri != m->dummytri) {
6820 sym(oppotri, oppooppotri);
6821 if ((triangleloop.tri != oppooppotri.tri)
6822 || (triangleloop.orient != oppooppotri.orient)) {
6823 printf(
" !! !! Asymmetric triangle-triangle bond:\n");
6824 if (triangleloop.tri == oppooppotri.tri) {
6825 printf(
" (Right triangle, wrong orientation)\n");
6828 printtriangle(m, b, &triangleloop);
6829 printf(
" Second (nonreciprocating) ");
6830 printtriangle(m, b, &oppotri);
6835 org(oppotri, oppoorg);
6836 dest(oppotri, oppodest);
6837 if ((triorg != oppodest) || (tridest != oppoorg)) {
6838 printf(
" !! !! Mismatched edge coordinates between two triangles:\n"
6840 printf(
" First mismatched ");
6841 printtriangle(m, b, &triangleloop);
6842 printf(
" Second mismatched ");
6843 printtriangle(m, b, &oppotri);
6848 triangleloop.tri = triangletraverse(m);
6852 printf(
" In my studied opinion, the mesh appears to be consistent.\n");
6854 }
else if (horrors == 1) {
6855 printf(
" !! !! !! !! Precisely one festering wound discovered.\n");
6857 printf(
" !! !! !! !! %d abominations witnessed.\n", horrors);
6860 b->noexact = saveexact;
6873 #ifdef ANSI_DECLARATORS
6876 void checkdelaunay(m, b)
6882 struct otri triangleloop;
6883 struct otri oppotri;
6884 struct osub opposubseg;
6885 vertex triorg, tridest, triapex;
6887 int shouldbedelaunay;
6894 saveexact = b->noexact;
6897 printf(
" Checking Delaunay property of mesh...\n");
6901 traversalinit(&m->triangles);
6902 triangleloop.tri = triangletraverse(m);
6903 while (triangleloop.tri != (triangle *) NULL) {
6905 for (triangleloop.orient = 0; triangleloop.orient < 3;
6906 triangleloop.orient++) {
6907 org(triangleloop, triorg);
6908 dest(triangleloop, tridest);
6909 apex(triangleloop, triapex);
6910 sym(triangleloop, oppotri);
6911 apex(oppotri, oppoapex);
6915 shouldbedelaunay = (oppotri.tri != m->dummytri) &&
6916 !deadtri(oppotri.tri) && (triangleloop.tri < oppotri.tri) &&
6917 (triorg != m->infvertex1) && (triorg != m->infvertex2) &&
6918 (triorg != m->infvertex3) &&
6919 (tridest != m->infvertex1) && (tridest != m->infvertex2) &&
6920 (tridest != m->infvertex3) &&
6921 (triapex != m->infvertex1) && (triapex != m->infvertex2) &&
6922 (triapex != m->infvertex3) &&
6923 (oppoapex != m->infvertex1) && (oppoapex != m->infvertex2) &&
6924 (oppoapex != m->infvertex3);
6925 if (m->checksegments && shouldbedelaunay) {
6928 tspivot(triangleloop, opposubseg);
6929 if (opposubseg.ss != m->dummysub){
6930 shouldbedelaunay = 0;
6933 if (shouldbedelaunay) {
6934 if (nonregular(m, b, triorg, tridest, triapex, oppoapex) > 0.0) {
6936 printf(
" !! !! Non-Delaunay pair of triangles:\n");
6937 printf(
" First non-Delaunay ");
6938 printtriangle(m, b, &triangleloop);
6939 printf(
" Second non-Delaunay ");
6941 printf(
" !! !! Non-regular pair of triangles:\n");
6942 printf(
" First non-regular ");
6943 printtriangle(m, b, &triangleloop);
6944 printf(
" Second non-regular ");
6946 printtriangle(m, b, &oppotri);
6951 triangleloop.tri = triangletraverse(m);
6956 " By virtue of my perceptive intelligence, I declare the mesh Delaunay.\n");
6958 }
else if (horrors == 1) {
6960 " !! !! !! !! Precisely one terrifying transgression identified.\n");
6962 printf(
" !! !! !! !! %d obscenities viewed with horror.\n", horrors);
6965 b->noexact = saveexact;