Chaste Commit::ffe82dac2b54c9cf0c85e574690a37a97450cb18
triangle.cpp
Go to the documentation of this file.
1
8/*****************************************************************************/
9/* */
10/* 888888888 ,o, / 888 */
11/* 888 88o88o " o8888o 88o8888o o88888o 888 o88888o */
12/* 888 888 888 88b 888 888 888 888 888 d888 88b */
13/* 888 888 888 o88^o888 888 888 "88888" 888 8888oo888 */
14/* 888 888 888 C888 888 888 888 / 888 q888 */
15/* 888 888 888 "88o^888 888 888 Cb 888 "88oooo" */
16/* "8oo8D */
17/* */
18/* A Two-Dimensional Quality Mesh Generator and Delaunay Triangulator. */
19/* (triangle.c) */
20/* */
21/* Version 1.6 */
22/* July 28, 2005 */
23/* */
24/* Copyright 1993, 1995, 1997, 1998, 2002, 2005 */
25/* Jonathan Richard Shewchuk */
26/* 2360 Woolsey #H */
27/* Berkeley, California 94705-1927 */
28/* jrs@cs.berkeley.edu */
29/* */
30/* This program may be freely redistributed under the condition that the */
31/* copyright notices (including this entire header and the copyright */
32/* notice printed when the `-h' switch is selected) are not removed, and */
33/* no compensation is received. Private, research, and institutional */
34/* use is free. You may distribute modified versions of this code UNDER */
35/* THE CONDITION THAT THIS CODE AND ANY MODIFICATIONS MADE TO IT IN THE */
36/* SAME FILE REMAIN UNDER COPYRIGHT OF THE ORIGINAL AUTHOR, BOTH SOURCE */
37/* AND OBJECT CODE ARE MADE FREELY AVAILABLE WITHOUT CHARGE, AND CLEAR */
38/* NOTICE IS GIVEN OF THE MODIFICATIONS. Distribution of this code as */
39/* part of a commercial system is permissible ONLY BY DIRECT ARRANGEMENT */
40/* WITH THE AUTHOR. (If you are not directly supplying this code to a */
41/* customer, and you are instead telling them how they can obtain it for */
42/* free, then you are not required to make any arrangement with me.) */
43/* */
44/* Hypertext instructions for Triangle are available on the Web at */
45/* */
46/* http://www.cs.cmu.edu/~quake/triangle.html */
47/* */
48/* Disclaimer: Neither I nor Carnegie Mellon warrant this code in any way */
49/* whatsoever. This code is provided "as-is". Use at your own risk. */
50/* */
51/* Some of the references listed below are marked with an asterisk. [*] */
52/* These references are available for downloading from the Web page */
53/* */
54/* http://www.cs.cmu.edu/~quake/triangle.research.html */
55/* */
56/* Three papers discussing aspects of Triangle are available. A short */
57/* overview appears in "Triangle: Engineering a 2D Quality Mesh */
58/* Generator and Delaunay Triangulator," in Applied Computational */
59/* Geometry: Towards Geometric Engineering, Ming C. Lin and Dinesh */
60/* Manocha, editors, Lecture Notes in Computer Science volume 1148, */
61/* pages 203-222, Springer-Verlag, Berlin, May 1996 (from the First ACM */
62/* Workshop on Applied Computational Geometry). [*] */
63/* */
64/* The algorithms are discussed in the greatest detail in "Delaunay */
65/* Refinement Algorithms for Triangular Mesh Generation," Computational */
66/* Geometry: Theory and Applications 22(1-3):21-74, May 2002. [*] */
67/* */
68/* More detail about the data structures may be found in my dissertation: */
69/* "Delaunay Refinement Mesh Generation," Ph.D. thesis, Technical Report */
70/* CMU-CS-97-137, School of Computer Science, Carnegie Mellon University, */
71/* Pittsburgh, Pennsylvania, 18 May 1997. [*] */
72/* */
73/* Triangle was created as part of the Quake Project in the School of */
74/* Computer Science at Carnegie Mellon University. For further */
75/* information, see Hesheng Bao, Jacobo Bielak, Omar Ghattas, Loukas F. */
76/* Kallivokas, David R. O'Hallaron, Jonathan R. Shewchuk, and Jifeng Xu, */
77/* "Large-scale Simulation of Elastic Wave Propagation in Heterogeneous */
78/* Media on Parallel Computers," Computer Methods in Applied Mechanics */
79/* and Engineering 152(1-2):85-102, 22 January 1998. */
80/* */
81/* Triangle's Delaunay refinement algorithm for quality mesh generation is */
82/* a hybrid of one due to Jim Ruppert, "A Delaunay Refinement Algorithm */
83/* for Quality 2-Dimensional Mesh Generation," Journal of Algorithms */
84/* 18(3):548-585, May 1995 [*], and one due to L. Paul Chew, "Guaranteed- */
85/* Quality Mesh Generation for Curved Surfaces," Proceedings of the Ninth */
86/* Annual Symposium on Computational Geometry (San Diego, California), */
87/* pages 274-280, Association for Computing Machinery, May 1993, */
88/* http://portal.acm.org/citation.cfm?id=161150 . */
89/* */
90/* The Delaunay refinement algorithm has been modified so that it meshes */
91/* domains with small input angles well, as described in Gary L. Miller, */
92/* Steven E. Pav, and Noel J. Walkington, "When and Why Ruppert's */
93/* Algorithm Works," Twelfth International Meshing Roundtable, pages */
94/* 91-102, Sandia National Laboratories, September 2003. [*] */
95/* */
96/* My implementation of the divide-and-conquer and incremental Delaunay */
97/* triangulation algorithms follows closely the presentation of Guibas */
98/* and Stolfi, even though I use a triangle-based data structure instead */
99/* of their quad-edge data structure. (In fact, I originally implemented */
100/* Triangle using the quad-edge data structure, but the switch to a */
101/* triangle-based data structure sped Triangle by a factor of two.) The */
102/* mesh manipulation primitives and the two aforementioned Delaunay */
103/* triangulation algorithms are described by Leonidas J. Guibas and Jorge */
104/* Stolfi, "Primitives for the Manipulation of General Subdivisions and */
105/* the Computation of Voronoi Diagrams," ACM Transactions on Graphics */
106/* 4(2):74-123, April 1985, http://portal.acm.org/citation.cfm?id=282923 .*/
107/* */
108/* Their O(n log n) divide-and-conquer algorithm is adapted from Der-Tsai */
109/* Lee and Bruce J. Schachter, "Two Algorithms for Constructing the */
110/* Delaunay Triangulation," International Journal of Computer and */
111/* Information Science 9(3):219-242, 1980. Triangle's improvement of the */
112/* divide-and-conquer algorithm by alternating between vertical and */
113/* horizontal cuts was introduced by Rex A. Dwyer, "A Faster Divide-and- */
114/* Conquer Algorithm for Constructing Delaunay Triangulations," */
115/* Algorithmica 2(2):137-151, 1987. */
116/* */
117/* The incremental insertion algorithm was first proposed by C. L. Lawson, */
118/* "Software for C1 Surface Interpolation," in Mathematical Software III, */
119/* John R. Rice, editor, Academic Press, New York, pp. 161-194, 1977. */
120/* For point location, I use the algorithm of Ernst P. Mucke, Isaac */
121/* Saias, and Binhai Zhu, "Fast Randomized Point Location Without */
122/* Preprocessing in Two- and Three-Dimensional Delaunay Triangulations," */
123/* Proceedings of the Twelfth Annual Symposium on Computational Geometry, */
124/* ACM, May 1996. [*] If I were to randomize the order of vertex */
125/* insertion (I currently don't bother), their result combined with the */
126/* result of Kenneth L. Clarkson and Peter W. Shor, "Applications of */
127/* Random Sampling in Computational Geometry II," Discrete & */
128/* Computational Geometry 4(1):387-421, 1989, would yield an expected */
129/* O(n^{4/3}) bound on running time. */
130/* */
131/* The O(n log n) sweepline Delaunay triangulation algorithm is taken from */
132/* Steven Fortune, "A Sweepline Algorithm for Voronoi Diagrams", */
133/* Algorithmica 2(2):153-174, 1987. A random sample of edges on the */
134/* boundary of the triangulation are maintained in a splay tree for the */
135/* purpose of point location. Splay trees are described by Daniel */
136/* Dominic Sleator and Robert Endre Tarjan, "Self-Adjusting Binary Search */
137/* Trees," Journal of the ACM 32(3):652-686, July 1985, */
138/* http://portal.acm.org/citation.cfm?id=3835 . */
139/* */
140/* The algorithms for exact computation of the signs of determinants are */
141/* described in Jonathan Richard Shewchuk, "Adaptive Precision Floating- */
142/* Point Arithmetic and Fast Robust Geometric Predicates," Discrete & */
143/* Computational Geometry 18(3):305-363, October 1997. (Also available */
144/* as Technical Report CMU-CS-96-140, School of Computer Science, */
145/* Carnegie Mellon University, Pittsburgh, Pennsylvania, May 1996.) [*] */
146/* An abbreviated version appears as Jonathan Richard Shewchuk, "Robust */
147/* Adaptive Floating-Point Geometric Predicates," Proceedings of the */
148/* Twelfth Annual Symposium on Computational Geometry, ACM, May 1996. [*] */
149/* Many of the ideas for my exact arithmetic routines originate with */
150/* Douglas M. Priest, "Algorithms for Arbitrary Precision Floating Point */
151/* Arithmetic," Tenth Symposium on Computer Arithmetic, pp. 132-143, IEEE */
152/* Computer Society Press, 1991. [*] Many of the ideas for the correct */
153/* evaluation of the signs of determinants are taken from Steven Fortune */
154/* and Christopher J. Van Wyk, "Efficient Exact Arithmetic for Computa- */
155/* tional Geometry," Proceedings of the Ninth Annual Symposium on */
156/* Computational Geometry, ACM, pp. 163-172, May 1993, and from Steven */
157/* Fortune, "Numerical Stability of Algorithms for 2D Delaunay Triangu- */
158/* lations," International Journal of Computational Geometry & Applica- */
159/* tions 5(1-2):193-213, March-June 1995. */
160/* */
161/* The method of inserting new vertices off-center (not precisely at the */
162/* circumcenter of every poor-quality triangle) is from Alper Ungor, */
163/* "Off-centers: A New Type of Steiner Points for Computing Size-Optimal */
164/* Quality-Guaranteed Delaunay Triangulations," Proceedings of LATIN */
165/* 2004 (Buenos Aires, Argentina), April 2004. */
166/* */
167/* For definitions of and results involving Delaunay triangulations, */
168/* constrained and conforming versions thereof, and other aspects of */
169/* triangular mesh generation, see the excellent survey by Marshall Bern */
170/* and David Eppstein, "Mesh Generation and Optimal Triangulation," in */
171/* Computing and Euclidean Geometry, Ding-Zhu Du and Frank Hwang, */
172/* editors, World Scientific, Singapore, pp. 23-90, 1992. [*] */
173/* */
174/* The time for incrementally adding PSLG (planar straight line graph) */
175/* segments to create a constrained Delaunay triangulation is probably */
176/* O(t^2) per segment in the worst case and O(t) per segment in the */
177/* common case, where t is the number of triangles that intersect the */
178/* segment before it is inserted. This doesn't count point location, */
179/* which can be much more expensive. I could improve this to O(d log d) */
180/* time, but d is usually quite small, so it's not worth the bother. */
181/* (This note does not apply when the -s switch is used, invoking a */
182/* different method is used to insert segments.) */
183/* */
184/* The time for deleting a vertex from a Delaunay triangulation is O(d^2) */
185/* in the worst case and O(d) in the common case, where d is the degree */
186/* of the vertex being deleted. I could improve this to O(d log d) time, */
187/* but d is usually quite small, so it's not worth the bother. */
188/* */
189/* Ruppert's Delaunay refinement algorithm typically generates triangles */
190/* at a linear rate (constant time per triangle) after the initial */
191/* triangulation is formed. There may be pathological cases where */
192/* quadratic time is required, but these never arise in practice. */
193/* */
194/* The geometric predicates (circumcenter calculations, segment */
195/* intersection formulae, etc.) appear in my "Lecture Notes on Geometric */
196/* Robustness" at http://www.cs.berkeley.edu/~jrs/mesh . */
197/* */
198/* If you make any improvements to this code, please please please let me */
199/* know, so that I may obtain the improvements. Even if you don't change */
200/* the code, I'd still love to hear what it's being used for. */
201/* */
202/*****************************************************************************/
203
204/* For single precision (which will save some memory and reduce paging), */
205/* define the symbol SINGLE by using the -DSINGLE compiler switch or by */
206/* writing "#define SINGLE" below. */
207/* */
208/* For double precision (which will allow you to refine meshes to a smaller */
209/* edge length), leave SINGLE undefined. */
210/* */
211/* Double precision uses more memory, but improves the resolution of the */
212/* meshes you can generate with Triangle. It also reduces the likelihood */
213/* of a floating exception due to overflow. Finally, it is much faster */
214/* than single precision on 64-bit architectures like the DEC Alpha. I */
215/* recommend double precision unless you want to generate a mesh for which */
216/* you do not have enough memory. */
217
218/* #define SINGLE */
219
220#ifdef SINGLE
221#define REAL float
222#else /* not SINGLE */
223#define REAL double
224#endif /* not SINGLE */
225
226/* If yours is not a Unix system, define the NO_TIMER compiler switch to */
227/* remove the Unix-specific timing code. */
228
229/* #define NO_TIMER */
230
231/* To insert lots of self-checks for internal errors, define the SELF_CHECK */
232/* symbol. This will slow down the program significantly. It is best to */
233/* define the symbol using the -DSELF_CHECK compiler switch, but you could */
234/* write "#define SELF_CHECK" below. If you are modifying this code, I */
235/* recommend you turn self-checks on until your work is debugged. */
236
237/* #define SELF_CHECK */
238
239/* To compile Triangle as a callable object library (triangle.o), define the */
240/* TRILIBRARY symbol. Read the file triangle.h for details on how to call */
241/* the procedure triangulate() that results. */
242
243/* #define TRILIBRARY */
244
245/* It is possible to generate a smaller version of Triangle using one or */
246/* both of the following symbols. Define the REDUCED symbol to eliminate */
247/* all features that are primarily of research interest; specifically, the */
248/* -i, -F, -s, and -C switches. Define the CDT_ONLY symbol to eliminate */
249/* all meshing algorithms above and beyond constrained Delaunay */
250/* triangulation; specifically, the -r, -q, -a, -u, -D, -S, and -s */
251/* switches. These reductions are most likely to be useful when */
252/* generating an object library (triangle.o) by defining the TRILIBRARY */
253/* symbol. */
254
255/* #define REDUCED */
256/* #define CDT_ONLY */
257
258/* On some machines, my exact arithmetic routines might be defeated by the */
259/* use of internal extended precision floating-point registers. The best */
260/* way to solve this problem is to set the floating-point registers to use */
261/* single or double precision internally. On 80x86 processors, this may */
262/* be accomplished by setting the CPU86 symbol for the Microsoft C */
263/* compiler, or the LINUX symbol for the gcc compiler running on Linux. */
264/* */
265/* An inferior solution is to declare certain values as `volatile', thus */
266/* forcing them to be stored to memory and rounded off. Unfortunately, */
267/* this solution might slow Triangle down quite a bit. To use volatile */
268/* values, write "#define INEXACT volatile" below. Normally, however, */
269/* INEXACT should be defined to be nothing. ("#define INEXACT".) */
270/* */
271/* For more discussion, see http://www.cs.cmu.edu/~quake/robust.pc.html . */
272/* For yet more discussion, see Section 5 of my paper, "Adaptive Precision */
273/* Floating-Point Arithmetic and Fast Robust Geometric Predicates" (also */
274/* available as Section 6.6 of my dissertation). */
275
276/* #define CPU86 */
277/* #define LINUX */
278
279#define INEXACT /* Nothing */
280/* #define INEXACT volatile */
281
282/* Maximum number of characters in a file name (including the null). */
283
284#define FILENAMESIZE 2048
285
286/* Maximum number of characters in a line read from a file (including the */
287/* null). */
288
289#define INPUTLINESIZE 1024
290
291/* For efficiency, a variety of data structures are allocated in bulk. The */
292/* following constants determine how many of each structure is allocated */
293/* at once. */
294
295#define TRIPERBLOCK 4092 /* Number of triangles allocated at once. */
296#define SUBSEGPERBLOCK 508 /* Number of subsegments allocated at once. */
297#define VERTEXPERBLOCK 4092 /* Number of vertices allocated at once. */
298#define VIRUSPERBLOCK 1020 /* Number of virus triangles allocated at once. */
299/* Number of encroached subsegments allocated at once. */
300#define BADSUBSEGPERBLOCK 252
301/* Number of skinny triangles allocated at once. */
302#define BADTRIPERBLOCK 4092
303/* Number of flipped triangles allocated at once. */
304#define FLIPSTACKERPERBLOCK 252
305/* Number of splay tree nodes allocated at once. */
306#define SPLAYNODEPERBLOCK 508
307
308/* The vertex types. A DEADVERTEX has been deleted entirely. An */
309/* UNDEADVERTEX is not part of the mesh, but is written to the output */
310/* .node file and affects the node indexing in the other output files. */
311
312#define INPUTVERTEX 0
313#define SEGMENTVERTEX 1
314#define FREEVERTEX 2
315#define DEADVERTEX -32768
316#define UNDEADVERTEX -32767
317
318/* The next line is used to outsmart some very stupid compilers. If your */
319/* compiler is smarter, feel free to replace the "int" with "void". */
320/* Not that it matters. */
321
322/*Windows_Port_Begins*/
323#ifdef _MSC_VER
324#else
325 /*This gave me a whole lot of grief with Windows SDK*/
326 #define VOID int
327#endif
328/*Windows_Port_Begins*/
329
330/* Two constants for algorithms based on random sampling. Both constants */
331/* have been chosen empirically to optimize their respective algorithms. */
332
333/* Used for the point location scheme of Mucke, Saias, and Zhu, to decide */
334/* how large a random sample of triangles to inspect. */
335
336#define SAMPLEFACTOR 11
337
338/* Used in Fortune's sweepline Delaunay algorithm to determine what fraction */
339/* of boundary edges should be maintained in the splay tree for point */
340/* location on the front. */
341
342#define SAMPLERATE 10
343
344/* A number that speaks for itself, every kissable digit. */
345
346#define PI 3.141592653589793238462643383279502884197169399375105820974944592308
347
348/* Another fave. */
349
350#define SQUAREROOTTWO 1.4142135623730950488016887242096980785696718753769480732
351
352/* And here's one for those of you who are intimidated by math. */
353
354#define ONETHIRD 0.333333333333333333333333333333333333333333333333333333333333
355
356#include <stdio.h>
357#include <stdlib.h>
358#include <string.h>
359#include <math.h>
360#ifndef NO_TIMER
361#ifdef _MSC_VER
362#include <time.h>
363#include <WinSock2.h>
364#include <windows.h>
365//The portion below is adapted from http://www.suacommunity.com/dictionary/gettimeofday-entry.php
366#if defined(_MSC_EXTENSIONS)
367 #define DELTA_EPOCH_IN_MICROSECS 11644473600000000Ui64
368#else
369 #define DELTA_EPOCH_IN_MICROSECS 11644473600000000ULL
370#endif
371
372struct timezone
373{
374 int tz_minuteswest; /* minutes W of Greenwich */
375 int tz_dsttime; /* type of dst correction */
376};
377
378// Definition of a gettimeofday function
379
380int gettimeofday(struct timeval *tv, struct timezone *tz)
381{
382 // Define a structure to receive the current Windows filetime
383 FILETIME ft;
384
385 // Initialize the present time to 0 and the timezone to UTC
386 unsigned __int64 tmpres = 0;
387 static int tzflag = 0;
388
389 if (NULL != tv)
390 {
391 GetSystemTimeAsFileTime(&ft);
392
393 // The GetSystemTimeAsFileTime returns the number of 100 nanosecond
394 // intervals since Jan 1, 1601 in a structure. Copy the high bits to
395 // the 64 bit tmpres, shift it left by 32 then or in the low 32 bits.
396 tmpres |= ft.dwHighDateTime;
397 tmpres <<= 32;
398 tmpres |= ft.dwLowDateTime;
399
400 // Convert to microseconds by dividing by 10
401 tmpres /= 10;
402
403 // The Unix epoch starts on Jan 1 1970. Need to subtract the difference
404 // in seconds from Jan 1 1601.
405 tmpres -= DELTA_EPOCH_IN_MICROSECS;
406
407 // Finally change microseconds to seconds and place in the seconds value.
408 // The modulus picks up the microseconds.
409 tv->tv_sec = (long)(tmpres / 1000000UL);
410 tv->tv_usec = (long)(tmpres % 1000000UL);
411 }
412
413 if (NULL != tz)
414 {
415 if (!tzflag)
416 {
417 _tzset();
418 tzflag++;
419 }
420
421 // Adjust for the timezone west of Greenwich
422 tz->tz_minuteswest = _timezone / 60;
423 tz->tz_dsttime = _daylight;
424 }
425 return 0;
426}
427#else
428#include <sys/time.h>
429#endif
430#endif /* not NO_TIMER */
431#ifdef CPU86
432#include <float.h>
433#endif /* CPU86 */
434#ifdef LINUX
435#include <fpu_control.h>
436#endif /* LINUX */
437#ifdef TRILIBRARY
438#include "triangle.h"
439#endif /* TRILIBRARY */
440
441/* A few forward declarations. */
442
443#ifndef TRILIBRARY
444char *readline();
445char *findfield();
446#endif /* not TRILIBRARY */
447
448/* Labels that signify the result of point location. The result of a */
449/* search indicates that the point falls in the interior of a triangle, on */
450/* an edge, on a vertex, or outside the mesh. */
451
452enum locateresult {INTRIANGLE, ONEDGE, ONVERTEX, OUTSIDE};
453
454/* Labels that signify the result of vertex insertion. The result indicates */
455/* that the vertex was inserted with complete success, was inserted but */
456/* encroaches upon a subsegment, was not inserted because it lies on a */
457/* segment, or was not inserted because another vertex occupies the same */
458/* location. */
459
460enum insertvertexresult {SUCCESSFULVERTEX, ENCROACHINGVERTEX, VIOLATINGVERTEX,
461 DUPLICATEVERTEX};
462
463/* Labels that signify the result of direction finding. The result */
464/* indicates that a segment connecting the two query points falls within */
465/* the direction triangle, along the left edge of the direction triangle, */
466/* or along the right edge of the direction triangle. */
467
468enum finddirectionresult {WITHIN, LEFTCOLLINEAR, RIGHTCOLLINEAR};
469
470/*****************************************************************************/
471/* */
472/* The basic mesh data structures */
473/* */
474/* There are three: vertices, triangles, and subsegments (abbreviated */
475/* `subseg'). These three data structures, linked by pointers, comprise */
476/* the mesh. A vertex simply represents a mesh vertex and its properties. */
477/* A triangle is a triangle. A subsegment is a special data structure used */
478/* to represent an impenetrable edge of the mesh (perhaps on the outer */
479/* boundary, on the boundary of a hole, or part of an internal boundary */
480/* separating two triangulated regions). Subsegments represent boundaries, */
481/* defined by the user, that triangles may not lie across. */
482/* */
483/* A triangle consists of a list of three vertices, a list of three */
484/* adjoining triangles, a list of three adjoining subsegments (when */
485/* segments exist), an arbitrary number of optional user-defined */
486/* floating-point attributes, and an optional area constraint. The latter */
487/* is an upper bound on the permissible area of each triangle in a region, */
488/* used for mesh refinement. */
489/* */
490/* For a triangle on a boundary of the mesh, some or all of the neighboring */
491/* triangles may not be present. For a triangle in the interior of the */
492/* mesh, often no neighboring subsegments are present. Such absent */
493/* triangles and subsegments are never represented by NULL pointers; they */
494/* are represented by two special records: `dummytri', the triangle that */
495/* fills "outer space", and `dummysub', the omnipresent subsegment. */
496/* `dummytri' and `dummysub' are used for several reasons; for instance, */
497/* they can be dereferenced and their contents examined without violating */
498/* protected memory. */
499/* */
500/* However, it is important to understand that a triangle includes other */
501/* information as well. The pointers to adjoining vertices, triangles, and */
502/* subsegments are ordered in a way that indicates their geometric relation */
503/* to each other. Furthermore, each of these pointers contains orientation */
504/* information. Each pointer to an adjoining triangle indicates which face */
505/* of that triangle is contacted. Similarly, each pointer to an adjoining */
506/* subsegment indicates which side of that subsegment is contacted, and how */
507/* the subsegment is oriented relative to the triangle. */
508/* */
509/* The data structure representing a subsegment may be thought to be */
510/* abutting the edge of one or two triangle data structures: either */
511/* sandwiched between two triangles, or resting against one triangle on an */
512/* exterior boundary or hole boundary. */
513/* */
514/* A subsegment consists of a list of four vertices--the vertices of the */
515/* subsegment, and the vertices of the segment it is a part of--a list of */
516/* two adjoining subsegments, and a list of two adjoining triangles. One */
517/* of the two adjoining triangles may not be present (though there should */
518/* always be one), and neighboring subsegments might not be present. */
519/* Subsegments also store a user-defined integer "boundary marker". */
520/* Typically, this integer is used to indicate what boundary conditions are */
521/* to be applied at that location in a finite element simulation. */
522/* */
523/* Like triangles, subsegments maintain information about the relative */
524/* orientation of neighboring objects. */
525/* */
526/* Vertices are relatively simple. A vertex is a list of floating-point */
527/* numbers, starting with the x, and y coordinates, followed by an */
528/* arbitrary number of optional user-defined floating-point attributes, */
529/* followed by an integer boundary marker. During the segment insertion */
530/* phase, there is also a pointer from each vertex to a triangle that may */
531/* contain it. Each pointer is not always correct, but when one is, it */
532/* speeds up segment insertion. These pointers are assigned values once */
533/* at the beginning of the segment insertion phase, and are not used or */
534/* updated except during this phase. Edge flipping during segment */
535/* insertion will render some of them incorrect. Hence, don't rely upon */
536/* them for anything. */
537/* */
538/* Other than the exception mentioned above, vertices have no information */
539/* about what triangles, subfacets, or subsegments they are linked to. */
540/* */
541/*****************************************************************************/
542
543/*****************************************************************************/
544/* */
545/* Handles */
546/* */
547/* The oriented triangle (`otri') and oriented subsegment (`osub') data */
548/* structures defined below do not themselves store any part of the mesh. */
549/* The mesh itself is made of `triangle's, `subseg's, and `vertex's. */
550/* */
551/* Oriented triangles and oriented subsegments will usually be referred to */
552/* as "handles." A handle is essentially a pointer into the mesh; it */
553/* allows you to "hold" one particular part of the mesh. Handles are used */
554/* to specify the regions in which one is traversing and modifying the mesh.*/
555/* A single `triangle' may be held by many handles, or none at all. (The */
556/* latter case is not a memory leak, because the triangle is still */
557/* connected to other triangles in the mesh.) */
558/* */
559/* An `otri' is a handle that holds a triangle. It holds a specific edge */
560/* of the triangle. An `osub' is a handle that holds a subsegment. It */
561/* holds either the left or right side of the subsegment. */
562/* */
563/* Navigation about the mesh is accomplished through a set of mesh */
564/* manipulation primitives, further below. Many of these primitives take */
565/* a handle and produce a new handle that holds the mesh near the first */
566/* handle. Other primitives take two handles and glue the corresponding */
567/* parts of the mesh together. The orientation of the handles is */
568/* important. For instance, when two triangles are glued together by the */
569/* bond() primitive, they are glued at the edges on which the handles lie. */
570/* */
571/* Because vertices have no information about which triangles they are */
572/* attached to, I commonly represent a vertex by use of a handle whose */
573/* origin is the vertex. A single handle can simultaneously represent a */
574/* triangle, an edge, and a vertex. */
575/* */
576/*****************************************************************************/
577
578/* The triangle data structure. Each triangle contains three pointers to */
579/* adjoining triangles, plus three pointers to vertices, plus three */
580/* pointers to subsegments (declared below; these pointers are usually */
581/* `dummysub'). It may or may not also contain user-defined attributes */
582/* and/or a floating-point "area constraint." It may also contain extra */
583/* pointers for nodes, when the user asks for high-order elements. */
584/* Because the size and structure of a `triangle' is not decided until */
585/* runtime, I haven't simply declared the type `triangle' as a struct. */
586
587typedef REAL **triangle; /* Really: typedef triangle *triangle */
588
589/* An oriented triangle: includes a pointer to a triangle and orientation. */
590/* The orientation denotes an edge of the triangle. Hence, there are */
591/* three possible orientations. By convention, each edge always points */
592/* counterclockwise about the corresponding triangle. */
593
594struct otri {
595 triangle *tri;
596 int orient; /* Ranges from 0 to 2. */
597};
598
599/* The subsegment data structure. Each subsegment contains two pointers to */
600/* adjoining subsegments, plus four pointers to vertices, plus two */
601/* pointers to adjoining triangles, plus one boundary marker, plus one */
602/* segment number. */
603
604typedef REAL **subseg; /* Really: typedef subseg *subseg */
605
606/* An oriented subsegment: includes a pointer to a subsegment and an */
607/* orientation. The orientation denotes a side of the edge. Hence, there */
608/* are two possible orientations. By convention, the edge is always */
609/* directed so that the "side" denoted is the right side of the edge. */
610
611struct osub {
612 subseg *ss;
613 int ssorient; /* Ranges from 0 to 1. */
614};
615
616/* The vertex data structure. Each vertex is actually an array of REALs. */
617/* The number of REALs is unknown until runtime. An integer boundary */
618/* marker, and sometimes a pointer to a triangle, is appended after the */
619/* REALs. */
620
621typedef REAL *vertex;
622
623/* A queue used to store encroached subsegments. Each subsegment's vertices */
624/* are stored so that we can check whether a subsegment is still the same. */
625
626struct badsubseg {
627 subseg encsubseg; /* An encroached subsegment. */
628 vertex subsegorg, subsegdest; /* Its two vertices. */
629};
630
631/* A queue used to store bad triangles. The key is the square of the cosine */
632/* of the smallest angle of the triangle. Each triangle's vertices are */
633/* stored so that one can check whether a triangle is still the same. */
634
635struct badtriang {
636 triangle poortri; /* A skinny or too-large triangle. */
637 REAL key; /* cos^2 of smallest (apical) angle. */
638 vertex triangorg, triangdest, triangapex; /* Its three vertices. */
639 struct badtriang *nexttriang; /* Pointer to next bad triangle. */
640};
641
642/* A stack of triangles flipped during the most recent vertex insertion. */
643/* The stack is used to undo the vertex insertion if the vertex encroaches */
644/* upon a subsegment. */
645
647 triangle flippedtri; /* A recently flipped triangle. */
648 struct flipstacker* prevflip; /* Previous flip in the stack. */
649};
650
651/* A node in a heap used to store events for the sweepline Delaunay */
652/* algorithm. Nodes do not point directly to their parents or children in */
653/* the heap. Instead, each node knows its position in the heap, and can */
654/* look up its parent and children in a separate array. The `eventptr' */
655/* points either to a `vertex' or to a triangle (in encoded format, so */
656/* that an orientation is included). In the latter case, the origin of */
657/* the oriented triangle is the apex of a "circle event" of the sweepline */
658/* algorithm. To distinguish site events from circle events, all circle */
659/* events are given an invalid (smaller than `xmin') x-coordinate `xkey'. */
660
661struct event {
662 REAL xkey, ykey; /* Coordinates of the event. */
663 VOID *eventptr; /* Can be a vertex or the location of a circle event. */
664 int heapposition; /* Marks this event's position in the heap. */
665};
666
667/* A node in the splay tree. Each node holds an oriented ghost triangle */
668/* that represents a boundary edge of the growing triangulation. When a */
669/* circle event covers two boundary edges with a triangle, so that they */
670/* are no longer boundary edges, those edges are not immediately deleted */
671/* from the tree; rather, they are lazily deleted when they are next */
672/* encountered. (Since only a random sample of boundary edges are kept */
673/* in the tree, lazy deletion is faster.) `keydest' is used to verify */
674/* that a triangle is still the same as when it entered the splay tree; if */
675/* it has been rotated (due to a circle event), it no longer represents a */
676/* boundary edge and should be deleted. */
677
678struct splaynode {
679 struct otri keyedge; /* Lprev of an edge on the front. */
680 vertex keydest; /* Used to verify that splay node is still live. */
681 struct splaynode *lchild, *rchild; /* Children in splay tree. */
682};
683
684/* A type used to allocate memory. firstblock is the first block of items. */
685/* nowblock is the block from which items are currently being allocated. */
686/* nextitem points to the next slab of free memory for an item. */
687/* deaditemstack is the head of a linked list (stack) of deallocated items */
688/* that can be recycled. unallocateditems is the number of items that */
689/* remain to be allocated from nowblock. */
690/* */
691/* Traversal is the process of walking through the entire list of items, and */
692/* is separate from allocation. Note that a traversal will visit items on */
693/* the "deaditemstack" stack as well as live items. pathblock points to */
694/* the block currently being traversed. pathitem points to the next item */
695/* to be traversed. pathitemsleft is the number of items that remain to */
696/* be traversed in pathblock. */
697/* */
698/* alignbytes determines how new records should be aligned in memory. */
699/* itembytes is the length of a record in bytes (after rounding up). */
700/* itemsperblock is the number of items allocated at once in a single */
701/* block. itemsfirstblock is the number of items in the first block, */
702/* which can vary from the others. items is the number of currently */
703/* allocated items. maxitems is the maximum number of items that have */
704/* been allocated at once; it is the current number of items plus the */
705/* number of records kept on deaditemstack. */
706
708 VOID **firstblock, **nowblock;
709 VOID *nextitem;
710 VOID *deaditemstack;
711 VOID **pathblock;
712 VOID* pathitem;
713 int alignbytes;
714 int itembytes;
715 int itemsperblock;
716 int itemsfirstblock;
717 long items, maxitems;
718 int unallocateditems;
719 int pathitemsleft;
720};
721
722
723/* Global constants. */
724
725REAL splitter; /* Used to split REAL factors for exact multiplication. */
726REAL epsilon; /* Floating-point machine epsilon. */
727REAL resulterrbound;
728REAL ccwerrboundA, ccwerrboundB, ccwerrboundC;
729REAL iccerrboundA, iccerrboundB, iccerrboundC;
730REAL o3derrboundA, o3derrboundB, o3derrboundC;
731
732/* Random number seed is not constant, but I've made it global anyway. */
733
734unsigned long randomseed; /* Current random number seed. */
735
736
737/* Mesh data structure. Triangle operates on only one mesh, but the mesh */
738/* structure is used (instead of global variables) to allow reentrancy. */
739
740struct mesh {
741
742/* Variables used to allocate memory for triangles, subsegments, vertices, */
743/* viri (triangles being eaten), encroached segments, bad (skinny or too */
744/* large) triangles, and splay tree nodes. */
745
746 struct memorypool triangles;
747 struct memorypool subsegs;
748 struct memorypool vertices;
749 struct memorypool viri;
750 struct memorypool badsubsegs;
751 struct memorypool badtriangles;
752 struct memorypool flipstackers;
753 struct memorypool splaynodes;
754
755/* Variables that maintain the bad triangle queues. The queues are */
756/* ordered from 4095 (highest priority) to 0 (lowest priority). */
757
758 struct badtriang *queuefront[4096];
759 struct badtriang *queuetail[4096];
760 int nextnonemptyq[4096];
761 int firstnonemptyq;
762
763/* Variable that maintains the stack of recently flipped triangles. */
764
765 struct flipstacker *lastflip;
766
767/* Other variables. */
768
769 REAL xmin, xmax, ymin, ymax; /* x and y bounds. */
770 REAL xminextreme; /* Nonexistent x value used as a flag in sweepline. */
771 int invertices; /* Number of input vertices. */
772 int inelements; /* Number of input triangles. */
773 int insegments; /* Number of input segments. */
774 int holes; /* Number of input holes. */
775 int regions; /* Number of input regions. */
776 int undeads; /* Number of input vertices that don't appear in the mesh. */
777 long edges; /* Number of output edges. */
778 int mesh_dim; /* Dimension (ought to be 2). */
779 int nextras; /* Number of attributes per vertex. */
780 int eextras; /* Number of attributes per triangle. */
781 long hullsize; /* Number of edges in convex hull. */
782 int steinerleft; /* Number of Steiner points not yet used. */
783 int vertexmarkindex; /* Index to find boundary marker of a vertex. */
784 int vertex2triindex; /* Index to find a triangle adjacent to a vertex. */
785 int highorderindex; /* Index to find extra nodes for high-order elements. */
786 int elemattribindex; /* Index to find attributes of a triangle. */
787 int areaboundindex; /* Index to find area bound of a triangle. */
788 int checksegments; /* Are there segments in the triangulation yet? */
789 int checkquality; /* Has quality triangulation begun yet? */
790 int readnodefile; /* Has a .node file been read? */
791 long samples; /* Number of random samples for point location. */
792
793 long incirclecount; /* Number of incircle tests performed. */
794 long counterclockcount; /* Number of counterclockwise tests performed. */
795 long orient3dcount; /* Number of 3D orientation tests performed. */
796 long hyperbolacount; /* Number of right-of-hyperbola tests performed. */
797 long circumcentercount; /* Number of circumcenter calculations performed. */
798 long circletopcount; /* Number of circle top calculations performed. */
799
800/* Triangular bounding box vertices. */
801
802 vertex infvertex1, infvertex2, infvertex3;
803
804/* Pointer to the `triangle' that occupies all of "outer space." */
805
806 triangle *dummytri;
807 triangle *dummytribase; /* Keep base address so we can free() it later. */
808
809/* Pointer to the omnipresent subsegment. Referenced by any triangle or */
810/* subsegment that isn't really connected to a subsegment at that */
811/* location. */
812
813 subseg *dummysub;
814 subseg *dummysubbase; /* Keep base address so we can free() it later. */
815
816/* Pointer to a recently visited triangle. Improves point location if */
817/* proximate vertices are inserted sequentially. */
818
819 struct otri recenttri;
820
821}; /* End of `struct mesh'. */
822
823
824/* Data structure for command line switches and file names. This structure */
825/* is used (instead of global variables) to allow reentrancy. */
826
827struct behavior {
828
829/* Switches for the triangulator. */
830/* poly: -p switch. refine: -r switch. */
831/* quality: -q switch. */
832/* minangle: minimum angle bound, specified after -q switch. */
833/* goodangle: cosine squared of minangle. */
834/* offconstant: constant used to place off-center Steiner points. */
835/* vararea: -a switch without number. */
836/* fixedarea: -a switch with number. */
837/* maxarea: maximum area bound, specified after -a switch. */
838/* usertest: -u switch. */
839/* regionattrib: -A switch. convex: -c switch. */
840/* weighted: 1 for -w switch, 2 for -W switch. jettison: -j switch */
841/* firstnumber: inverse of -z switch. All items are numbered starting */
842/* from `firstnumber'. */
843/* edgesout: -e switch. voronoi: -v switch. */
844/* neighbors: -n switch. geomview: -g switch. */
845/* nobound: -B switch. nopolywritten: -P switch. */
846/* nonodewritten: -N switch. noelewritten: -E switch. */
847/* noiterationnum: -I switch. noholes: -O switch. */
848/* noexact: -X switch. */
849/* order: element order, specified after -o switch. */
850/* nobisect: count of how often -Y switch is selected. */
851/* steiner: maximum number of Steiner points, specified after -S switch. */
852/* incremental: -i switch. sweepline: -F switch. */
853/* dwyer: inverse of -l switch. */
854/* splitseg: -s switch. */
855/* conformdel: -D switch. docheck: -C switch. */
856/* quiet: -Q switch. verbose: count of how often -V switch is selected. */
857/* usesegments: -p, -r, -q, or -c switch; determines whether segments are */
858/* used at all. */
859/* */
860/* Read the instructions to find out the meaning of these switches. */
861
862 int poly, refine, quality, vararea, fixedarea, usertest;
863 int regionattrib, convex, weighted, jettison;
864 int firstnumber;
865 int edgesout, voronoi, neighbors, geomview;
866 int nobound, nopolywritten, nonodewritten, noelewritten, noiterationnum;
867 int noholes, noexact, conformdel;
868 int incremental, sweepline, dwyer;
869 int splitseg;
870 int docheck;
871 int quiet, verbose;
872 int usesegments;
873 int order;
874 int nobisect;
875 int steiner;
876 REAL minangle, goodangle, offconstant;
877 REAL maxarea;
878
879/* Variables for file names. */
880
881#ifndef TRILIBRARY
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];
894#endif /* not TRILIBRARY */
895
896}; /* End of `struct behavior'. */
897
898
899/*****************************************************************************/
900/* */
901/* Mesh manipulation primitives. Each triangle contains three pointers to */
902/* other triangles, with orientations. Each pointer points not to the */
903/* first byte of a triangle, but to one of the first three bytes of a */
904/* triangle. It is necessary to extract both the triangle itself and the */
905/* orientation. To save memory, I keep both pieces of information in one */
906/* pointer. To make this possible, I assume that all triangles are aligned */
907/* to four-byte boundaries. The decode() routine below decodes a pointer, */
908/* extracting an orientation (in the range 0 to 2) and a pointer to the */
909/* beginning of a triangle. The encode() routine compresses a pointer to a */
910/* triangle and an orientation into a single pointer. My assumptions that */
911/* triangles are four-byte-aligned and that the `unsigned long' type is */
912/* long enough to hold a pointer are two of the few kludges in this program.*/
913/* */
914/* Subsegments are manipulated similarly. A pointer to a subsegment */
915/* carries both an address and an orientation in the range 0 to 1. */
916/* */
917/* The other primitives take an oriented triangle or oriented subsegment, */
918/* and return an oriented triangle or oriented subsegment or vertex; or */
919/* they change the connections in the data structure. */
920/* */
921/* Below, triangles and subsegments are denoted by their vertices. The */
922/* triangle abc has origin (org) a, destination (dest) b, and apex (apex) */
923/* c. These vertices occur in counterclockwise order about the triangle. */
924/* The handle abc may simultaneously denote vertex a, edge ab, and triangle */
925/* abc. */
926/* */
927/* Similarly, the subsegment ab has origin (sorg) a and destination (sdest) */
928/* b. If ab is thought to be directed upward (with b directly above a), */
929/* then the handle ab is thought to grasp the right side of ab, and may */
930/* simultaneously denote vertex a and edge ab. */
931/* */
932/* An asterisk (*) denotes a vertex whose identity is unknown. */
933/* */
934/* Given this notation, a partial list of mesh manipulation primitives */
935/* follows. */
936/* */
937/* */
938/* For triangles: */
939/* */
940/* sym: Find the abutting triangle; same edge. */
941/* sym(abc) -> ba* */
942/* */
943/* lnext: Find the next edge (counterclockwise) of a triangle. */
944/* lnext(abc) -> bca */
945/* */
946/* lprev: Find the previous edge (clockwise) of a triangle. */
947/* lprev(abc) -> cab */
948/* */
949/* onext: Find the next edge counterclockwise with the same origin. */
950/* onext(abc) -> ac* */
951/* */
952/* oprev: Find the next edge clockwise with the same origin. */
953/* oprev(abc) -> a*b */
954/* */
955/* dnext: Find the next edge counterclockwise with the same destination. */
956/* dnext(abc) -> *ba */
957/* */
958/* dprev: Find the next edge clockwise with the same destination. */
959/* dprev(abc) -> cb* */
960/* */
961/* rnext: Find the next edge (counterclockwise) of the adjacent triangle. */
962/* rnext(abc) -> *a* */
963/* */
964/* rprev: Find the previous edge (clockwise) of the adjacent triangle. */
965/* rprev(abc) -> b** */
966/* */
967/* org: Origin dest: Destination apex: Apex */
968/* org(abc) -> a dest(abc) -> b apex(abc) -> c */
969/* */
970/* bond: Bond two triangles together at the resepective handles. */
971/* bond(abc, bad) */
972/* */
973/* */
974/* For subsegments: */
975/* */
976/* ssym: Reverse the orientation of a subsegment. */
977/* ssym(ab) -> ba */
978/* */
979/* spivot: Find adjoining subsegment with the same origin. */
980/* spivot(ab) -> a* */
981/* */
982/* snext: Find next subsegment in sequence. */
983/* snext(ab) -> b* */
984/* */
985/* sorg: Origin sdest: Destination */
986/* sorg(ab) -> a sdest(ab) -> b */
987/* */
988/* sbond: Bond two subsegments together at the respective origins. */
989/* sbond(ab, ac) */
990/* */
991/* */
992/* For interacting tetrahedra and subfacets: */
993/* */
994/* tspivot: Find a subsegment abutting a triangle. */
995/* tspivot(abc) -> ba */
996/* */
997/* stpivot: Find a triangle abutting a subsegment. */
998/* stpivot(ab) -> ba* */
999/* */
1000/* tsbond: Bond a triangle to a subsegment. */
1001/* tsbond(abc, ba) */
1002/* */
1003/*****************************************************************************/
1004
1005/********* Mesh manipulation primitives begin here *********/
1009/* Fast lookup arrays to speed some of the mesh manipulation primitives. */
1010
1011int plus1mod3[3] = {1, 2, 0};
1012int minus1mod3[3] = {2, 0, 1};
1013
1014/********* Primitives for triangles *********/
1015/* */
1016/* */
1017
1018/* decode() converts a pointer to an oriented triangle. The orientation is */
1019/* extracted from the two least significant bits of the pointer. */
1020
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)
1025
1026/* encode() compresses an oriented triangle into a single pointer. It */
1027/* relies on the assumption that all triangles are aligned to four-byte */
1028/* boundaries, so the two least significant bits of (otri).tri are zero. */
1029
1030#define encode(otri) \
1031 (triangle) ((unsigned long) (otri).tri | (unsigned long) (otri).orient)
1032
1033/* The following handle manipulation primitives are all described by Guibas */
1034/* and Stolfi. However, Guibas and Stolfi use an edge-based data */
1035/* structure, whereas I use a triangle-based data structure. */
1036
1037/* sym() finds the abutting triangle, on the same edge. Note that the edge */
1038/* direction is necessarily reversed, because the handle specified by an */
1039/* oriented triangle is directed counterclockwise around the triangle. */
1040
1041#define sym(otri1, otri2) \
1042 ptr = (otri1).tri[(otri1).orient]; \
1043 decode(ptr, otri2);
1044
1045#define symself(otri) \
1046 ptr = (otri).tri[(otri).orient]; \
1047 decode(ptr, otri);
1048
1049/* lnext() finds the next edge (counterclockwise) of a triangle. */
1050
1051#define lnext(otri1, otri2) \
1052 (otri2).tri = (otri1).tri; \
1053 (otri2).orient = plus1mod3[(otri1).orient]
1054
1055#define lnextself(otri) \
1056 (otri).orient = plus1mod3[(otri).orient]
1057
1058/* lprev() finds the previous edge (clockwise) of a triangle. */
1059
1060#define lprev(otri1, otri2) \
1061 (otri2).tri = (otri1).tri; \
1062 (otri2).orient = minus1mod3[(otri1).orient]
1063
1064#define lprevself(otri) \
1065 (otri).orient = minus1mod3[(otri).orient]
1066
1067/* onext() spins counterclockwise around a vertex; that is, it finds the */
1068/* next edge with the same origin in the counterclockwise direction. This */
1069/* edge is part of a different triangle. */
1070
1071#define onext(otri1, otri2) \
1072 lprev(otri1, otri2); \
1073 symself(otri2);
1074
1075#define onextself(otri) \
1076 lprevself(otri); \
1077 symself(otri);
1078
1079/* oprev() spins clockwise around a vertex; that is, it finds the next edge */
1080/* with the same origin in the clockwise direction. This edge is part of */
1081/* a different triangle. */
1082
1083#define oprev(otri1, otri2) \
1084 sym(otri1, otri2); \
1085 lnextself(otri2);
1086
1087#define oprevself(otri) \
1088 symself(otri); \
1089 lnextself(otri);
1090
1091/* dnext() spins counterclockwise around a vertex; that is, it finds the */
1092/* next edge with the same destination in the counterclockwise direction. */
1093/* This edge is part of a different triangle. */
1094
1095#define dnext(otri1, otri2) \
1096 sym(otri1, otri2); \
1097 lprevself(otri2);
1098
1099#define dnextself(otri) \
1100 symself(otri); \
1101 lprevself(otri);
1102
1103/* dprev() spins clockwise around a vertex; that is, it finds the next edge */
1104/* with the same destination in the clockwise direction. This edge is */
1105/* part of a different triangle. */
1106
1107#define dprev(otri1, otri2) \
1108 lnext(otri1, otri2); \
1109 symself(otri2);
1110
1111#define dprevself(otri) \
1112 lnextself(otri); \
1113 symself(otri);
1114
1115/* rnext() moves one edge counterclockwise about the adjacent triangle. */
1116/* (It's best understood by reading Guibas and Stolfi. It involves */
1117/* changing triangles twice.) */
1118
1119#define rnext(otri1, otri2) \
1120 sym(otri1, otri2); \
1121 lnextself(otri2); \
1122 symself(otri2);
1123
1124#define rnextself(otri) \
1125 symself(otri); \
1126 lnextself(otri); \
1127 symself(otri);
1128
1129/* rprev() moves one edge clockwise about the adjacent triangle. */
1130/* (It's best understood by reading Guibas and Stolfi. It involves */
1131/* changing triangles twice.) */
1132
1133#define rprev(otri1, otri2) \
1134 sym(otri1, otri2); \
1135 lprevself(otri2); \
1136 symself(otri2);
1137
1138#define rprevself(otri) \
1139 symself(otri); \
1140 lprevself(otri); \
1141 symself(otri);
1142
1143/* These primitives determine or set the origin, destination, or apex of a */
1144/* triangle. */
1145
1146#define org(otri, vertexptr) \
1147 vertexptr = (vertex) (otri).tri[plus1mod3[(otri).orient] + 3]
1148
1149#define dest(otri, vertexptr) \
1150 vertexptr = (vertex) (otri).tri[minus1mod3[(otri).orient] + 3]
1151
1152#define apex(otri, vertexptr) \
1153 vertexptr = (vertex) (otri).tri[(otri).orient + 3]
1154
1155#define setorg(otri, vertexptr) \
1156 (otri).tri[plus1mod3[(otri).orient] + 3] = (triangle) vertexptr
1157
1158#define setdest(otri, vertexptr) \
1159 (otri).tri[minus1mod3[(otri).orient] + 3] = (triangle) vertexptr
1160
1161#define setapex(otri, vertexptr) \
1162 (otri).tri[(otri).orient + 3] = (triangle) vertexptr
1163
1164/* Bond two triangles together. */
1165
1166#define bond(otri1, otri2) \
1167 (otri1).tri[(otri1).orient] = encode(otri2); \
1168 (otri2).tri[(otri2).orient] = encode(otri1)
1169
1170/* Dissolve a bond (from one side). Note that the other triangle will still */
1171/* think it's connected to this triangle. Usually, however, the other */
1172/* triangle is being deleted entirely, or bonded to another triangle, so */
1173/* it doesn't matter. */
1174
1175#define dissolve(otri) \
1176 (otri).tri[(otri).orient] = (triangle) m->dummytri
1177
1178/* Copy an oriented triangle. */
1179
1180#define otricopy(otri1, otri2) \
1181 (otri2).tri = (otri1).tri; \
1182 (otri2).orient = (otri1).orient
1183
1184/* Test for equality of oriented triangles. */
1185
1186#define otriequal(otri1, otri2) \
1187 (((otri1).tri == (otri2).tri) && \
1188 ((otri1).orient == (otri2).orient))
1189
1190/* Primitives to infect or cure a triangle with the virus. These rely on */
1191/* the assumption that all subsegments are aligned to four-byte boundaries.*/
1192
1193#define infect(otri) \
1194 (otri).tri[6] = (triangle) \
1195 ((unsigned long) (otri).tri[6] | (unsigned long) 2l)
1196
1197#define uninfect(otri) \
1198 (otri).tri[6] = (triangle) \
1199 ((unsigned long) (otri).tri[6] & ~ (unsigned long) 2l)
1200
1201/* Test a triangle for viral infection. */
1202
1203#define infected(otri) \
1204 (((unsigned long) (otri).tri[6] & (unsigned long) 2l) != 0l)
1205
1206/* Check or set a triangle's attributes. */
1207
1208#define elemattribute(otri, attnum) \
1209 ((REAL *) (otri).tri)[m->elemattribindex + (attnum)]
1210
1211#define setelemattribute(otri, attnum, value) \
1212 ((REAL *) (otri).tri)[m->elemattribindex + (attnum)] = value
1213
1214/* Check or set a triangle's maximum area bound. */
1215
1216#define areabound(otri) ((REAL *) (otri).tri)[m->areaboundindex]
1217
1218#define setareabound(otri, value) \
1219 ((REAL *) (otri).tri)[m->areaboundindex] = value
1220
1221/* Check or set a triangle's deallocation. Its second pointer is set to */
1222/* NULL to indicate that it is not allocated. (Its first pointer is used */
1223/* for the stack of dead items.) Its fourth pointer (its first vertex) */
1224/* is set to NULL in case a `badtriang' structure points to it. */
1225
1226#define deadtri(tria) ((tria)[1] == (triangle) NULL)
1227
1228#define killtri(tria) \
1229 (tria)[1] = (triangle) NULL; \
1230 (tria)[3] = (triangle) NULL
1231
1232/********* Primitives for subsegments *********/
1233/* */
1234/* */
1235
1236/* sdecode() converts a pointer to an oriented subsegment. The orientation */
1237/* is extracted from the least significant bit of the pointer. The two */
1238/* least significant bits (one for orientation, one for viral infection) */
1239/* are masked out to produce the real pointer. */
1240
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)
1245
1246/* sencode() compresses an oriented subsegment into a single pointer. It */
1247/* relies on the assumption that all subsegments are aligned to two-byte */
1248/* boundaries, so the least significant bit of (osub).ss is zero. */
1249
1250#define sencode(osub) \
1251 (subseg) ((unsigned long) (osub).ss | (unsigned long) (osub).ssorient)
1252
1253/* ssym() toggles the orientation of a subsegment. */
1254
1255#define ssym(osub1, osub2) \
1256 (osub2).ss = (osub1).ss; \
1257 (osub2).ssorient = 1 - (osub1).ssorient
1258
1259#define ssymself(osub) \
1260 (osub).ssorient = 1 - (osub).ssorient
1261
1262/* spivot() finds the other subsegment (from the same segment) that shares */
1263/* the same origin. */
1264
1265#define spivot(osub1, osub2) \
1266 sptr = (osub1).ss[(osub1).ssorient]; \
1267 sdecode(sptr, osub2)
1268
1269#define spivotself(osub) \
1270 sptr = (osub).ss[(osub).ssorient]; \
1271 sdecode(sptr, osub)
1272
1273/* snext() finds the next subsegment (from the same segment) in sequence; */
1274/* one whose origin is the input subsegment's destination. */
1275
1276#define snext(osub1, osub2) \
1277 sptr = (osub1).ss[1 - (osub1).ssorient]; \
1278 sdecode(sptr, osub2)
1279
1280#define snextself(osub) \
1281 sptr = (osub).ss[1 - (osub).ssorient]; \
1282 sdecode(sptr, osub)
1283
1284/* These primitives determine or set the origin or destination of a */
1285/* subsegment or the segment that includes it. */
1286
1287#define sorg(osub, vertexptr) \
1288 vertexptr = (vertex) (osub).ss[2 + (osub).ssorient]
1289
1290#define sdest(osub, vertexptr) \
1291 vertexptr = (vertex) (osub).ss[3 - (osub).ssorient]
1292
1293#define setsorg(osub, vertexptr) \
1294 (osub).ss[2 + (osub).ssorient] = (subseg) vertexptr
1295
1296#define setsdest(osub, vertexptr) \
1297 (osub).ss[3 - (osub).ssorient] = (subseg) vertexptr
1298
1299#define segorg(osub, vertexptr) \
1300 vertexptr = (vertex) (osub).ss[4 + (osub).ssorient]
1301
1302#define segdest(osub, vertexptr) \
1303 vertexptr = (vertex) (osub).ss[5 - (osub).ssorient]
1304
1305#define setsegorg(osub, vertexptr) \
1306 (osub).ss[4 + (osub).ssorient] = (subseg) vertexptr
1307
1308#define setsegdest(osub, vertexptr) \
1309 (osub).ss[5 - (osub).ssorient] = (subseg) vertexptr
1310
1311/* These primitives read or set a boundary marker. Boundary markers are */
1312/* used to hold user-defined tags for setting boundary conditions in */
1313/* finite element solvers. */
1314
1315#define mark(osub) (* (int *) ((osub).ss + 8))
1316
1317#define setmark(osub, value) \
1318 * (int *) ((osub).ss + 8) = value
1319
1320/* Bond two subsegments together. */
1321
1322#define sbond(osub1, osub2) \
1323 (osub1).ss[(osub1).ssorient] = sencode(osub2); \
1324 (osub2).ss[(osub2).ssorient] = sencode(osub1)
1325
1326/* Dissolve a subsegment bond (from one side). Note that the other */
1327/* subsegment will still think it's connected to this subsegment. */
1328
1329#define sdissolve(osub) \
1330 (osub).ss[(osub).ssorient] = (subseg) m->dummysub
1331
1332/* Copy a subsegment. */
1333
1334#define subsegcopy(osub1, osub2) \
1335 (osub2).ss = (osub1).ss; \
1336 (osub2).ssorient = (osub1).ssorient
1337
1338/* Test for equality of subsegments. */
1339
1340#define subsegequal(osub1, osub2) \
1341 (((osub1).ss == (osub2).ss) && \
1342 ((osub1).ssorient == (osub2).ssorient))
1343
1344/* Check or set a subsegment's deallocation. Its second pointer is set to */
1345/* NULL to indicate that it is not allocated. (Its first pointer is used */
1346/* for the stack of dead items.) Its third pointer (its first vertex) */
1347/* is set to NULL in case a `badsubseg' structure points to it. */
1348
1349#define deadsubseg(sub) ((sub)[1] == (subseg) NULL)
1350
1351#define killsubseg(sub) \
1352 (sub)[1] = (subseg) NULL; \
1353 (sub)[2] = (subseg) NULL
1354
1355/********* Primitives for interacting triangles and subsegments *********/
1356/* */
1357/* */
1358
1359/* tspivot() finds a subsegment abutting a triangle. */
1360
1361#define tspivot(otri, osub) \
1362 sptr = (subseg) (otri).tri[6 + (otri).orient]; \
1363 sdecode(sptr, osub)
1364
1365/* stpivot() finds a triangle abutting a subsegment. It requires that the */
1366/* variable `ptr' of type `triangle' be defined. */
1367
1368#define stpivot(osub, otri) \
1369 ptr = (triangle) (osub).ss[6 + (osub).ssorient]; \
1370 decode(ptr, otri)
1371
1372/* Bond a triangle to a subsegment. */
1373
1374#define tsbond(otri, osub) \
1375 (otri).tri[6 + (otri).orient] = (triangle) sencode(osub); \
1376 (osub).ss[6 + (osub).ssorient] = (subseg) encode(otri)
1377
1378/* Dissolve a bond (from the triangle side). */
1379
1380#define tsdissolve(otri) \
1381 (otri).tri[6 + (otri).orient] = (triangle) m->dummysub
1382
1383/* Dissolve a bond (from the subsegment side). */
1384
1385#define stdissolve(osub) \
1386 (osub).ss[6 + (osub).ssorient] = (subseg) m->dummytri
1387
1388/********* Primitives for vertices *********/
1389/* */
1390/* */
1391
1392#define vertexmark(vx) ((int *) (vx))[m->vertexmarkindex]
1393
1394#define setvertexmark(vx, value) \
1395 ((int *) (vx))[m->vertexmarkindex] = value
1396
1397#define vertextype(vx) ((int *) (vx))[m->vertexmarkindex + 1]
1398
1399#define setvertextype(vx, value) \
1400 ((int *) (vx))[m->vertexmarkindex + 1] = value
1401
1402#define vertex2tri(vx) ((triangle *) (vx))[m->vertex2triindex]
1403
1404#define setvertex2tri(vx, value) \
1405 ((triangle *) (vx))[m->vertex2triindex] = value
1406
1409/********* Mesh manipulation primitives end here *********/
1410
1411/********* User-defined triangle evaluation routine begins here *********/
1415/*****************************************************************************/
1416/* */
1417/* triunsuitable() Determine if a triangle is unsuitable, and thus must */
1418/* be further refined. */
1419/* */
1420/* You may write your own procedure that decides whether or not a selected */
1421/* triangle is too big (and needs to be refined). There are two ways to do */
1422/* this. */
1423/* */
1424/* (1) Modify the procedure `triunsuitable' below, then recompile */
1425/* Triangle. */
1426/* */
1427/* (2) Define the symbol EXTERNAL_TEST (either by adding the definition */
1428/* to this file, or by using the appropriate compiler switch). This way, */
1429/* you can compile triangle.c separately from your test. Write your own */
1430/* `triunsuitable' procedure in a separate C file (using the same prototype */
1431/* as below). Compile it and link the object code with triangle.o. */
1432/* */
1433/* This procedure returns 1 if the triangle is too large and should be */
1434/* refined; 0 otherwise. */
1435/* */
1436/*****************************************************************************/
1437
1438#ifdef EXTERNAL_TEST
1439
1440int triunsuitable();
1441
1442#else /* not EXTERNAL_TEST */
1443
1444#ifdef ANSI_DECLARATORS
1445int triunsuitable(vertex triorg, vertex tridest, vertex triapex, REAL area)
1446#else /* not ANSI_DECLARATORS */
1447int triunsuitable(triorg, tridest, triapex, area)
1448vertex triorg; /* The triangle's origin vertex. */
1449vertex tridest; /* The triangle's destination vertex. */
1450vertex triapex; /* The triangle's apex vertex. */
1451REAL area; /* The area of the triangle. */
1452#endif /* not ANSI_DECLARATORS */
1453
1454{
1455 REAL dxoa, dxda, dxod;
1456 REAL dyoa, dyda, dyod;
1457 REAL oalen, dalen, odlen;
1458 REAL maxlen;
1459
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];
1466 /* Find the squares of the lengths of the triangle's three edges. */
1467 oalen = dxoa * dxoa + dyoa * dyoa;
1468 dalen = dxda * dxda + dyda * dyda;
1469 odlen = dxod * dxod + dyod * dyod;
1470 /* Find the square of the length of the longest edge. */
1471 maxlen = (dalen > oalen) ? dalen : oalen;
1472 maxlen = (odlen > maxlen) ? odlen : maxlen;
1473
1474 if (maxlen > 0.05 * (triorg[0] * triorg[0] + triorg[1] * triorg[1]) + 0.02) {
1475 return 1;
1476 } else {
1477 return 0;
1478 }
1479}
1480
1481#endif /* not EXTERNAL_TEST */
1482
1485/********* User-defined triangle evaluation routine ends here *********/
1486
1487/********* Memory allocation and program exit wrappers begin here *********/
1491#ifdef ANSI_DECLARATORS
1492void triexit(int status)
1493#else /* not ANSI_DECLARATORS */
1494void triexit(status)
1495int status;
1496#endif /* not ANSI_DECLARATORS */
1497
1498{
1499 exit(status);
1500}
1501
1502#ifdef ANSI_DECLARATORS
1503VOID *trimalloc(int size)
1504#else /* not ANSI_DECLARATORS */
1505VOID *trimalloc(size)
1506int size;
1507#endif /* not ANSI_DECLARATORS */
1508
1509{
1510 VOID *memptr;
1511
1512 memptr = (VOID *) malloc((unsigned int) size);
1513 if (memptr == (VOID *) NULL) {
1514 printf("Error: Out of memory.\n");
1515 triexit(1);
1516 }
1517 return(memptr);
1518}
1519
1520#ifdef ANSI_DECLARATORS
1521void trifree(VOID *memptr)
1522#else /* not ANSI_DECLARATORS */
1523void trifree(memptr)
1524VOID *memptr;
1525#endif /* not ANSI_DECLARATORS */
1526
1527{
1528 free(memptr);
1529}
1530
1533/********* Memory allocation and program exit wrappers end here *********/
1534
1535/********* User interaction routines begin here *********/
1539/*****************************************************************************/
1540/* */
1541/* syntax() Print list of command line switches. */
1542/* */
1543/*****************************************************************************/
1544
1545#ifndef TRILIBRARY
1546
1547void syntax()
1548{
1549#ifdef CDT_ONLY
1550#ifdef REDUCED
1551 printf("triangle [-pAcjevngBPNEIOXzo_lQVh] input_file\n");
1552#else /* not REDUCED */
1553 printf("triangle [-pAcjevngBPNEIOXzo_iFlCQVh] input_file\n");
1554#endif /* not REDUCED */
1555#else /* not CDT_ONLY */
1556#ifdef REDUCED
1557 printf("triangle [-prq__a__uAcDjevngBPNEIOXzo_YS__lQVh] input_file\n");
1558#else /* not REDUCED */
1559 printf("triangle [-prq__a__uAcDjevngBPNEIOXzo_YS__iFlsCQVh] input_file\n");
1560#endif /* not REDUCED */
1561#endif /* not CDT_ONLY */
1562
1563 printf(" -p Triangulates a Planar Straight Line Graph (.poly file).\n");
1564#ifndef CDT_ONLY
1565 printf(" -r Refines a previously generated mesh.\n");
1566 printf(
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");
1570#endif /* not CDT_ONLY */
1571 printf(
1572 " -A Applies attributes to identify triangles in certain regions.\n");
1573 printf(" -c Encloses the convex hull with segments.\n");
1574#ifndef CDT_ONLY
1575 printf(" -D Conforming Delaunay: all triangles are truly Delaunay.\n");
1576#endif /* not CDT_ONLY */
1577/*
1578 printf(" -w Weighted Delaunay triangulation.\n");
1579 printf(" -W Regular triangulation (lower hull of a height field).\n");
1580*/
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");
1595#ifndef CDT_ONLY
1596 printf(" -Y Suppresses boundary segment splitting.\n");
1597 printf(" -S Specifies maximum number of added Steiner points.\n");
1598#endif /* not CDT_ONLY */
1599#ifndef REDUCED
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");
1602#endif /* not REDUCED */
1603 printf(" -l Uses vertical cuts only, rather than alternating cuts.\n");
1604#ifndef REDUCED
1605#ifndef CDT_ONLY
1606 printf(
1607 " -s Force segments into mesh by splitting (instead of using CDT).\n");
1608#endif /* not CDT_ONLY */
1609 printf(" -C Check consistency of final mesh.\n");
1610#endif /* not REDUCED */
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");
1614 triexit(0);
1615}
1616
1617#endif /* not TRILIBRARY */
1618
1619/*****************************************************************************/
1620/* */
1621/* info() Print out complete instructions. */
1622/* */
1623/*****************************************************************************/
1624
1625#ifndef TRILIBRARY
1626
1627void info()
1628{
1629 printf("Triangle\n");
1630 printf(
1631"A Two-Dimensional Quality Mesh Generator and Delaunay Triangulator.\n");
1632 printf("Version 1.6\n\n");
1633 printf(
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");
1637 printf(
1638"Created as part of the Quake project (tools for earthquake simulation).\n");
1639 printf(
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");
1642#ifdef SINGLE
1643 printf("This executable is compiled for single precision arithmetic.\n\n\n");
1644#else /* not SINGLE */
1645 printf("This executable is compiled for double precision arithmetic.\n\n\n");
1646#endif /* not SINGLE */
1647 printf(
1648"Triangle generates exact Delaunay triangulations, constrained Delaunay\n");
1649 printf(
1650"triangulations, conforming Delaunay triangulations, Voronoi diagrams, and\n");
1651 printf(
1652"high-quality triangular meshes. The latter can be generated with no small\n"
1653);
1654 printf(
1655"or large angles, and are thus suitable for finite element analysis. If no\n"
1656);
1657 printf(
1658"command line switch is specified, your .node input file is read, and the\n");
1659 printf(
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");
1663 printf(
1664"Underscores indicate that numbers may optionally follow certain switches.\n");
1665 printf(
1666"Do not leave any space between a switch and its numeric parameter.\n");
1667 printf(
1668"input_file must be a file with extension .node, or extension .poly if the\n");
1669 printf(
1670"-p switch is used. If -r is used, you must supply .node and .ele files,\n");
1671 printf(
1672"and possibly a .poly file and an .area file as well. The formats of these\n"
1673);
1674 printf("files are described below.\n\n");
1675 printf("Command Line Switches:\n\n");
1676 printf(
1677" -p Reads a Planar Straight Line Graph (.poly file), which can specify\n"
1678);
1679 printf(
1680" vertices, segments, holes, regional attributes, and regional area\n");
1681 printf(
1682" constraints. Generates a constrained Delaunay triangulation (CDT)\n"
1683);
1684 printf(
1685" fitting the input; or, if -s, -q, -a, or -u is used, a conforming\n");
1686 printf(
1687" constrained Delaunay triangulation (CCDT). If you want a truly\n");
1688 printf(
1689" Delaunay (not just constrained Delaunay) triangulation, use -D as\n");
1690 printf(
1691" well. When -p is not used, Triangle reads a .node file by default.\n"
1692);
1693 printf(
1694" -r Refines a previously generated mesh. The mesh is read from a .node\n"
1695);
1696 printf(
1697" file and an .ele file. If -p is also used, a .poly file is read\n");
1698 printf(
1699" and used to constrain segments in the mesh. If -a is also used\n");
1700 printf(
1701" (with no number following), an .area file is read and used to\n");
1702 printf(
1703" impose area constraints on the mesh. Further details on refinement\n"
1704);
1705 printf(" appear below.\n");
1706 printf(
1707" -q Quality mesh generation by Delaunay refinement (a hybrid of Paul\n");
1708 printf(
1709" Chew's and Jim Ruppert's algorithms). Adds vertices to the mesh to\n"
1710);
1711 printf(
1712" ensure that all angles are between 20 and 140 degrees. An\n");
1713 printf(
1714" alternative bound on the minimum angle, replacing 20 degrees, may\n");
1715 printf(
1716" be specified after the `q'. The specified angle may include a\n");
1717 printf(
1718" decimal point, but not exponential notation. Note that a bound of\n"
1719);
1720 printf(
1721" theta degrees on the smallest angle also implies a bound of\n");
1722 printf(
1723" (180 - 2 theta) on the largest angle. If the minimum angle is 28.6\n"
1724);
1725 printf(
1726" degrees or smaller, Triangle is mathematically guaranteed to\n");
1727 printf(
1728" terminate (assuming infinite precision arithmetic--Triangle may\n");
1729 printf(
1730" fail to terminate if you run out of precision). In practice,\n");
1731 printf(
1732" Triangle often succeeds for minimum angles up to 34 degrees. For\n");
1733 printf(
1734" some meshes, however, you might need to reduce the minimum angle to\n"
1735);
1736 printf(
1737" avoid problems associated with insufficient floating-point\n");
1738 printf(" precision.\n");
1739 printf(
1740" -a Imposes a maximum triangle area. If a number follows the `a', no\n");
1741 printf(
1742" triangle is generated whose area is larger than that number. If no\n"
1743);
1744 printf(
1745" number is specified, an .area file (if -r is used) or .poly file\n");
1746 printf(
1747" (if -r is not used) specifies a set of maximum area constraints.\n");
1748 printf(
1749" An .area file contains a separate area constraint for each\n");
1750 printf(
1751" triangle, and is useful for refining a finite element mesh based on\n"
1752);
1753 printf(
1754" a posteriori error estimates. A .poly file can optionally contain\n"
1755);
1756 printf(
1757" an area constraint for each segment-bounded region, thereby\n");
1758 printf(
1759" controlling triangle densities in a first triangulation of a PSLG.\n"
1760);
1761 printf(
1762" You can impose both a fixed area constraint and a varying area\n");
1763 printf(
1764" constraint by invoking the -a switch twice, once with and once\n");
1765 printf(
1766" without a number following. Each area specified may include a\n");
1767 printf(" decimal point.\n");
1768 printf(
1769" -u Imposes a user-defined constraint on triangle size. There are two\n"
1770);
1771 printf(
1772" ways to use this feature. One is to edit the triunsuitable()\n");
1773 printf(
1774" procedure in triangle.c to encode any constraint you like, then\n");
1775 printf(
1776" recompile Triangle. The other is to compile triangle.c with the\n");
1777 printf(
1778" EXTERNAL_TEST symbol set (compiler switch -DEXTERNAL_TEST), then\n");
1779 printf(
1780" link Triangle with a separate object file that implements\n");
1781 printf(
1782" triunsuitable(). In either case, the -u switch causes the user-\n");
1783 printf(" defined test to be applied to every triangle.\n");
1784 printf(
1785" -A Assigns an additional floating-point attribute to each triangle\n");
1786 printf(
1787" that identifies what segment-bounded region each triangle belongs\n");
1788 printf(
1789" to. Attributes are assigned to regions by the .poly file. If a\n");
1790 printf(
1791" region is not explicitly marked by the .poly file, triangles in\n");
1792 printf(
1793" that region are assigned an attribute of zero. The -A switch has\n");
1794 printf(
1795" an effect only when the -p switch is used and the -r switch is not.\n"
1796);
1797 printf(
1798" -c Creates segments on the convex hull of the triangulation. If you\n");
1799 printf(
1800" are triangulating a vertex set, this switch causes a .poly file to\n"
1801);
1802 printf(
1803" be written, containing all edges of the convex hull. If you are\n");
1804 printf(
1805" triangulating a PSLG, this switch specifies that the whole convex\n");
1806 printf(
1807" hull of the PSLG should be triangulated, regardless of what\n");
1808 printf(
1809" segments the PSLG has. If you do not use this switch when\n");
1810 printf(
1811" triangulating a PSLG, Triangle assumes that you have identified the\n"
1812);
1813 printf(
1814" region to be triangulated by surrounding it with segments of the\n");
1815 printf(
1816" input PSLG. Beware: if you are not careful, this switch can cause\n"
1817);
1818 printf(
1819" the introduction of an extremely thin angle between a PSLG segment\n"
1820);
1821 printf(
1822" and a convex hull segment, which can cause overrefinement (and\n");
1823 printf(
1824" possibly failure if Triangle runs out of precision). If you are\n");
1825 printf(
1826" refining a mesh, the -c switch works differently: it causes a\n");
1827 printf(
1828" .poly file to be written containing the boundary edges of the mesh\n"
1829);
1830 printf(" (useful if no .poly file was read).\n");
1831 printf(
1832" -D Conforming Delaunay triangulation: use this switch if you want to\n"
1833);
1834 printf(
1835" ensure that all the triangles in the mesh are Delaunay, and not\n");
1836 printf(
1837" merely constrained Delaunay; or if you want to ensure that all the\n"
1838);
1839 printf(
1840" Voronoi vertices lie within the triangulation. (Some finite volume\n"
1841);
1842 printf(
1843" methods have this requirement.) This switch invokes Ruppert's\n");
1844 printf(
1845" original algorithm, which splits every subsegment whose diametral\n");
1846 printf(
1847" circle is encroached. It usually increases the number of vertices\n"
1848);
1849 printf(" and triangles.\n");
1850 printf(
1851" -j Jettisons vertices that are not part of the final triangulation\n");
1852 printf(
1853" from the output .node file. By default, Triangle copies all\n");
1854 printf(
1855" vertices in the input .node file to the output .node file, in the\n");
1856 printf(
1857" same order, so their indices do not change. The -j switch prevents\n"
1858);
1859 printf(
1860" duplicated input vertices, or vertices `eaten' by holes, from\n");
1861 printf(
1862" appearing in the output .node file. Thus, if two input vertices\n");
1863 printf(
1864" have exactly the same coordinates, only the first appears in the\n");
1865 printf(
1866" output. If any vertices are jettisoned, the vertex numbering in\n");
1867 printf(
1868" the output .node file differs from that of the input .node file.\n");
1869 printf(
1870" -e Outputs (to an .edge file) a list of edges of the triangulation.\n");
1871 printf(
1872" -v Outputs the Voronoi diagram associated with the triangulation.\n");
1873 printf(
1874" Does not attempt to detect degeneracies, so some Voronoi vertices\n");
1875 printf(
1876" may be duplicated. See the discussion of Voronoi diagrams below.\n");
1877 printf(
1878" -n Outputs (to a .neigh file) a list of triangles neighboring each\n");
1879 printf(" triangle.\n");
1880 printf(
1881" -g Outputs the mesh to an Object File Format (.off) file, suitable for\n"
1882);
1883 printf(" viewing with the Geometry Center's Geomview package.\n");
1884 printf(
1885" -B No boundary markers in the output .node, .poly, and .edge output\n");
1886 printf(
1887" files. See the detailed discussion of boundary markers below.\n");
1888 printf(
1889" -P No output .poly file. Saves disk space, but you lose the ability\n");
1890 printf(
1891" to maintain constraining segments on later refinements of the mesh.\n"
1892);
1893 printf(" -N No output .node file.\n");
1894 printf(" -E No output .ele file.\n");
1895 printf(
1896" -I No iteration numbers. Suppresses the output of .node and .poly\n");
1897 printf(
1898" files, so your input files won't be overwritten. (If your input is\n"
1899);
1900 printf(
1901" a .poly file only, a .node file is written.) Cannot be used with\n");
1902 printf(
1903" the -r switch, because that would overwrite your input .ele file.\n");
1904 printf(
1905" Shouldn't be used with the -q, -a, -u, or -s switch if you are\n");
1906 printf(
1907" using a .node file for input, because no .node file is written, so\n"
1908);
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");
1911 printf(
1912" -X No exact arithmetic. Normally, Triangle uses exact floating-point\n"
1913);
1914 printf(
1915" arithmetic for certain tests if it thinks the inexact tests are not\n"
1916);
1917 printf(
1918" accurate enough. Exact arithmetic ensures the robustness of the\n");
1919 printf(
1920" triangulation algorithms, despite floating-point roundoff error.\n");
1921 printf(
1922" Disabling exact arithmetic with the -X switch causes a small\n");
1923 printf(
1924" improvement in speed and creates the possibility that Triangle will\n"
1925);
1926 printf(" fail to produce a valid mesh. Not recommended.\n");
1927 printf(
1928" -z Numbers all items starting from zero (rather than one). Note that\n"
1929);
1930 printf(
1931" this switch is normally overridden by the value used to number the\n"
1932);
1933 printf(
1934" first vertex of the input .node or .poly file. However, this\n");
1935 printf(
1936" switch is useful when calling Triangle from another program.\n");
1937 printf(
1938" -o2 Generates second-order subparametric elements with six nodes each.\n"
1939);
1940 printf(
1941" -Y No new vertices on the boundary. This switch is useful when the\n");
1942 printf(
1943" mesh boundary must be preserved so that it conforms to some\n");
1944 printf(
1945" adjacent mesh. Be forewarned that you will probably sacrifice much\n"
1946);
1947 printf(
1948" of the quality of the mesh; Triangle will try, but the resulting\n");
1949 printf(
1950" mesh may contain poorly shaped triangles. Works well if all the\n");
1951 printf(
1952" boundary vertices are closely spaced. Specify this switch twice\n");
1953 printf(
1954" (`-YY') to prevent all segment splitting, including internal\n");
1955 printf(" boundaries.\n");
1956 printf(
1957" -S Specifies the maximum number of Steiner points (vertices that are\n");
1958 printf(
1959" not in the input, but are added to meet the constraints on minimum\n"
1960);
1961 printf(
1962" angle and maximum area). The default is to allow an unlimited\n");
1963 printf(
1964" number. If you specify this switch with no number after it,\n");
1965 printf(
1966" the limit is set to zero. Triangle always adds vertices at segment\n"
1967);
1968 printf(
1969" intersections, even if it needs to use more vertices than the limit\n"
1970);
1971 printf(
1972" you set. When Triangle inserts segments by splitting (-s), it\n");
1973 printf(
1974" always adds enough vertices to ensure that all the segments of the\n"
1975);
1976 printf(" PLSG are recovered, ignoring the limit if necessary.\n");
1977 printf(
1978" -i Uses an incremental rather than a divide-and-conquer algorithm to\n");
1979 printf(
1980" construct a Delaunay triangulation. Try it if the divide-and-\n");
1981 printf(" conquer algorithm fails.\n");
1982 printf(
1983" -F Uses Steven Fortune's sweepline algorithm to construct a Delaunay\n");
1984 printf(
1985" triangulation. Warning: does not use exact arithmetic for all\n");
1986 printf(" calculations. An exact result is not guaranteed.\n");
1987 printf(
1988" -l Uses only vertical cuts in the divide-and-conquer algorithm. By\n");
1989 printf(
1990" default, Triangle alternates between vertical and horizontal cuts,\n"
1991);
1992 printf(
1993" which usually improve the speed except with vertex sets that are\n");
1994 printf(
1995" small or short and wide. This switch is primarily of theoretical\n");
1996 printf(" interest.\n");
1997 printf(
1998" -s Specifies that segments should be forced into the triangulation by\n"
1999);
2000 printf(
2001" recursively splitting them at their midpoints, rather than by\n");
2002 printf(
2003" generating a constrained Delaunay triangulation. Segment splitting\n"
2004);
2005 printf(
2006" is true to Ruppert's original algorithm, but can create needlessly\n"
2007);
2008 printf(
2009" small triangles. This switch is primarily of theoretical interest.\n"
2010);
2011 printf(
2012" -C Check the consistency of the final mesh. Uses exact arithmetic for\n"
2013);
2014 printf(
2015" checking, even if the -X switch is used. Useful if you suspect\n");
2016 printf(" Triangle is buggy.\n");
2017 printf(
2018" -Q Quiet: Suppresses all explanation of what Triangle is doing,\n");
2019 printf(" unless an error occurs.\n");
2020 printf(
2021" -V Verbose: Gives detailed information about what Triangle is doing.\n"
2022);
2023 printf(
2024" Add more `V's for increasing amount of detail. `-V' is most\n");
2025 printf(
2026" useful; itgives information on algorithmic progress and much more\n");
2027 printf(
2028" detailed statistics. `-VV' gives vertex-by-vertex details, and\n");
2029 printf(
2030" prints so much that Triangle runs much more slowly. `-VVVV' gives\n"
2031);
2032 printf(" information only a debugger could love.\n");
2033 printf(" -h Help: Displays these instructions.\n");
2034 printf("\n");
2035 printf("Definitions:\n");
2036 printf("\n");
2037 printf(
2038" A Delaunay triangulation of a vertex set is a triangulation whose\n");
2039 printf(
2040" vertices are the vertex set, that covers the convex hull of the vertex\n");
2041 printf(
2042" set. A Delaunay triangulation has the property that no vertex lies\n");
2043 printf(
2044" inside the circumscribing circle (circle that passes through all three\n");
2045 printf(" vertices) of any triangle in the triangulation.\n\n");
2046 printf(
2047" A Voronoi diagram of a vertex set is a subdivision of the plane into\n");
2048 printf(
2049" polygonal cells (some of which may be unbounded, meaning infinitely\n");
2050 printf(
2051" large), where each cell is the set of points in the plane that are closer\n"
2052);
2053 printf(
2054" to some input vertex than to any other input vertex. The Voronoi diagram\n"
2055);
2056 printf(" is a geometric dual of the Delaunay triangulation.\n\n");
2057 printf(
2058" A Planar Straight Line Graph (PSLG) is a set of vertices and segments.\n");
2059 printf(
2060" Segments are simply edges, whose endpoints are all vertices in the PSLG.\n"
2061);
2062 printf(
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");
2065 printf(
2066" A constrained Delaunay triangulation (CDT) of a PSLG is similar to a\n");
2067 printf(
2068" Delaunay triangulation, but each PSLG segment is present as a single edge\n"
2069);
2070 printf(
2071" of the CDT. (A constrained Delaunay triangulation is not truly a\n");
2072 printf(
2073" Delaunay triangulation, because some of its triangles might not be\n");
2074 printf(
2075" Delaunay.) By definition, a CDT does not have any vertices other than\n");
2076 printf(
2077" those specified in the input PSLG. Depending on context, a CDT might\n");
2078 printf(
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");
2081 printf(
2082" A conforming Delaunay triangulation of a PSLG is a triangulation in which\n"
2083);
2084 printf(
2085" each triangle is truly Delaunay, and each PSLG segment is represented by\n"
2086);
2087 printf(
2088" a linear contiguous sequence of edges of the triangulation. New vertices\n"
2089);
2090 printf(
2091" (not part of the PSLG) may appear, and each input segment may have been\n");
2092 printf(
2093" subdivided into shorter edges (subsegments) by these additional vertices.\n"
2094);
2095 printf(
2096" The new vertices are frequently necessary to maintain the Delaunay\n");
2097 printf(" property while ensuring that every segment is represented.\n\n");
2098 printf(
2099" A conforming constrained Delaunay triangulation (CCDT) of a PSLG is a\n");
2100 printf(
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");
2103 printf(
2104" subsegments, but not to guarantee that segments are respected; rather, to\n"
2105);
2106 printf(
2107" improve the quality of the triangles. The high-quality meshes produced\n");
2108 printf(
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");
2112 printf(
2113" All files may contain comments prefixed by the character '#'. Vertices,\n"
2114);
2115 printf(
2116" triangles, edges, holes, and maximum area constraints must be numbered\n");
2117 printf(
2118" consecutively, starting from either 1 or 0. Whichever you choose, all\n");
2119 printf(
2120" input files must be consistent; if the vertices are numbered from 1, so\n");
2121 printf(
2122" must be all other objects. Triangle automatically detects your choice\n");
2123 printf(
2124" while reading the .node (or .poly) file. (When calling Triangle from\n");
2125 printf(
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");
2129 printf(
2130" First line: <# of vertices> <dimension (must be 2)> <# of attributes>\n"
2131);
2132 printf(
2133" <# of boundary markers (0 or 1)>\n"
2134);
2135 printf(
2136" Remaining lines: <vertex #> <x> <y> [attributes] [boundary marker]\n");
2137 printf("\n");
2138 printf(
2139" The attributes, which are typically floating-point values of physical\n");
2140 printf(
2141" quantities (such as mass or conductivity) associated with the nodes of\n"
2142);
2143 printf(
2144" a finite element mesh, are copied unchanged to the output mesh. If -q,\n"
2145);
2146 printf(
2147" -a, -u, -D, or -s is selected, each new Steiner point added to the mesh\n"
2148);
2149 printf(" has attributes assigned to it by linear interpolation.\n\n");
2150 printf(
2151" If the fourth entry of the first line is `1', the last column of the\n");
2152 printf(
2153" remainder of the file is assumed to contain boundary markers. Boundary\n"
2154);
2155 printf(
2156" markers are used to identify boundary vertices and vertices resting on\n"
2157);
2158 printf(
2159" PSLG segments; a complete description appears in a section below. The\n"
2160);
2161 printf(
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");
2165 printf(
2166" First line: <# of triangles> <nodes per triangle> <# of attributes>\n");
2167 printf(
2168" Remaining lines: <triangle #> <node> <node> <node> ... [attributes]\n");
2169 printf("\n");
2170 printf(
2171" Nodes are indices into the corresponding .node file. The first three\n");
2172 printf(
2173" nodes are the corner vertices, and are listed in counterclockwise order\n"
2174);
2175 printf(
2176" around each triangle. (The remaining nodes, if any, depend on the type\n"
2177);
2178 printf(" of finite element used.)\n\n");
2179 printf(
2180" The attributes are just like those of .node files. Because there is no\n"
2181);
2182 printf(
2183" simple mapping from input to output triangles, Triangle attempts to\n");
2184 printf(
2185" interpolate attributes, and may cause a lot of diffusion of attributes\n"
2186);
2187 printf(
2188" among nearby triangles as the triangulation is refined. Attributes do\n"
2189);
2190 printf(" not diffuse across segments, so attributes used to identify\n");
2191 printf(" segment-bounded regions remain intact.\n\n");
2192 printf(
2193" In .ele files produced by Triangle, each triangular element has three\n");
2194 printf(
2195" nodes (vertices) unless the -o2 switch is used, in which case\n");
2196 printf(
2197" subparametric quadratic elements with six nodes each are generated.\n");
2198 printf(
2199" The first three nodes are the corners in counterclockwise order, and\n");
2200 printf(
2201" the fourth, fifth, and sixth nodes lie on the midpoints of the edges\n");
2202 printf(
2203" opposite the first, second, and third vertices, respectively.\n");
2204 printf("\n");
2205 printf(" .poly files:\n");
2206 printf(
2207" First line: <# of vertices> <dimension (must be 2)> <# of attributes>\n"
2208);
2209 printf(
2210" <# of boundary markers (0 or 1)>\n"
2211);
2212 printf(
2213" Following lines: <vertex #> <x> <y> [attributes] [boundary marker]\n");
2214 printf(" One line: <# of segments> <# of boundary markers (0 or 1)>\n");
2215 printf(
2216" Following lines: <segment #> <endpoint> <endpoint> [boundary marker]\n");
2217 printf(" One line: <# of holes>\n");
2218 printf(" Following lines: <hole #> <x> <y>\n");
2219 printf(
2220" Optional line: <# of regional attributes and/or area constraints>\n");
2221 printf(
2222" Optional following lines: <region #> <x> <y> <attribute> <max area>\n");
2223 printf("\n");
2224 printf(
2225" A .poly file represents a PSLG, as well as some additional information.\n"
2226);
2227 printf(
2228" The first section lists all the vertices, and is identical to the\n");
2229 printf(
2230" format of .node files. <# of vertices> may be set to zero to indicate\n"
2231);
2232 printf(
2233" that the vertices are listed in a separate .node file; .poly files\n");
2234 printf(
2235" produced by Triangle always have this format. A vertex set represented\n"
2236);
2237 printf(
2238" this way has the advantage that it may easily be triangulated with or\n");
2239 printf(
2240" without segments (depending on whether the -p switch is invoked).\n");
2241 printf("\n");
2242 printf(
2243" The second section lists the segments. Segments are edges whose\n");
2244 printf(
2245" presence in the triangulation is enforced. (Depending on the choice of\n"
2246);
2247 printf(
2248" switches, segment might be subdivided into smaller edges). Each\n");
2249 printf(
2250" segment is specified by listing the indices of its two endpoints. This\n"
2251);
2252 printf(
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");
2255 printf(
2256" If -q, -a, -u, and -s are not selected, Triangle produces a constrained\n"
2257);
2258 printf(
2259" Delaunay triangulation (CDT), in which each segment appears as a single\n"
2260);
2261 printf(
2262" edge in the triangulation. If -q, -a, -u, or -s is selected, Triangle\n"
2263);
2264 printf(
2265" produces a conforming constrained Delaunay triangulation (CCDT), in\n");
2266 printf(
2267" which segments may be subdivided into smaller edges. If -D is\n");
2268 printf(
2269" selected, Triangle produces a conforming Delaunay triangulation, so\n");
2270 printf(
2271" that every triangle is Delaunay, and not just constrained Delaunay.\n");
2272 printf("\n");
2273 printf(
2274" The third section lists holes (and concavities, if -c is selected) in\n");
2275 printf(
2276" the triangulation. Holes are specified by identifying a point inside\n");
2277 printf(
2278" each hole. After the triangulation is formed, Triangle creates holes\n");
2279 printf(
2280" by eating triangles, spreading out from each hole point until its\n");
2281 printf(
2282" progress is blocked by segments in the PSLG. You must be careful to\n");
2283 printf(
2284" enclose each hole in segments, or your whole triangulation might be\n");
2285 printf(
2286" eaten away. If the two triangles abutting a segment are eaten, the\n");
2287 printf(
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");
2291 printf(
2292" The optional fourth section lists regional attributes (to be assigned\n");
2293 printf(
2294" to all triangles in a region) and regional constraints on the maximum\n");
2295 printf(
2296" triangle area. Triangle reads this section only if the -A switch is\n");
2297 printf(
2298" used or the -a switch is used without a number following it, and the -r\n"
2299);
2300 printf(
2301" switch is not used. Regional attributes and area constraints are\n");
2302 printf(
2303" propagated in the same manner as holes: you specify a point for each\n");
2304 printf(
2305" attribute and/or constraint, and the attribute and/or constraint\n");
2306 printf(
2307" affects the whole region (bounded by segments) containing the point.\n");
2308 printf(
2309" If two values are written on a line after the x and y coordinate, the\n");
2310 printf(
2311" first such value is assumed to be a regional attribute (but is only\n");
2312 printf(
2313" applied if the -A switch is selected), and the second value is assumed\n"
2314);
2315 printf(
2316" to be a regional area constraint (but is only applied if the -a switch\n"
2317);
2318 printf(
2319" is selected). You may specify just one value after the coordinates,\n");
2320 printf(
2321" which can serve as both an attribute and an area constraint, depending\n"
2322);
2323 printf(
2324" on the choice of switches. If you are using the -A and -a switches\n");
2325 printf(
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");
2328 printf(
2329" When a triangulation is created from a .poly file, you must either\n");
2330 printf(
2331" enclose the entire region to be triangulated in PSLG segments, or\n");
2332 printf(
2333" use the -c switch, which automatically creates extra segments that\n");
2334 printf(
2335" enclose the convex hull of the PSLG. If you do not use the -c switch,\n"
2336);
2337 printf(
2338" Triangle eats all triangles that are not enclosed by segments; if you\n");
2339 printf(
2340" are not careful, your whole triangulation may be eaten away. If you do\n"
2341);
2342 printf(
2343" use the -c switch, you can still produce concavities by the appropriate\n"
2344);
2345 printf(
2346" placement of holes just inside the boundary of the convex hull.\n");
2347 printf("\n");
2348 printf(
2349" An ideal PSLG has no intersecting segments, nor any vertices that lie\n");
2350 printf(
2351" upon segments (except, of course, the endpoints of each segment). You\n"
2352);
2353 printf(
2354" aren't required to make your .poly files ideal, but you should be aware\n"
2355);
2356 printf(
2357" of what can go wrong. Segment intersections are relatively safe--\n");
2358 printf(
2359" Triangle calculates the intersection points for you and adds them to\n");
2360 printf(
2361" the triangulation--as long as your machine's floating-point precision\n");
2362 printf(
2363" doesn't become a problem. You are tempting the fates if you have three\n"
2364);
2365 printf(
2366" segments that cross at the same location, and expect Triangle to figure\n"
2367);
2368 printf(
2369" out where the intersection point is. Thanks to floating-point roundoff\n"
2370);
2371 printf(
2372" error, Triangle will probably decide that the three segments intersect\n"
2373);
2374 printf(
2375" at three different points, and you will find a minuscule triangle in\n");
2376 printf(
2377" your output--unless Triangle tries to refine the tiny triangle, uses\n");
2378 printf(
2379" up the last bit of machine precision, and fails to terminate at all.\n");
2380 printf(
2381" You're better off putting the intersection point in the input files,\n");
2382 printf(
2383" and manually breaking up each segment into two. Similarly, if you\n");
2384 printf(
2385" place a vertex at the middle of a segment, and hope that Triangle will\n"
2386);
2387 printf(
2388" break up the segment at that vertex, you might get lucky. On the other\n"
2389);
2390 printf(
2391" hand, Triangle might decide that the vertex doesn't lie precisely on\n");
2392 printf(
2393" the segment, and you'll have a needle-sharp triangle in your output--or\n"
2394);
2395 printf(" a lot of tiny triangles if you're generating a quality mesh.\n");
2396 printf("\n");
2397 printf(
2398" When Triangle reads a .poly file, it also writes a .poly file, which\n");
2399 printf(
2400" includes all the subsegments--the edges that are parts of input\n");
2401 printf(
2402" segments. If the -c switch is used, the output .poly file also\n");
2403 printf(
2404" includes all of the edges on the convex hull. Hence, the output .poly\n"
2405);
2406 printf(
2407" file is useful for finding edges associated with input segments and for\n"
2408);
2409 printf(
2410" setting boundary conditions in finite element simulations. Moreover,\n");
2411 printf(
2412" you will need the output .poly file if you plan to refine the output\n");
2413 printf(
2414" mesh, and don't want segments to be missing in later triangulations.\n");
2415 printf("\n");
2416 printf(" .area files:\n");
2417 printf(" First line: <# of triangles>\n");
2418 printf(" Following lines: <triangle #> <maximum area>\n");
2419 printf("\n");
2420 printf(
2421" An .area file associates with each triangle a maximum area that is used\n"
2422);
2423 printf(
2424" for mesh refinement. As with other file formats, every triangle must\n");
2425 printf(
2426" be represented, and the triangles must be numbered consecutively. A\n");
2427 printf(
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");
2432 printf(
2433" Following lines: <edge #> <endpoint> <endpoint> [boundary marker]\n");
2434 printf("\n");
2435 printf(
2436" Endpoints are indices into the corresponding .node file. Triangle can\n"
2437);
2438 printf(
2439" produce .edge files (use the -e switch), but cannot read them. The\n");
2440 printf(
2441" optional column of boundary markers is suppressed by the -B switch.\n");
2442 printf("\n");
2443 printf(
2444" In Voronoi diagrams, one also finds a special kind of edge that is an\n");
2445 printf(
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");
2449 printf(
2450" The `direction' is a floating-point vector that indicates the direction\n"
2451);
2452 printf(" of the infinite ray.\n\n");
2453 printf(" .neigh files:\n");
2454 printf(
2455" First line: <# of triangles> <# of neighbors per triangle (always 3)>\n"
2456);
2457 printf(
2458" Following lines: <triangle #> <neighbor> <neighbor> <neighbor>\n");
2459 printf("\n");
2460 printf(
2461" Neighbors are indices into the corresponding .ele file. An index of -1\n"
2462);
2463 printf(
2464" indicates no neighbor (because the triangle is on an exterior\n");
2465 printf(
2466" boundary). The first neighbor of triangle i is opposite the first\n");
2467 printf(" corner of triangle i, and so on.\n\n");
2468 printf(
2469" Triangle can produce .neigh files (use the -n switch), but cannot read\n"
2470);
2471 printf(" them.\n\n");
2472 printf("Boundary Markers:\n\n");
2473 printf(
2474" Boundary markers are tags used mainly to identify which output vertices\n");
2475 printf(
2476" and edges are associated with which PSLG segment, and to identify which\n");
2477 printf(
2478" vertices and edges occur on a boundary of the triangulation. A common\n");
2479 printf(
2480" use is to determine where boundary conditions should be applied to a\n");
2481 printf(
2482" finite element mesh. You can prevent boundary markers from being written\n"
2483);
2484 printf(" into files produced by Triangle by using the -B switch.\n\n");
2485 printf(
2486" The boundary marker associated with each segment in an output .poly file\n"
2487);
2488 printf(" and each edge in an output .edge file is chosen as follows:\n");
2489 printf(
2490" - If an output edge is part or all of a PSLG segment with a nonzero\n");
2491 printf(
2492" boundary marker, then the edge is assigned the same marker.\n");
2493 printf(
2494" - Otherwise, if the edge lies on a boundary of the triangulation\n");
2495 printf(
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");
2499 printf(
2500" The boundary marker associated with each vertex in an output .node file\n");
2501 printf(" is chosen as follows:\n");
2502 printf(
2503" - If a vertex is assigned a nonzero boundary marker in the input file,\n"
2504);
2505 printf(
2506" then it is assigned the same marker in the output .node file.\n");
2507 printf(
2508" - Otherwise, if the vertex lies on a PSLG segment (even if it is an\n");
2509 printf(
2510" endpoint of the segment) with a nonzero boundary marker, then the\n");
2511 printf(
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");
2514 printf(
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");
2518 printf("\n");
2519 printf(
2520" If you want Triangle to determine for you which vertices and edges are on\n"
2521);
2522 printf(
2523" the boundary, assign them the boundary marker zero (or use no markers at\n"
2524);
2525 printf(
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");
2529 printf(
2530" Because Triangle can read and refine its own triangulations, input\n");
2531 printf(
2532" and output files have iteration numbers. For instance, Triangle might\n");
2533 printf(
2534" read the files mesh.3.node, mesh.3.ele, and mesh.3.poly, refine the\n");
2535 printf(
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");
2538 printf(
2539" their iteration number is zero; hence, Triangle might read the file\n");
2540 printf(
2541" points.node, triangulate it, and produce the files points.1.node and\n");
2542 printf(" points.1.ele.\n\n");
2543 printf(
2544" Iteration numbers allow you to create a sequence of successively finer\n");
2545 printf(
2546" meshes suitable for multigrid methods. They also allow you to produce a\n"
2547);
2548 printf(
2549" sequence of meshes using error estimate-driven mesh refinement.\n");
2550 printf("\n");
2551 printf(
2552" If you're not using refinement or quality meshing, and you don't like\n");
2553 printf(
2554" iteration numbers, use the -I switch to disable them. This switch also\n");
2555 printf(
2556" disables output of .node and .poly files to prevent your input files from\n"
2557);
2558 printf(
2559" being overwritten. (If the input is a .poly file that contains its own\n");
2560 printf(
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");
2564 printf(
2565" `triangle dots' reads vertices from dots.node, and writes their Delaunay\n"
2566);
2567 printf(
2568" triangulation to dots.1.node and dots.1.ele. (dots.1.node is identical\n");
2569 printf(
2570" to dots.node.) `triangle -I dots' writes the triangulation to dots.ele\n");
2571 printf(
2572" instead. (No additional .node file is needed, so none is written.)\n");
2573 printf("\n");
2574 printf(
2575" `triangle -pe object.1' reads a PSLG from object.1.poly (and possibly\n");
2576 printf(
2577" object.1.node, if the vertices are omitted from object.1.poly) and writes\n"
2578);
2579 printf(
2580" its constrained Delaunay triangulation to object.2.node and object.2.ele.\n"
2581);
2582 printf(
2583" The segments are copied to object.2.poly, and all edges are written to\n");
2584 printf(" object.2.edge.\n\n");
2585 printf(
2586" `triangle -pq31.5a.1 object' reads a PSLG from object.poly (and possibly\n"
2587);
2588 printf(
2589" object.node), generates a mesh whose angles are all between 31.5 and 117\n"
2590);
2591 printf(
2592" degrees and whose triangles all have areas of 0.1 or less, and writes the\n"
2593);
2594 printf(
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");
2597 printf("\n");
2598 printf(
2599" Here is a sample file `box.poly' describing a square with a square hole:\n"
2600);
2601 printf("\n");
2602 printf(
2603" # A box with eight vertices in 2D, no attributes, one boundary marker.\n"
2604);
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");
2617 printf(" 5 1\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");
2625 printf(" 1\n");
2626 printf(" 1 1.5 1.5\n");
2627 printf("\n");
2628 printf(
2629" Note that some segments are missing from the outer square, so you must\n");
2630 printf(
2631" use the `-c' switch. After `triangle -pqc box.poly', here is the output\n"
2632);
2633 printf(
2634" file `box.1.node', with twelve vertices. The last four vertices were\n");
2635 printf(
2636" added to meet the angle constraint. Vertices 1, 2, and 9 have markers\n");
2637 printf(
2638" from segment 1. Vertices 6 and 8 have markers from segment 4. All the\n");
2639 printf(
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");
2656 printf("\n");
2657 printf(" Here is the output file `box.1.ele', with twelve triangles.\n");
2658 printf("\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");
2673 printf(
2674" Here is the output file `box.1.poly'. Note that segments have been added\n"
2675);
2676 printf(
2677" to represent the convex hull, and some segments have been subdivided by\n");
2678 printf(
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");
2681 printf("\n");
2682 printf(" 0 2 0 1\n");
2683 printf(" 12 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");
2696 printf(" 1\n");
2697 printf(" 1 1.5 1.5\n");
2698 printf(" # Generated by triangle -pqc box.poly\n");
2699 printf("\n");
2700 printf("Refinement and Area Constraints:\n");
2701 printf("\n");
2702 printf(
2703" The -r switch causes a mesh (.node and .ele files) to be read and\n");
2704 printf(
2705" refined. If the -p switch is also used, a .poly file is read and used to\n"
2706);
2707 printf(
2708" specify edges that are constrained and cannot be eliminated (although\n");
2709 printf(
2710" they can be subdivided into smaller edges) by the refinement process.\n");
2711 printf("\n");
2712 printf(
2713" When you refine a mesh, you generally want to impose tighter constraints.\n"
2714);
2715 printf(
2716" One way to accomplish this is to use -q with a larger angle, or -a\n");
2717 printf(
2718" followed by a smaller area than you used to generate the mesh you are\n");
2719 printf(
2720" refining. Another way to do this is to create an .area file, which\n");
2721 printf(
2722" specifies a maximum area for each triangle, and use the -a switch\n");
2723 printf(
2724" (without a number following). Each triangle's area constraint is applied\n"
2725);
2726 printf(
2727" to that triangle. Area constraints tend to diffuse as the mesh is\n");
2728 printf(
2729" refined, so if there are large variations in area constraint between\n");
2730 printf(
2731" adjacent triangles, you may not get the results you want. In that case,\n"
2732);
2733 printf(
2734" consider instead using the -u switch and writing a C procedure that\n");
2735 printf(" determines which triangles are too large.\n\n");
2736 printf(
2737" If you are refining a mesh composed of linear (three-node) elements, the\n"
2738);
2739 printf(
2740" output mesh contains all the nodes present in the input mesh, in the same\n"
2741);
2742 printf(
2743" order, with new nodes added at the end of the .node file. However, the\n");
2744 printf(
2745" refinement is not hierarchical: there is no guarantee that each output\n");
2746 printf(
2747" element is contained in a single input element. Often, an output element\n"
2748);
2749 printf(
2750" can overlap two or three input elements, and some input edges are not\n");
2751 printf(
2752" present in the output mesh. Hence, a sequence of refined meshes forms a\n"
2753);
2754 printf(
2755" hierarchy of nodes, but not a hierarchy of elements. If you refine a\n");
2756 printf(
2757" mesh of higher-order elements, the hierarchical property applies only to\n"
2758);
2759 printf(
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");
2762 printf(
2763" Maximum area constraints in .poly files operate differently from those in\n"
2764);
2765 printf(
2766" .area files. A maximum area in a .poly file applies to the whole\n");
2767 printf(
2768" (segment-bounded) region in which a point falls, whereas a maximum area\n");
2769 printf(
2770" in an .area file applies to only one triangle. Area constraints in .poly\n"
2771);
2772 printf(
2773" files are used only when a mesh is first generated, whereas area\n");
2774 printf(
2775" constraints in .area files are used only to refine an existing mesh, and\n"
2776);
2777 printf(
2778" are typically based on a posteriori error estimates resulting from a\n");
2779 printf(" finite element simulation on that mesh.\n\n");
2780 printf(
2781" `triangle -rq25 object.1' reads object.1.node and object.1.ele, then\n");
2782 printf(
2783" refines the triangulation to enforce a 25 degree minimum angle, and then\n"
2784);
2785 printf(
2786" writes the refined triangulation to object.2.node and object.2.ele.\n");
2787 printf("\n");
2788 printf(
2789" `triangle -rpaa6.2 z.3' reads z.3.node, z.3.ele, z.3.poly, and z.3.area.\n"
2790);
2791 printf(
2792" After reconstructing the mesh and its subsegments, Triangle refines the\n");
2793 printf(
2794" mesh so that no triangle has area greater than 6.2, and furthermore the\n");
2795 printf(
2796" triangles satisfy the maximum area constraints in z.3.area. No angle\n");
2797 printf(
2798" bound is imposed at all. The output is written to z.4.node, z.4.ele, and\n"
2799);
2800 printf(" z.4.poly.\n\n");
2801 printf(
2802" The sequence `triangle -qa1 x', `triangle -rqa.3 x.1', `triangle -rqa.1\n");
2803 printf(
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");
2807 printf(
2808" If the input is a vertex set (not a PSLG), Triangle produces its convex\n");
2809 printf(
2810" hull as a by-product in the output .poly file if you use the -c switch.\n");
2811 printf(
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");
2814 printf(
2815" If the input is an unconstrained mesh (you are using the -r switch but\n");
2816 printf(
2817" not the -p switch), Triangle produces a list of its boundary edges\n");
2818 printf(
2819" (including hole boundaries) as a by-product when you use the -c switch.\n");
2820 printf(
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");
2824 printf(
2825" The -v switch produces a Voronoi diagram, in files suffixed .v.node and\n");
2826 printf(
2827" .v.edge. For example, `triangle -v points' reads points.node, produces\n");
2828 printf(
2829" its Delaunay triangulation in points.1.node and points.1.ele, and\n");
2830 printf(
2831" produces its Voronoi diagram in points.1.v.node and points.1.v.edge. The\n"
2832);
2833 printf(
2834" .v.node file contains a list of all Voronoi vertices, and the .v.edge\n");
2835 printf(
2836" file contains a list of all Voronoi edges, some of which may be infinite\n"
2837);
2838 printf(
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");
2841 printf(
2842" This implementation does not use exact arithmetic to compute the Voronoi\n"
2843);
2844 printf(
2845" vertices, and does not check whether neighboring vertices are identical.\n"
2846);
2847 printf(
2848" Be forewarned that if the Delaunay triangulation is degenerate or\n");
2849 printf(
2850" near-degenerate, the Voronoi diagram may have duplicate vertices or\n");
2851 printf(" crossing edges.\n\n");
2852 printf(
2853" The result is a valid Voronoi diagram only if Triangle's output is a true\n"
2854);
2855 printf(
2856" Delaunay triangulation. The Voronoi output is usually meaningless (and\n");
2857 printf(
2858" may contain crossing edges and other pathology) if the output is a CDT or\n"
2859);
2860 printf(
2861" CCDT, or if it has holes or concavities. If the triangulated domain is\n");
2862 printf(
2863" convex and has no holes, you can use -D switch to force Triangle to\n");
2864 printf(
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");
2868 printf(
2869" You may wish to know which triangles are adjacent to a certain Delaunay\n");
2870 printf(
2871" edge in an .edge file, which Voronoi cells are adjacent to a certain\n");
2872 printf(
2873" Voronoi edge in a .v.edge file, or which Voronoi cells are adjacent to\n");
2874 printf(
2875" each other. All of this information can be found by cross-referencing\n");
2876 printf(
2877" output files with the recollection that the Delaunay triangulation and\n");
2878 printf(" the Voronoi diagram are planar duals.\n\n");
2879 printf(
2880" Specifically, edge i of an .edge file is the dual of Voronoi edge i of\n");
2881 printf(
2882" the corresponding .v.edge file, and is rotated 90 degrees counterclock-\n");
2883 printf(
2884" wise from the Voronoi edge. Triangle j of an .ele file is the dual of\n");
2885 printf(
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");
2888 printf(
2889" Hence, to find the triangles adjacent to a Delaunay edge, look at the\n");
2890 printf(
2891" vertices of the corresponding Voronoi edge. If the endpoints of a\n");
2892 printf(
2893" Voronoi edge are Voronoi vertices 2 and 6 respectively, then triangles 2\n"
2894);
2895 printf(
2896" and 6 adjoin the left and right sides of the corresponding Delaunay edge,\n"
2897);
2898 printf(
2899" respectively. To find the Voronoi cells adjacent to a Voronoi edge, look\n"
2900);
2901 printf(
2902" at the endpoints of the corresponding Delaunay edge. If the endpoints of\n"
2903);
2904 printf(
2905" a Delaunay edge are input vertices 7 and 12, then Voronoi cells 7 and 12\n"
2906);
2907 printf(
2908" adjoin the right and left sides of the corresponding Voronoi edge,\n");
2909 printf(
2910" respectively. To find which Voronoi cells are adjacent to each other,\n");
2911 printf(" just read the list of Delaunay edges.\n\n");
2912 printf(
2913" Triangle does not write a list of the edges adjoining each Voronoi cell,\n"
2914);
2915 printf(
2916" but you can reconstructed it straightforwardly. For instance, to find\n");
2917 printf(
2918" all the edges of Voronoi cell 1, search the output .edge file for every\n");
2919 printf(
2920" edge that has input vertex 1 as an endpoint. The corresponding dual\n");
2921 printf(
2922" edges in the output .v.edge file form the boundary of Voronoi cell 1.\n");
2923 printf("\n");
2924 printf(
2925" For each Voronoi vertex, the .neigh file gives a list of the three\n");
2926 printf(
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");
2930 printf(
2931" Triangle generates meshes with subparametric quadratic elements if the\n");
2932 printf(
2933" -o2 switch is specified. Quadratic elements have six nodes per element,\n"
2934);
2935 printf(
2936" rather than three. `Subparametric' means that the edges of the triangles\n"
2937);
2938 printf(
2939" are always straight, so that subparametric quadratic elements are\n");
2940 printf(
2941" geometrically identical to linear elements, even though they can be used\n"
2942);
2943 printf(
2944" with quadratic interpolating functions. The three extra nodes of an\n");
2945 printf(
2946" element fall at the midpoints of the three edges, with the fourth, fifth,\n"
2947);
2948 printf(
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");
2952 printf(
2953" If two input segments adjoin each other at a small angle, clearly the -q\n"
2954);
2955 printf(
2956" switch cannot remove the small angle. Moreover, Triangle may have no\n");
2957 printf(
2958" choice but to generate additional triangles whose smallest angles are\n");
2959 printf(
2960" smaller than the specified bound. However, these triangles only appear\n");
2961 printf(
2962" between input segments separated by small angles. Moreover, if you\n");
2963 printf(
2964" request a minimum angle of theta degrees, Triangle will generally produce\n"
2965);
2966 printf(
2967" no angle larger than 180 - 2 theta, even if it is forced to compromise on\n"
2968);
2969 printf(" the minimum angle.\n\n");
2970 printf("Statistics:\n\n");
2971 printf(
2972" After generating a mesh, Triangle prints a count of entities in the\n");
2973 printf(
2974" output mesh, including the number of vertices, triangles, edges, exterior\n"
2975);
2976 printf(
2977" boundary edges (i.e. subsegments on the boundary of the triangulation,\n");
2978 printf(
2979" including hole boundaries), interior boundary edges (i.e. subsegments of\n"
2980);
2981 printf(
2982" input segments not on the boundary), and total subsegments. If you've\n");
2983 printf(
2984" forgotten the statistics for an existing mesh, run Triangle on that mesh\n"
2985);
2986 printf(
2987" with the -rNEP switches to read the mesh and print the statistics without\n"
2988);
2989 printf(
2990" writing any files. Use -rpNEP if you've got a .poly file for the mesh.\n");
2991 printf("\n");
2992 printf(
2993" The -V switch produces extended statistics, including a rough estimate\n");
2994 printf(
2995" of memory use, the number of calls to geometric predicates, and\n");
2996 printf(
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");
3000 printf(
3001" Triangle uses adaptive exact arithmetic to perform what computational\n");
3002 printf(
3003" geometers call the `orientation' and `incircle' tests. If the floating-\n"
3004);
3005 printf(
3006" point arithmetic of your machine conforms to the IEEE 754 standard (as\n");
3007 printf(
3008" most workstations do), and does not use extended precision internal\n");
3009 printf(
3010" floating-point registers, then your output is guaranteed to be an\n");
3011 printf(
3012" absolutely true Delaunay or constrained Delaunay triangulation, roundoff\n"
3013);
3014 printf(
3015" error notwithstanding. The word `adaptive' implies that these arithmetic\n"
3016);
3017 printf(
3018" routines compute the result only to the precision necessary to guarantee\n"
3019);
3020 printf(
3021" correctness, so they are usually nearly as fast as their approximate\n");
3022 printf(" counterparts.\n\n");
3023 printf(
3024" May CPUs, including Intel x86 processors, have extended precision\n");
3025 printf(
3026" floating-point registers. These must be reconfigured so their precision\n"
3027);
3028 printf(
3029" is reduced to memory precision. Triangle does this if it is compiled\n");
3030 printf(" correctly. See the makefile for details.\n\n");
3031 printf(
3032" The exact tests can be disabled with the -X switch. On most inputs, this\n"
3033);
3034 printf(
3035" switch reduces the computation time by about eight percent--it's not\n");
3036 printf(
3037" worth the risk. There are rare difficult inputs (having many collinear\n");
3038 printf(
3039" and cocircular vertices), however, for which the difference in speed\n");
3040 printf(
3041" could be a factor of two. Be forewarned that these are precisely the\n");
3042 printf(
3043" inputs most likely to cause errors if you use the -X switch. Hence, the\n"
3044);
3045 printf(" -X switch is not recommended.\n\n");
3046 printf(
3047" Unfortunately, the exact tests don't solve every numerical problem.\n");
3048 printf(
3049" Exact arithmetic is not used to compute the positions of new vertices,\n");
3050 printf(
3051" because the bit complexity of vertex coordinates would grow without\n");
3052 printf(
3053" bound. Hence, segment intersections aren't computed exactly; in very\n");
3054 printf(
3055" unusual cases, roundoff error in computing an intersection point might\n");
3056 printf(
3057" actually lead to an inverted triangle and an invalid triangulation.\n");
3058 printf(
3059" (This is one reason to specify your own intersection points in your .poly\n"
3060);
3061 printf(
3062" files.) Similarly, exact arithmetic is not used to compute the vertices\n"
3063);
3064 printf(" of the Voronoi diagram.\n\n");
3065 printf(
3066" Another pair of problems not solved by the exact arithmetic routines is\n");
3067 printf(
3068" underflow and overflow. If Triangle is compiled for double precision\n");
3069 printf(
3070" arithmetic, I believe that Triangle's geometric predicates work correctly\n"
3071);
3072 printf(
3073" if the exponent of every input coordinate falls in the range [-148, 201].\n"
3074);
3075 printf(
3076" Underflow can silently prevent the orientation and incircle tests from\n");
3077 printf(
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");
3085 printf(
3086" If you're using a PSLG, you've probably failed to specify a proper set\n"
3087);
3088 printf(
3089" of bounding segments, or forgotten to use the -c switch. Or you may\n");
3090 printf(
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");
3093 printf(
3094" Alternatively, all your input vertices may be collinear, in which case\n"
3095);
3096 printf(" you can hardly expect to triangulate them.\n\n");
3097 printf(" `Triangle doesn't terminate, or just crashes.'\n\n");
3098 printf(
3099" Bad things can happen when triangles get so small that the distance\n");
3100 printf(
3101" between their vertices isn't much larger than the precision of your\n");
3102 printf(
3103" machine's arithmetic. If you've compiled Triangle for single-precision\n"
3104);
3105 printf(
3106" arithmetic, you might do better by recompiling it for double-precision.\n"
3107);
3108 printf(
3109" Then again, you might just have to settle for more lenient constraints\n"
3110);
3111 printf(
3112" on the minimum angle and the maximum area than you had planned.\n");
3113 printf("\n");
3114 printf(
3115" You can minimize precision problems by ensuring that the origin lies\n");
3116 printf(
3117" inside your vertex set, or even inside the densest part of your\n");
3118 printf(
3119" mesh. If you're triangulating an object whose x-coordinates all fall\n");
3120 printf(
3121" between 6247133 and 6247134, you're not leaving much floating-point\n");
3122 printf(" precision for Triangle to work with.\n\n");
3123 printf(
3124" Precision problems can occur covertly if the input PSLG contains two\n");
3125 printf(
3126" segments that meet (or intersect) at an extremely small angle, or if\n");
3127 printf(
3128" such an angle is introduced by the -c switch. If you don't realize\n");
3129 printf(
3130" that a tiny angle is being formed, you might never discover why\n");
3131 printf(
3132" Triangle is crashing. To check for this possibility, use the -S switch\n"
3133);
3134 printf(
3135" (with an appropriate limit on the number of Steiner points, found by\n");
3136 printf(
3137" trial-and-error) to stop Triangle early, and view the output .poly file\n"
3138);
3139 printf(
3140" with Show Me (described below). Look carefully for regions where dense\n"
3141);
3142 printf(
3143" clusters of vertices are forming and for small angles between segments.\n"
3144);
3145 printf(
3146" Zoom in closely, as such segments might look like a single segment from\n"
3147);
3148 printf(" a distance.\n\n");
3149 printf(
3150" If some of the input values are too large, Triangle may suffer a\n");
3151 printf(
3152" floating exception due to overflow when attempting to perform an\n");
3153 printf(
3154" orientation or incircle test. (Read the section on exact arithmetic\n");
3155 printf(
3156" above.) Again, I recommend compiling Triangle for double (rather\n");
3157 printf(" than single) precision arithmetic.\n\n");
3158 printf(
3159" Unexpected problems can arise if you use quality meshing (-q, -a, or\n");
3160 printf(
3161" -u) with an input that is not segment-bounded--that is, if your input\n");
3162 printf(
3163" is a vertex set, or you're using the -c switch. If the convex hull of\n"
3164);
3165 printf(
3166" your input vertices has collinear vertices on its boundary, an input\n");
3167 printf(
3168" vertex that you think lies on the convex hull might actually lie just\n");
3169 printf(
3170" inside the convex hull. If so, the vertex and the nearby convex hull\n");
3171 printf(
3172" edge form an extremely thin triangle. When Triangle tries to refine\n");
3173 printf(
3174" the mesh to enforce angle and area constraints, Triangle might generate\n"
3175);
3176 printf(
3177" extremely tiny triangles, or it might fail because of insufficient\n");
3178 printf(" floating-point precision.\n\n");
3179 printf(
3180" `The numbering of the output vertices doesn't match the input vertices.'\n"
3181);
3182 printf("\n");
3183 printf(
3184" You may have had duplicate input vertices, or you may have eaten some\n");
3185 printf(
3186" of your input vertices with a hole, or by placing them outside the area\n"
3187);
3188 printf(
3189" enclosed by segments. In any case, you can solve the problem by not\n");
3190 printf(" using the -j switch.\n\n");
3191 printf(
3192" `Triangle executes without incident, but when I look at the resulting\n");
3193 printf(
3194" mesh, it has overlapping triangles or other geometric inconsistencies.'\n");
3195 printf("\n");
3196 printf(
3197" If you select the -X switch, Triangle occasionally makes mistakes due\n");
3198 printf(
3199" to floating-point roundoff error. Although these errors are rare,\n");
3200 printf(
3201" don't use the -X switch. If you still have problems, please report the\n"
3202);
3203 printf(" bug.\n\n");
3204 printf(
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");
3208 printf("\n");
3209 printf(
3210" If your input is a PSLG (-p), you can only expect a meaningful Voronoi\n"
3211);
3212 printf(
3213" diagram if the domain you are triangulating is convex and free of\n");
3214 printf(
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");
3217 printf(
3218" Strange things can happen if you've taken liberties with your PSLG. Do\n");
3219 printf(
3220" you have a vertex lying in the middle of a segment? Triangle sometimes\n");
3221 printf(
3222" copes poorly with that sort of thing. Do you want to lay out a collinear\n"
3223);
3224 printf(
3225" row of evenly spaced, segment-connected vertices? Have you simply\n");
3226 printf(
3227" defined one long segment connecting the leftmost vertex to the rightmost\n"
3228);
3229 printf(
3230" vertex, and a bunch of vertices lying along it? This method occasionally\n"
3231);
3232 printf(
3233" works, especially with horizontal and vertical lines, but often it\n");
3234 printf(
3235" doesn't, and you'll have to connect each adjacent pair of vertices with a\n"
3236);
3237 printf(" separate segment. If you don't like it, tough.\n\n");
3238 printf(
3239" Furthermore, if you have segments that intersect other than at their\n");
3240 printf(
3241" endpoints, try not to let the intersections fall extremely close to PSLG\n"
3242);
3243 printf(" vertices or each other.\n\n");
3244 printf(
3245" If you have problems refining a triangulation not produced by Triangle:\n");
3246 printf(
3247" Are you sure the triangulation is geometrically valid? Is it formatted\n");
3248 printf(
3249" correctly for Triangle? Are the triangles all listed so the first three\n"
3250);
3251 printf(
3252" vertices are their corners in counterclockwise order? Are all of the\n");
3253 printf(
3254" triangles constrained Delaunay? Triangle's Delaunay refinement algorithm\n"
3255);
3256 printf(" assumes that it starts with a CDT.\n\n");
3257 printf("Show Me:\n\n");
3258 printf(
3259" Triangle comes with a separate program named `Show Me', whose primary\n");
3260 printf(
3261" purpose is to draw meshes on your screen or in PostScript. Its secondary\n"
3262);
3263 printf(
3264" purpose is to check the validity of your input files, and do so more\n");
3265 printf(
3266" thoroughly than Triangle does. Unlike Triangle, Show Me requires that\n");
3267 printf(
3268" you have the X Windows system. Sorry, Microsoft Windows users.\n");
3269 printf("\n");
3270 printf("Triangle on the Web:\n");
3271 printf("\n");
3272 printf(" To see an illustrated version of these instructions, check out\n");
3273 printf("\n");
3274 printf(" http://www.cs.cmu.edu/~quake/triangle.html\n");
3275 printf("\n");
3276 printf("A Brief Plea:\n");
3277 printf("\n");
3278 printf(
3279" If you use Triangle, and especially if you use it to accomplish real\n");
3280 printf(
3281" work, I would like very much to hear from you. A short letter or email\n");
3282 printf(
3283" (to jrs@cs.berkeley.edu) describing how you use Triangle will mean a lot\n"
3284);
3285 printf(
3286" to me. The more people I know are using this program, the more easily I\n"
3287);
3288 printf(
3289" can justify spending time on improvements, which in turn will benefit\n");
3290 printf(
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");
3293 printf(
3294" If you use a mesh generated by Triangle in a publication, please include\n"
3295);
3296 printf(
3297" an acknowledgment as well. And please spell Triangle with a capital `T'!\n"
3298);
3299 printf(
3300" If you want to include a citation, use `Jonathan Richard Shewchuk,\n");
3301 printf(
3302" ``Triangle: Engineering a 2D Quality Mesh Generator and Delaunay\n");
3303 printf(
3304" Triangulator,'' in Applied Computational Geometry: Towards Geometric\n");
3305 printf(
3306" Engineering (Ming C. Lin and Dinesh Manocha, editors), volume 1148 of\n");
3307 printf(
3308" Lecture Notes in Computer Science, pages 203-222, Springer-Verlag,\n");
3309 printf(
3310" Berlin, May 1996. (From the First ACM Workshop on Applied Computational\n"
3311);
3312 printf(" Geometry.)'\n\n");
3313 printf("Research credit:\n\n");
3314 printf(
3315" Of course, I can take credit for only a fraction of the ideas that made\n");
3316 printf(
3317" this mesh generator possible. Triangle owes its existence to the efforts\n"
3318);
3319 printf(
3320" of many fine computational geometers and other researchers, including\n");
3321 printf(
3322" Marshall Bern, L. Paul Chew, Kenneth L. Clarkson, Boris Delaunay, Rex A.\n"
3323);
3324 printf(
3325" Dwyer, David Eppstein, Steven Fortune, Leonidas J. Guibas, Donald E.\n");
3326 printf(
3327" Knuth, Charles L. Lawson, Der-Tsai Lee, Gary L. Miller, Ernst P. Mucke,\n");
3328 printf(
3329" Steven E. Pav, Douglas M. Priest, Jim Ruppert, Isaac Saias, Bruce J.\n");
3330 printf(
3331" Schachter, Micha Sharir, Peter W. Shor, Daniel D. Sleator, Jorge Stolfi,\n"
3332);
3333 printf(" Robert E. Tarjan, Alper Ungor, Christopher J. Van Wyk, Noel J.\n");
3334 printf(
3335" Walkington, and Binhai Zhu. See the comments at the beginning of the\n");
3336 printf(" source code for references.\n\n");
3337 triexit(0);
3338}
3339
3340#endif /* not TRILIBRARY */
3341
3342/*****************************************************************************/
3343/* */
3344/* internalerror() Ask the user to send me the defective product. Exit. */
3345/* */
3346/*****************************************************************************/
3347
3348void internalerror()
3349{
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");
3353 triexit(1);
3354}
3355
3356/*****************************************************************************/
3357/* */
3358/* parsecommandline() Read the command line, identify switches, and set */
3359/* up options and file names. */
3360/* */
3361/*****************************************************************************/
3362
3363#ifdef ANSI_DECLARATORS
3364void parsecommandline(int argc, char **argv, struct behavior *b)
3365#else /* not ANSI_DECLARATORS */
3366void parsecommandline(argc, argv, b)
3367int argc;
3368char **argv;
3369struct behavior *b;
3370#endif /* not ANSI_DECLARATORS */
3371
3372{
3373#ifdef TRILIBRARY
3374#define STARTINDEX 0
3375#else /* not TRILIBRARY */
3376#define STARTINDEX 1
3377 int increment;
3378 int meshnumber;
3379#endif /* not TRILIBRARY */
3380 int i, j, k;
3381 char workstring[FILENAMESIZE];
3382
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;
3386 b->firstnumber = 1;
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;
3392 b->dwyer = 1;
3393 b->splitseg = 0;
3394 b->docheck = 0;
3395 b->nobisect = 0;
3396 b->conformdel = 0;
3397 b->steiner = -1;
3398 b->order = 1;
3399 b->minangle = 0.0;
3400 b->maxarea = -1.0;
3401 b->quiet = b->verbose = 0;
3402#ifndef TRILIBRARY
3403 b->innodefilename[0] = '\0';
3404#endif /* not TRILIBRARY */
3405
3406 for (i = STARTINDEX; i < argc; i++) {
3407#ifndef TRILIBRARY
3408 if (argv[i][0] == '-') {
3409#endif /* not TRILIBRARY */
3410 for (j = STARTINDEX; argv[i][j] != '\0'; j++) {
3411 if (argv[i][j] == 'p') {
3412 b->poly = 1;
3413 }
3414#ifndef CDT_ONLY
3415 if (argv[i][j] == 'r') {
3416 b->refine = 1;
3417 }
3418 if (argv[i][j] == 'q') {
3419 b->quality = 1;
3420 if (((argv[i][j + 1] >= '0') && (argv[i][j + 1] <= '9')) ||
3421 (argv[i][j + 1] == '.')) {
3422 k = 0;
3423 while (((argv[i][j + 1] >= '0') && (argv[i][j + 1] <= '9')) ||
3424 (argv[i][j + 1] == '.')) {
3425 j++;
3426 workstring[k] = argv[i][j];
3427 k++;
3428 }
3429 workstring[k] = '\0';
3430 b->minangle = (REAL) strtod(workstring, (char **) NULL);
3431 } else {
3432 b->minangle = 20.0;
3433 }
3434 }
3435 if (argv[i][j] == 'a') {
3436 b->quality = 1;
3437 if (((argv[i][j + 1] >= '0') && (argv[i][j + 1] <= '9')) ||
3438 (argv[i][j + 1] == '.')) {
3439 b->fixedarea = 1;
3440 k = 0;
3441 while (((argv[i][j + 1] >= '0') && (argv[i][j + 1] <= '9')) ||
3442 (argv[i][j + 1] == '.')) {
3443 j++;
3444 workstring[k] = argv[i][j];
3445 k++;
3446 }
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");
3451 triexit(1);
3452 }
3453 } else {
3454 b->vararea = 1;
3455 }
3456 }
3457 if (argv[i][j] == 'u') {
3458 b->quality = 1;
3459 b->usertest = 1;
3460 }
3461#endif /* not CDT_ONLY */
3462 if (argv[i][j] == 'A') {
3463 b->regionattrib = 1;
3464 }
3465 if (argv[i][j] == 'c') {
3466 b->convex = 1;
3467 }
3468 if (argv[i][j] == 'w') {
3469 b->weighted = 1;
3470 }
3471 if (argv[i][j] == 'W') {
3472 b->weighted = 2;
3473 }
3474 if (argv[i][j] == 'j') {
3475 b->jettison = 1;
3476 }
3477 if (argv[i][j] == 'z') {
3478 b->firstnumber = 0;
3479 }
3480 if (argv[i][j] == 'e') {
3481 b->edgesout = 1;
3482 }
3483 if (argv[i][j] == 'v') {
3484 b->voronoi = 1;
3485 }
3486 if (argv[i][j] == 'n') {
3487 b->neighbors = 1;
3488 }
3489 if (argv[i][j] == 'g') {
3490 b->geomview = 1;
3491 }
3492 if (argv[i][j] == 'B') {
3493 b->nobound = 1;
3494 }
3495 if (argv[i][j] == 'P') {
3496 b->nopolywritten = 1;
3497 }
3498 if (argv[i][j] == 'N') {
3499 b->nonodewritten = 1;
3500 }
3501 if (argv[i][j] == 'E') {
3502 b->noelewritten = 1;
3503 }
3504#ifndef TRILIBRARY
3505 if (argv[i][j] == 'I') {
3506 b->noiterationnum = 1;
3507 }
3508#endif /* not TRILIBRARY */
3509 if (argv[i][j] == 'O') {
3510 b->noholes = 1;
3511 }
3512 if (argv[i][j] == 'X') {
3513 b->noexact = 1;
3514 }
3515 if (argv[i][j] == 'o') {
3516 if (argv[i][j + 1] == '2') {
3517 j++;
3518 b->order = 2;
3519 }
3520 }
3521#ifndef CDT_ONLY
3522 if (argv[i][j] == 'Y') {
3523 b->nobisect++;
3524 }
3525 if (argv[i][j] == 'S') {
3526 b->steiner = 0;
3527 while ((argv[i][j + 1] >= '0') && (argv[i][j + 1] <= '9')) {
3528 j++;
3529 b->steiner = b->steiner * 10 + (int) (argv[i][j] - '0');
3530 }
3531 }
3532#endif /* not CDT_ONLY */
3533#ifndef REDUCED
3534 if (argv[i][j] == 'i') {
3535 b->incremental = 1;
3536 }
3537 if (argv[i][j] == 'F') {
3538 b->sweepline = 1;
3539 }
3540#endif /* not REDUCED */
3541 if (argv[i][j] == 'l') {
3542 b->dwyer = 0;
3543 }
3544#ifndef REDUCED
3545#ifndef CDT_ONLY
3546 if (argv[i][j] == 's') {
3547 b->splitseg = 1;
3548 }
3549 if ((argv[i][j] == 'D') || (argv[i][j] == 'L')) {
3550 b->quality = 1;
3551 b->conformdel = 1;
3552 }
3553#endif /* not CDT_ONLY */
3554 if (argv[i][j] == 'C') {
3555 b->docheck = 1;
3556 }
3557#endif /* not REDUCED */
3558 if (argv[i][j] == 'Q') {
3559 b->quiet = 1;
3560 }
3561 if (argv[i][j] == 'V') {
3562 b->verbose++;
3563 }
3564#ifndef TRILIBRARY
3565 if ((argv[i][j] == 'h') || (argv[i][j] == 'H') ||
3566 (argv[i][j] == '?')) {
3567 info();
3568 }
3569#endif /* not TRILIBRARY */
3570 }
3571#ifndef TRILIBRARY
3572 } else {
3573 strncpy(b->innodefilename, argv[i], FILENAMESIZE - 1);
3574 b->innodefilename[FILENAMESIZE - 1] = '\0';
3575 }
3576#endif /* not TRILIBRARY */
3577 }
3578#ifndef TRILIBRARY
3579 if (b->innodefilename[0] == '\0') {
3580 syntax();
3581 }
3582 if (!strcmp(&b->innodefilename[strlen(b->innodefilename) - 5], ".node")) {
3583 b->innodefilename[strlen(b->innodefilename) - 5] = '\0';
3584 }
3585 if (!strcmp(&b->innodefilename[strlen(b->innodefilename) - 5], ".poly")) {
3586 b->innodefilename[strlen(b->innodefilename) - 5] = '\0';
3587 b->poly = 1;
3588 }
3589#ifndef CDT_ONLY
3590 if (!strcmp(&b->innodefilename[strlen(b->innodefilename) - 4], ".ele")) {
3591 b->innodefilename[strlen(b->innodefilename) - 4] = '\0';
3592 b->refine = 1;
3593 }
3594 if (!strcmp(&b->innodefilename[strlen(b->innodefilename) - 5], ".area")) {
3595 b->innodefilename[strlen(b->innodefilename) - 5] = '\0';
3596 b->refine = 1;
3597 b->quality = 1;
3598 b->vararea = 1;
3599 }
3600#endif /* not CDT_ONLY */
3601#endif /* not TRILIBRARY */
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;
3606 } else {
3607 b->offconstant = 0.475 * sqrt((1.0 + b->goodangle) / (1.0 - b->goodangle));
3608 }
3609 b->goodangle *= b->goodangle;
3610 if (b->refine && b->noiterationnum) {
3611 printf(
3612 "Error: You cannot use the -I switch when refining a triangulation.\n");
3613 triexit(1);
3614 }
3615 /* Be careful not to allocate space for element area constraints that */
3616 /* will never be assigned any value (other than the default -1.0). */
3617 if (!b->refine && !b->poly) {
3618 b->vararea = 0;
3619 }
3620 /* Be careful not to add an extra attribute to each element unless the */
3621 /* input supports it (PSLG in, but not refining a preexisting mesh). */
3622 if (b->refine || !b->poly) {
3623 b->regionattrib = 0;
3624 }
3625 /* Regular/weighted triangulations are incompatible with PSLGs */
3626 /* and meshing. */
3627 if (b->weighted && (b->poly || b->quality)) {
3628 b->weighted = 0;
3629 if (!b->quiet) {
3630 printf("Warning: weighted triangulations (-w, -W) are incompatible\n");
3631 printf(" with PSLGs (-p) and meshing (-q, -a, -u). Weights ignored.\n"
3632 );
3633 }
3634 }
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.");
3639 }
3640
3641#ifndef TRILIBRARY
3642 strcpy(b->inpolyfilename, b->innodefilename);
3643 strcpy(b->inelefilename, b->innodefilename);
3644 strcpy(b->areafilename, b->innodefilename);
3645 increment = 0;
3646 strcpy(workstring, b->innodefilename);
3647 j = 1;
3648 while (workstring[j] != '\0') {
3649 if ((workstring[j] == '.') && (workstring[j + 1] != '\0')) {
3650 increment = j + 1;
3651 }
3652 j++;
3653 }
3654 meshnumber = 0;
3655 if (increment > 0) {
3656 j = increment;
3657 do {
3658 if ((workstring[j] >= '0') && (workstring[j] <= '9')) {
3659 meshnumber = meshnumber * 10 + (int) (workstring[j] - '0');
3660 } else {
3661 increment = 0;
3662 }
3663 j++;
3664 } while (workstring[j] != '\0');
3665 }
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");
3698 } else {
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");
3718 }
3719 strcat(b->innodefilename, ".node");
3720 strcat(b->inpolyfilename, ".poly");
3721 strcat(b->inelefilename, ".ele");
3722 strcat(b->areafilename, ".area");
3723#endif /* not TRILIBRARY */
3724}
3725
3728/********* User interaction routines begin here *********/
3729
3730/********* Debugging routines begin here *********/
3734/*****************************************************************************/
3735/* */
3736/* printtriangle() Print out the details of an oriented triangle. */
3737/* */
3738/* I originally wrote this procedure to simplify debugging; it can be */
3739/* called directly from the debugger, and presents information about an */
3740/* oriented triangle in digestible form. It's also used when the */
3741/* highest level of verbosity (`-VVV') is specified. */
3742/* */
3743/*****************************************************************************/
3744
3745#ifdef ANSI_DECLARATORS
3746void printtriangle(struct mesh *m, struct behavior *b, struct otri *t)
3747#else /* not ANSI_DECLARATORS */
3748void printtriangle(m, b, t)
3749struct mesh *m;
3750struct behavior *b;
3751struct otri *t;
3752#endif /* not ANSI_DECLARATORS */
3753
3754{
3755 struct otri printtri;
3756 struct osub printsh;
3757 vertex printvertex;
3758
3759 printf("triangle x%lx with orientation %d:\n", (unsigned long) t->tri,
3760 t->orient);
3761 decode(t->tri[0], printtri);
3762 if (printtri.tri == m->dummytri) {
3763 printf(" [0] = Outer space\n");
3764 } else {
3765 printf(" [0] = x%lx %d\n", (unsigned long) printtri.tri,
3766 printtri.orient);
3767 }
3768 decode(t->tri[1], printtri);
3769 if (printtri.tri == m->dummytri) {
3770 printf(" [1] = Outer space\n");
3771 } else {
3772 printf(" [1] = x%lx %d\n", (unsigned long) printtri.tri,
3773 printtri.orient);
3774 }
3775 decode(t->tri[2], printtri);
3776 if (printtri.tri == m->dummytri) {
3777 printf(" [2] = Outer space\n");
3778 } else {
3779 printf(" [2] = x%lx %d\n", (unsigned long) printtri.tri,
3780 printtri.orient);
3781 }
3782
3783 org(*t, printvertex);
3784 if (printvertex == (vertex) NULL)
3785 printf(" Origin[%d] = NULL\n", (t->orient + 1) % 3 + 3);
3786 else
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);
3793 else
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);
3800 else
3801 printf(" Apex [%d] = x%lx (%.12g, %.12g)\n",
3802 t->orient + 3, (unsigned long) printvertex,
3803 printvertex[0], printvertex[1]);
3804
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,
3809 printsh.ssorient);
3810 }
3811 sdecode(t->tri[7], printsh);
3812 if (printsh.ss != m->dummysub) {
3813 printf(" [7] = x%lx %d\n", (unsigned long) printsh.ss,
3814 printsh.ssorient);
3815 }
3816 sdecode(t->tri[8], printsh);
3817 if (printsh.ss != m->dummysub) {
3818 printf(" [8] = x%lx %d\n", (unsigned long) printsh.ss,
3819 printsh.ssorient);
3820 }
3821 }
3822
3823 if (b->vararea) {
3824 printf(" Area constraint: %.4g\n", areabound(*t));
3825 }
3826}
3827
3828/*****************************************************************************/
3829/* */
3830/* printsubseg() Print out the details of an oriented subsegment. */
3831/* */
3832/* I originally wrote this procedure to simplify debugging; it can be */
3833/* called directly from the debugger, and presents information about an */
3834/* oriented subsegment in digestible form. It's also used when the highest */
3835/* level of verbosity (`-VVV') is specified. */
3836/* */
3837/*****************************************************************************/
3838
3839#ifdef ANSI_DECLARATORS
3840void printsubseg(struct mesh *m, struct behavior *b, struct osub *s)
3841#else /* not ANSI_DECLARATORS */
3842void printsubseg(m, b, s)
3843struct mesh *m;
3844struct behavior *b;
3845struct osub *s;
3846#endif /* not ANSI_DECLARATORS */
3847
3848{
3849 struct osub printsh;
3850 struct otri printtri;
3851 vertex printvertex;
3852
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");
3858 } else {
3859 printf(" [0] = x%lx %d\n", (unsigned long) printsh.ss,
3860 printsh.ssorient);
3861 }
3862 sdecode(s->ss[1], printsh);
3863 if (printsh.ss == m->dummysub) {
3864 printf(" [1] = No subsegment\n");
3865 } else {
3866 printf(" [1] = x%lx %d\n", (unsigned long) printsh.ss,
3867 printsh.ssorient);
3868 }
3869
3870 sorg(*s, printvertex);
3871 if (printvertex == (vertex) NULL)
3872 printf(" Origin[%d] = NULL\n", 2 + s->ssorient);
3873 else
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);
3880 else
3881 printf(" Dest [%d] = x%lx (%.12g, %.12g)\n",
3882 3 - s->ssorient, (unsigned long) printvertex,
3883 printvertex[0], printvertex[1]);
3884
3885 decode(s->ss[6], printtri);
3886 if (printtri.tri == m->dummytri) {
3887 printf(" [6] = Outer space\n");
3888 } else {
3889 printf(" [6] = x%lx %d\n", (unsigned long) printtri.tri,
3890 printtri.orient);
3891 }
3892 decode(s->ss[7], printtri);
3893 if (printtri.tri == m->dummytri) {
3894 printf(" [7] = Outer space\n");
3895 } else {
3896 printf(" [7] = x%lx %d\n", (unsigned long) printtri.tri,
3897 printtri.orient);
3898 }
3899
3900 segorg(*s, printvertex);
3901 if (printvertex == (vertex) NULL)
3902 printf(" Segment origin[%d] = NULL\n", 4 + s->ssorient);
3903 else
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);
3910 else
3911 printf(" Segment dest [%d] = x%lx (%.12g, %.12g)\n",
3912 5 - s->ssorient, (unsigned long) printvertex,
3913 printvertex[0], printvertex[1]);
3914}
3915
3918/********* Debugging routines end here *********/
3919
3920/********* Memory management routines begin here *********/
3924/*****************************************************************************/
3925/* */
3926/* poolzero() Set all of a pool's fields to zero. */
3927/* */
3928/* This procedure should never be called on a pool that has any memory */
3929/* allocated to it, as that memory would leak. */
3930/* */
3931/*****************************************************************************/
3932
3933#ifdef ANSI_DECLARATORS
3934void poolzero(struct memorypool* pool)
3935#else /* not ANSI_DECLARATORS */
3936void poolzero(pool)
3937struct memorypool* pool;
3938#endif /* not ANSI_DECLARATORS */
3939
3940{
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;
3951 pool->items = 0;
3952 pool->maxitems = 0;
3953 pool->unallocateditems = 0;
3954 pool->pathitemsleft = 0;
3955}
3956
3957/*****************************************************************************/
3958/* */
3959/* poolrestart() Deallocate all items in a pool. */
3960/* */
3961/* The pool is returned to its starting state, except that no memory is */
3962/* freed to the operating system. Rather, the previously allocated blocks */
3963/* are ready to be reused. */
3964/* */
3965/*****************************************************************************/
3966
3967#ifdef ANSI_DECLARATORS
3968void poolrestart(struct memorypool* pool)
3969#else /* not ANSI_DECLARATORS */
3970void poolrestart(pool)
3971struct memorypool* pool;
3972#endif /* not ANSI_DECLARATORS */
3973
3974{
3975 unsigned long alignptr;
3976
3977 pool->items = 0;
3978 pool->maxitems = 0;
3979
3980 /* Set the currently active block. */
3981 pool->nowblock = pool->firstblock;
3982 /* Find the first item in the pool. Increment by the size of (VOID *). */
3983 alignptr = (unsigned long) (pool->nowblock + 1);
3984 /* Align the item on an `alignbytes'-byte boundary. */
3985 pool->nextitem = (VOID *)
3986 (alignptr + (unsigned long) pool->alignbytes -
3987 (alignptr % (unsigned long) pool->alignbytes));
3988 /* There are lots of unallocated items left in this block. */
3989 pool->unallocateditems = pool->itemsfirstblock;
3990 /* The stack of deallocated items is empty. */
3991 pool->deaditemstack = (VOID *) NULL;
3992}
3993
3994/*****************************************************************************/
3995/* */
3996/* poolinit() Initialize a pool of memory for allocation of items. */
3997/* */
3998/* This routine initializes the machinery for allocating items. A `pool' */
3999/* is created whose records have size at least `bytecount'. Items will be */
4000/* allocated in `itemcount'-item blocks. Each item is assumed to be a */
4001/* collection of words, and either pointers or floating-point values are */
4002/* assumed to be the "primary" word type. (The "primary" word type is used */
4003/* to determine alignment of items.) If `alignment' isn't zero, all items */
4004/* will be `alignment'-byte aligned in memory. `alignment' must be either */
4005/* a multiple or a factor of the primary word size; powers of two are safe. */
4006/* `alignment' is normally used to create a few unused bits at the bottom */
4007/* of each item's pointer, in which information may be stored. */
4008/* */
4009/* Don't change this routine unless you understand it. */
4010/* */
4011/*****************************************************************************/
4012
4013#ifdef ANSI_DECLARATORS
4014void poolinit(struct memorypool* pool, int bytecount, int itemcount,
4015 int firstitemcount, unsigned alignment)
4016#else /* not ANSI_DECLARATORS */
4017void poolinit(pool, bytecount, itemcount, firstitemcount, alignment)
4018struct memorypool* pool;
4019int bytecount;
4020int itemcount;
4021int firstitemcount;
4022unsigned alignment;
4023#endif /* not ANSI_DECLARATORS */
4024
4025{
4026 /* Find the proper alignment, which must be at least as large as: */
4027 /* - The parameter `alignment'. */
4028 /* - sizeof(VOID *), so the stack of dead items can be maintained */
4029 /* without unaligned accesses. */
4030 if (alignment > sizeof(VOID *)) {
4031 pool->alignbytes = alignment;
4032 } else {
4033 pool->alignbytes = sizeof(VOID *);
4034 }
4035 pool->itembytes = ((bytecount - 1) / pool->alignbytes + 1) *
4036 pool->alignbytes;
4037 pool->itemsperblock = itemcount;
4038 if (firstitemcount == 0) {
4039 pool->itemsfirstblock = itemcount;
4040 } else {
4041 pool->itemsfirstblock = firstitemcount;
4042 }
4043
4044 /* Allocate a block of items. Space for `itemsfirstblock' items and one */
4045 /* pointer (to point to the next block) are allocated, as well as space */
4046 /* to ensure alignment of the items. */
4047 pool->firstblock = (VOID **)
4048 trimalloc(pool->itemsfirstblock * pool->itembytes + (int) sizeof(VOID *) +
4049 pool->alignbytes);
4050 /* Set the next block pointer to NULL. */
4051 *(pool->firstblock) = (VOID *) NULL;
4052 poolrestart(pool);
4053}
4054
4055/*****************************************************************************/
4056/* */
4057/* pooldeinit() Free to the operating system all memory taken by a pool. */
4058/* */
4059/*****************************************************************************/
4060
4061#ifdef ANSI_DECLARATORS
4062void pooldeinit(struct memorypool* pool)
4063#else /* not ANSI_DECLARATORS */
4064void pooldeinit(pool)
4065struct memorypool* pool;
4066#endif /* not ANSI_DECLARATORS */
4067
4068{
4069 while (pool->firstblock != (VOID **) NULL) {
4070 pool->nowblock = (VOID **) *(pool->firstblock);
4071 trifree((VOID *) pool->firstblock);
4072 pool->firstblock = pool->nowblock;
4073 }
4074}
4075
4076/*****************************************************************************/
4077/* */
4078/* poolalloc() Allocate space for an item. */
4079/* */
4080/*****************************************************************************/
4081
4082#ifdef ANSI_DECLARATORS
4083VOID* poolalloc(struct memorypool* pool)
4084#else /* not ANSI_DECLARATORS */
4085VOID* poolalloc(pool)
4086struct memorypool* pool;
4087#endif /* not ANSI_DECLARATORS */
4088
4089{
4090 VOID *newitem;
4091 VOID **newblock;
4092 unsigned long alignptr;
4093
4094 /* First check the linked list of dead items. If the list is not */
4095 /* empty, allocate an item from the list rather than a fresh one. */
4096 if (pool->deaditemstack != (VOID *) NULL) {
4097 newitem = pool->deaditemstack; /* Take first item in list. */
4098 pool->deaditemstack = * (VOID **) pool->deaditemstack;
4099 } else {
4100 /* Check if there are any free items left in the current block. */
4101 if (pool->unallocateditems == 0) {
4102 /* Check if another block must be allocated. */
4103 if (*(pool->nowblock) == (VOID *) NULL) {
4104 /* Allocate a new block of items, pointed to by the previous block. */
4105 newblock = (VOID **) trimalloc(pool->itemsperblock * pool->itembytes +
4106 (int) sizeof(VOID *) +
4107 pool->alignbytes);
4108 *(pool->nowblock) = (VOID *) newblock;
4109 /* The next block pointer is NULL. */
4110 *newblock = (VOID *) NULL;
4111 }
4112
4113 /* Move to the new block. */
4114 pool->nowblock = (VOID **) *(pool->nowblock);
4115 /* Find the first item in the block. */
4116 /* Increment by the size of (VOID *). */
4117 alignptr = (unsigned long) (pool->nowblock + 1);
4118 /* Align the item on an `alignbytes'-byte boundary. */
4119 pool->nextitem = (VOID *)
4120 (alignptr + (unsigned long) pool->alignbytes -
4121 (alignptr % (unsigned long) pool->alignbytes));
4122 /* There are lots of unallocated items left in this block. */
4123 pool->unallocateditems = pool->itemsperblock;
4124 }
4125
4126 /* Allocate a new item. */
4127 newitem = pool->nextitem;
4128 /* Advance `nextitem' pointer to next free item in block. */
4129 pool->nextitem = (VOID *) ((char *) pool->nextitem + pool->itembytes);
4130 pool->unallocateditems--;
4131 pool->maxitems++;
4132 }
4133 pool->items++;
4134 return newitem;
4135}
4136
4137/*****************************************************************************/
4138/* */
4139/* pooldealloc() Deallocate space for an item. */
4140/* */
4141/* The deallocated space is stored in a queue for later reuse. */
4142/* */
4143/*****************************************************************************/
4144
4145#ifdef ANSI_DECLARATORS
4146void pooldealloc(struct memorypool* pool, VOID *dyingitem)
4147#else /* not ANSI_DECLARATORS */
4148void pooldealloc(pool, dyingitem)
4149struct memorypool* pool;
4150VOID *dyingitem;
4151#endif /* not ANSI_DECLARATORS */
4152
4153{
4154 /* Push freshly killed item onto stack. */
4155 *((VOID **) dyingitem) = pool->deaditemstack;
4156 pool->deaditemstack = dyingitem;
4157 pool->items--;
4158}
4159
4160/*****************************************************************************/
4161/* */
4162/* traversalinit() Prepare to traverse the entire list of items. */
4163/* */
4164/* This routine is used in conjunction with traverse(). */
4165/* */
4166/*****************************************************************************/
4167
4168#ifdef ANSI_DECLARATORS
4169void traversalinit(struct memorypool* pool)
4170#else /* not ANSI_DECLARATORS */
4171void traversalinit(pool)
4172struct memorypool* pool;
4173#endif /* not ANSI_DECLARATORS */
4174
4175{
4176 unsigned long alignptr;
4177
4178 /* Begin the traversal in the first block. */
4179 pool->pathblock = pool->firstblock;
4180 /* Find the first item in the block. Increment by the size of (VOID *). */
4181 alignptr = (unsigned long) (pool->pathblock + 1);
4182 /* Align with item on an `alignbytes'-byte boundary. */
4183 pool->pathitem = (VOID *)
4184 (alignptr + (unsigned long) pool->alignbytes -
4185 (alignptr % (unsigned long) pool->alignbytes));
4186 /* Set the number of items left in the current block. */
4187 pool->pathitemsleft = pool->itemsfirstblock;
4188}
4189
4190/*****************************************************************************/
4191/* */
4192/* traverse() Find the next item in the list. */
4193/* */
4194/* This routine is used in conjunction with traversalinit(). Be forewarned */
4195/* that this routine successively returns all items in the list, including */
4196/* deallocated ones on the deaditemqueue. It's up to you to figure out */
4197/* which ones are actually dead. Why? I don't want to allocate extra */
4198/* space just to demarcate dead items. It can usually be done more */
4199/* space-efficiently by a routine that knows something about the structure */
4200/* of the item. */
4201/* */
4202/*****************************************************************************/
4203
4204#ifdef ANSI_DECLARATORS
4205VOID *traverse(struct memorypool* pool)
4206#else /* not ANSI_DECLARATORS */
4207VOID *traverse(pool)
4208struct memorypool* pool;
4209#endif /* not ANSI_DECLARATORS */
4210
4211{
4212 VOID *newitem;
4213 unsigned long alignptr;
4214
4215 /* Stop upon exhausting the list of items. */
4216 if (pool->pathitem == pool->nextitem) {
4217 return (VOID *) NULL;
4218 }
4219
4220 /* Check whether any untraversed items remain in the current block. */
4221 if (pool->pathitemsleft == 0) {
4222 /* Find the next block. */
4223 pool->pathblock = (VOID **) *(pool->pathblock);
4224 /* Find the first item in the block. Increment by the size of (VOID *). */
4225 alignptr = (unsigned long) (pool->pathblock + 1);
4226 /* Align with item on an `alignbytes'-byte boundary. */
4227 pool->pathitem = (VOID *)
4228 (alignptr + (unsigned long) pool->alignbytes -
4229 (alignptr % (unsigned long) pool->alignbytes));
4230 /* Set the number of items left in the current block. */
4231 pool->pathitemsleft = pool->itemsperblock;
4232 }
4233
4234 newitem = pool->pathitem;
4235 /* Find the next item in the block. */
4236 pool->pathitem = (VOID *) ((char *) pool->pathitem + pool->itembytes);
4237 pool->pathitemsleft--;
4238 return newitem;
4239}
4240
4241/*****************************************************************************/
4242/* */
4243/* dummyinit() Initialize the triangle that fills "outer space" and the */
4244/* omnipresent subsegment. */
4245/* */
4246/* The triangle that fills "outer space," called `dummytri', is pointed to */
4247/* by every triangle and subsegment on a boundary (be it outer or inner) of */
4248/* the triangulation. Also, `dummytri' points to one of the triangles on */
4249/* the convex hull (until the holes and concavities are carved), making it */
4250/* possible to find a starting triangle for point location. */
4251/* */
4252/* The omnipresent subsegment, `dummysub', is pointed to by every triangle */
4253/* or subsegment that doesn't have a full complement of real subsegments */
4254/* to point to. */
4255/* */
4256/* `dummytri' and `dummysub' are generally required to fulfill only a few */
4257/* invariants: their vertices must remain NULL and `dummytri' must always */
4258/* be bonded (at offset zero) to some triangle on the convex hull of the */
4259/* mesh, via a boundary edge. Otherwise, the connections of `dummytri' and */
4260/* `dummysub' may change willy-nilly. This makes it possible to avoid */
4261/* writing a good deal of special-case code (in the edge flip, for example) */
4262/* for dealing with the boundary of the mesh, places where no subsegment is */
4263/* present, and so forth. Other entities are frequently bonded to */
4264/* `dummytri' and `dummysub' as if they were real mesh entities, with no */
4265/* harm done. */
4266/* */
4267/*****************************************************************************/
4268
4269#ifdef ANSI_DECLARATORS
4270void dummyinit(struct mesh *m, struct behavior *b, int trianglebytes,
4271 int subsegbytes)
4272#else /* not ANSI_DECLARATORS */
4273void dummyinit(m, b, trianglebytes, subsegbytes)
4274struct mesh *m;
4275struct behavior *b;
4276int trianglebytes;
4277int subsegbytes;
4278#endif /* not ANSI_DECLARATORS */
4279
4280{
4281 unsigned long alignptr;
4282
4283 /* Set up `dummytri', the `triangle' that occupies "outer space." */
4284 m->dummytribase = (triangle *) trimalloc(trianglebytes +
4285 m->triangles.alignbytes);
4286 /* Align `dummytri' on a `triangles.alignbytes'-byte boundary. */
4287 alignptr = (unsigned long) m->dummytribase;
4288 m->dummytri = (triangle *)
4289 (alignptr + (unsigned long) m->triangles.alignbytes -
4290 (alignptr % (unsigned long) m->triangles.alignbytes));
4291 /* Initialize the three adjoining triangles to be "outer space." These */
4292 /* will eventually be changed by various bonding operations, but their */
4293 /* values don't really matter, as long as they can legally be */
4294 /* dereferenced. */
4295 m->dummytri[0] = (triangle) m->dummytri;
4296 m->dummytri[1] = (triangle) m->dummytri;
4297 m->dummytri[2] = (triangle) m->dummytri;
4298 /* Three NULL vertices. */
4299 m->dummytri[3] = (triangle) NULL;
4300 m->dummytri[4] = (triangle) NULL;
4301 m->dummytri[5] = (triangle) NULL;
4302
4303 if (b->usesegments) {
4304 /* Set up `dummysub', the omnipresent subsegment pointed to by any */
4305 /* triangle side or subsegment end that isn't attached to a real */
4306 /* subsegment. */
4307 m->dummysubbase = (subseg *) trimalloc(subsegbytes +
4308 m->subsegs.alignbytes);
4309 /* Align `dummysub' on a `subsegs.alignbytes'-byte boundary. */
4310 alignptr = (unsigned long) m->dummysubbase;
4311 m->dummysub = (subseg *)
4312 (alignptr + (unsigned long) m->subsegs.alignbytes -
4313 (alignptr % (unsigned long) m->subsegs.alignbytes));
4314 /* Initialize the two adjoining subsegments to be the omnipresent */
4315 /* subsegment. These will eventually be changed by various bonding */
4316 /* operations, but their values don't really matter, as long as they */
4317 /* can legally be dereferenced. */
4318 m->dummysub[0] = (subseg) m->dummysub;
4319 m->dummysub[1] = (subseg) m->dummysub;
4320 /* Four NULL vertices. */
4321 m->dummysub[2] = (subseg) NULL;
4322 m->dummysub[3] = (subseg) NULL;
4323 m->dummysub[4] = (subseg) NULL;
4324 m->dummysub[5] = (subseg) NULL;
4325 /* Initialize the two adjoining triangles to be "outer space." */
4326 m->dummysub[6] = (subseg) m->dummytri;
4327 m->dummysub[7] = (subseg) m->dummytri;
4328 /* Set the boundary marker to zero. */
4329 * (int *) (m->dummysub + 8) = 0;
4330
4331 /* Initialize the three adjoining subsegments of `dummytri' to be */
4332 /* the omnipresent subsegment. */
4333 m->dummytri[6] = (triangle) m->dummysub;
4334 m->dummytri[7] = (triangle) m->dummysub;
4335 m->dummytri[8] = (triangle) m->dummysub;
4336 }
4337}
4338
4339/*****************************************************************************/
4340/* */
4341/* initializevertexpool() Calculate the size of the vertex data structure */
4342/* and initialize its memory pool. */
4343/* */
4344/* This routine also computes the `vertexmarkindex' and `vertex2triindex' */
4345/* indices used to find values within each vertex. */
4346/* */
4347/*****************************************************************************/
4348
4349#ifdef ANSI_DECLARATORS
4350void initializevertexpool(struct mesh *m, struct behavior *b)
4351#else /* not ANSI_DECLARATORS */
4352void initializevertexpool(m, b)
4353struct mesh *m;
4354struct behavior *b;
4355#endif /* not ANSI_DECLARATORS */
4356
4357{
4358 int vertexsize;
4359
4360 /* The index within each vertex at which the boundary marker is found, */
4361 /* followed by the vertex type. Ensure the vertex marker is aligned to */
4362 /* a sizeof(int)-byte address. */
4363 m->vertexmarkindex = ((m->mesh_dim + m->nextras) * sizeof(REAL) +
4364 sizeof(int) - 1) /
4365 sizeof(int);
4366 vertexsize = (m->vertexmarkindex + 2) * sizeof(int);
4367 if (b->poly) {
4368 /* The index within each vertex at which a triangle pointer is found. */
4369 /* Ensure the pointer is aligned to a sizeof(triangle)-byte address. */
4370 m->vertex2triindex = (vertexsize + sizeof(triangle) - 1) /
4371 sizeof(triangle);
4372 vertexsize = (m->vertex2triindex + 1) * sizeof(triangle);
4373 }
4374
4375 /* Initialize the pool of vertices. */
4376 poolinit(&m->vertices, vertexsize, VERTEXPERBLOCK,
4377 m->invertices > VERTEXPERBLOCK ? m->invertices : VERTEXPERBLOCK,
4378 sizeof(REAL));
4379}
4380
4381/*****************************************************************************/
4382/* */
4383/* initializetrisubpools() Calculate the sizes of the triangle and */
4384/* subsegment data structures and initialize */
4385/* their memory pools. */
4386/* */
4387/* This routine also computes the `highorderindex', `elemattribindex', and */
4388/* `areaboundindex' indices used to find values within each triangle. */
4389/* */
4390/*****************************************************************************/
4391
4392#ifdef ANSI_DECLARATORS
4393void initializetrisubpools(struct mesh *m, struct behavior *b)
4394#else /* not ANSI_DECLARATORS */
4395void initializetrisubpools(m, b)
4396struct mesh *m;
4397struct behavior *b;
4398#endif /* not ANSI_DECLARATORS */
4399
4400{
4401 unsigned trisize;
4402
4403 /* The index within each triangle at which the extra nodes (above three) */
4404 /* associated with high order elements are found. There are three */
4405 /* pointers to other triangles, three pointers to corners, and possibly */
4406 /* three pointers to subsegments before the extra nodes. */
4407 m->highorderindex = 6 + (b->usesegments * 3);
4408 /* The number of bytes occupied by a triangle. */
4409 trisize = ((b->order + 1) * (b->order + 2) / 2 + (m->highorderindex - 3)) *
4410 sizeof(triangle);
4411 /* The index within each triangle at which its attributes are found, */
4412 /* where the index is measured in REALs. */
4413 m->elemattribindex = (trisize + sizeof(REAL) - 1) / sizeof(REAL);
4414 /* The index within each triangle at which the maximum area constraint */
4415 /* is found, where the index is measured in REALs. Note that if the */
4416 /* `regionattrib' flag is set, an additional attribute will be added. */
4417 m->areaboundindex = m->elemattribindex + m->eextras + b->regionattrib;
4418 /* If triangle attributes or an area bound are needed, increase the number */
4419 /* of bytes occupied by a triangle. */
4420 if (b->vararea) {
4421 trisize = (m->areaboundindex + 1) * sizeof(REAL);
4422 } else if (m->eextras + b->regionattrib > 0) {
4423 trisize = m->areaboundindex * sizeof(REAL);
4424 }
4425 /* If a Voronoi diagram or triangle neighbor graph is requested, make */
4426 /* sure there's room to store an integer index in each triangle. This */
4427 /* integer index can occupy the same space as the subsegment pointers */
4428 /* or attributes or area constraint or extra nodes. */
4429 if ((b->voronoi || b->neighbors) &&
4430 (trisize < 6 * sizeof(triangle) + sizeof(int))) {
4431 trisize = 6 * sizeof(triangle) + sizeof(int);
4432 }
4433
4434 /* Having determined the memory size of a triangle, initialize the pool. */
4435 poolinit(&m->triangles, trisize, TRIPERBLOCK,
4436 (2 * m->invertices - 2) > TRIPERBLOCK ? (2 * m->invertices - 2) :
4437 TRIPERBLOCK, 4);
4438
4439 if (b->usesegments) {
4440 /* Initialize the pool of subsegments. Take into account all eight */
4441 /* pointers and one boundary marker. */
4442 poolinit(&m->subsegs, 8 * sizeof(triangle) + sizeof(int),
4443 SUBSEGPERBLOCK, SUBSEGPERBLOCK, 4);
4444
4445 /* Initialize the "outer space" triangle and omnipresent subsegment. */
4446 dummyinit(m, b, m->triangles.itembytes, m->subsegs.itembytes);
4447 } else {
4448 /* Initialize the "outer space" triangle. */
4449 dummyinit(m, b, m->triangles.itembytes, 0);
4450 }
4451}
4452
4453/*****************************************************************************/
4454/* */
4455/* triangledealloc() Deallocate space for a triangle, marking it dead. */
4456/* */
4457/*****************************************************************************/
4458
4459#ifdef ANSI_DECLARATORS
4460void triangledealloc(struct mesh *m, triangle *dyingtriangle)
4461#else /* not ANSI_DECLARATORS */
4462void triangledealloc(m, dyingtriangle)
4463struct mesh *m;
4464triangle *dyingtriangle;
4465#endif /* not ANSI_DECLARATORS */
4466
4467{
4468 /* Mark the triangle as dead. This makes it possible to detect dead */
4469 /* triangles when traversing the list of all triangles. */
4470 killtri(dyingtriangle);
4471 pooldealloc(&m->triangles, (VOID *) dyingtriangle);
4472}
4473
4474/*****************************************************************************/
4475/* */
4476/* triangletraverse() Traverse the triangles, skipping dead ones. */
4477/* */
4478/*****************************************************************************/
4479
4480#ifdef ANSI_DECLARATORS
4481triangle *triangletraverse(struct mesh *m)
4482#else /* not ANSI_DECLARATORS */
4483triangle *triangletraverse(m)
4484struct mesh *m;
4485#endif /* not ANSI_DECLARATORS */
4486
4487{
4488 triangle *newtriangle;
4489
4490 do {
4491 newtriangle = (triangle *) traverse(&m->triangles);
4492 if (newtriangle == (triangle *) NULL) {
4493 return (triangle *) NULL;
4494 }
4495 } while (deadtri(newtriangle)); /* Skip dead ones. */
4496 return newtriangle;
4497}
4498
4499/*****************************************************************************/
4500/* */
4501/* subsegdealloc() Deallocate space for a subsegment, marking it dead. */
4502/* */
4503/*****************************************************************************/
4504
4505#ifdef ANSI_DECLARATORS
4506void subsegdealloc(struct mesh *m, subseg *dyingsubseg)
4507#else /* not ANSI_DECLARATORS */
4508void subsegdealloc(m, dyingsubseg)
4509struct mesh *m;
4510subseg *dyingsubseg;
4511#endif /* not ANSI_DECLARATORS */
4512
4513{
4514 /* Mark the subsegment as dead. This makes it possible to detect dead */
4515 /* subsegments when traversing the list of all subsegments. */
4516 killsubseg(dyingsubseg);
4517 pooldealloc(&m->subsegs, (VOID *) dyingsubseg);
4518}
4519
4520/*****************************************************************************/
4521/* */
4522/* subsegtraverse() Traverse the subsegments, skipping dead ones. */
4523/* */
4524/*****************************************************************************/
4525
4526#ifdef ANSI_DECLARATORS
4527subseg *subsegtraverse(struct mesh *m)
4528#else /* not ANSI_DECLARATORS */
4529subseg *subsegtraverse(m)
4530struct mesh *m;
4531#endif /* not ANSI_DECLARATORS */
4532
4533{
4534 subseg *newsubseg;
4535
4536 do {
4537 newsubseg = (subseg *) traverse(&m->subsegs);
4538 if (newsubseg == (subseg *) NULL) {
4539 return (subseg *) NULL;
4540 }
4541 } while (deadsubseg(newsubseg)); /* Skip dead ones. */
4542 return newsubseg;
4543}
4544
4545/*****************************************************************************/
4546/* */
4547/* vertexdealloc() Deallocate space for a vertex, marking it dead. */
4548/* */
4549/*****************************************************************************/
4550
4551#ifdef ANSI_DECLARATORS
4552void vertexdealloc(struct mesh *m, vertex dyingvertex)
4553#else /* not ANSI_DECLARATORS */
4554void vertexdealloc(m, dyingvertex)
4555struct mesh *m;
4556vertex dyingvertex;
4557#endif /* not ANSI_DECLARATORS */
4558
4559{
4560 /* Mark the vertex as dead. This makes it possible to detect dead */
4561 /* vertices when traversing the list of all vertices. */
4562 setvertextype(dyingvertex, DEADVERTEX);
4563 pooldealloc(&m->vertices, (VOID *) dyingvertex);
4564}
4565
4566/*****************************************************************************/
4567/* */
4568/* vertextraverse() Traverse the vertices, skipping dead ones. */
4569/* */
4570/*****************************************************************************/
4571
4572#ifdef ANSI_DECLARATORS
4573vertex vertextraverse(struct mesh *m)
4574#else /* not ANSI_DECLARATORS */
4575vertex vertextraverse(m)
4576struct mesh *m;
4577#endif /* not ANSI_DECLARATORS */
4578
4579{
4580 vertex newvertex;
4581
4582 do {
4583 newvertex = (vertex) traverse(&m->vertices);
4584 if (newvertex == (vertex) NULL) {
4585 return (vertex) NULL;
4586 }
4587 } while (vertextype(newvertex) == DEADVERTEX); /* Skip dead ones. */
4588 return newvertex;
4589}
4590
4591/*****************************************************************************/
4592/* */
4593/* badsubsegdealloc() Deallocate space for a bad subsegment, marking it */
4594/* dead. */
4595/* */
4596/*****************************************************************************/
4597
4598#ifndef CDT_ONLY
4599
4600#ifdef ANSI_DECLARATORS
4601void badsubsegdealloc(struct mesh *m, struct badsubseg *dyingseg)
4602#else /* not ANSI_DECLARATORS */
4603void badsubsegdealloc(m, dyingseg)
4604struct mesh *m;
4605struct badsubseg *dyingseg;
4606#endif /* not ANSI_DECLARATORS */
4607
4608{
4609 /* Set subsegment's origin to NULL. This makes it possible to detect dead */
4610 /* badsubsegs when traversing the list of all badsubsegs . */
4611 dyingseg->subsegorg = (vertex) NULL;
4612 pooldealloc(&m->badsubsegs, (VOID *) dyingseg);
4613}
4614
4615#endif /* not CDT_ONLY */
4616
4617/*****************************************************************************/
4618/* */
4619/* badsubsegtraverse() Traverse the bad subsegments, skipping dead ones. */
4620/* */
4621/*****************************************************************************/
4622
4623#ifndef CDT_ONLY
4624
4625#ifdef ANSI_DECLARATORS
4626struct badsubseg *badsubsegtraverse(struct mesh *m)
4627#else /* not ANSI_DECLARATORS */
4628struct badsubseg *badsubsegtraverse(m)
4629struct mesh *m;
4630#endif /* not ANSI_DECLARATORS */
4631
4632{
4633 struct badsubseg *newseg;
4634
4635 do {
4636 newseg = (struct badsubseg *) traverse(&m->badsubsegs);
4637 if (newseg == (struct badsubseg *) NULL) {
4638 return (struct badsubseg *) NULL;
4639 }
4640 } while (newseg->subsegorg == (vertex) NULL); /* Skip dead ones. */
4641 return newseg;
4642}
4643
4644#endif /* not CDT_ONLY */
4645
4646/*****************************************************************************/
4647/* */
4648/* getvertex() Get a specific vertex, by number, from the list. */
4649/* */
4650/* The first vertex is number 'firstnumber'. */
4651/* */
4652/* Note that this takes O(n) time (with a small constant, if VERTEXPERBLOCK */
4653/* is large). I don't care to take the trouble to make it work in constant */
4654/* time. */
4655/* */
4656/*****************************************************************************/
4657
4658#ifdef ANSI_DECLARATORS
4659vertex getvertex(struct mesh *m, struct behavior *b, int number)
4660#else /* not ANSI_DECLARATORS */
4661vertex getvertex(m, b, number)
4662struct mesh *m;
4663struct behavior *b;
4664int number;
4665#endif /* not ANSI_DECLARATORS */
4666
4667{
4668 VOID **getblock;
4669 char *foundvertex;
4670 unsigned long alignptr;
4671 int current;
4672
4673 getblock = m->vertices.firstblock;
4674 current = b->firstnumber;
4675
4676 /* Find the right block. */
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;
4683 }
4684 }
4685
4686 /* Now find the right vertex. */
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));
4691}
4692
4693/*****************************************************************************/
4694/* */
4695/* triangledeinit() Free all remaining allocated memory. */
4696/* */
4697/*****************************************************************************/
4698
4699#ifdef ANSI_DECLARATORS
4700void triangledeinit(struct mesh *m, struct behavior *b)
4701#else /* not ANSI_DECLARATORS */
4702void triangledeinit(m, b)
4703struct mesh *m;
4704struct behavior *b;
4705#endif /* not ANSI_DECLARATORS */
4706
4707{
4708 pooldeinit(&m->triangles);
4709 trifree((VOID *) m->dummytribase);
4710 if (b->usesegments) {
4711 pooldeinit(&m->subsegs);
4712 trifree((VOID *) m->dummysubbase);
4713 }
4714 pooldeinit(&m->vertices);
4715#ifndef CDT_ONLY
4716 if (b->quality) {
4717 pooldeinit(&m->badsubsegs);
4718 if ((b->minangle > 0.0) || b->vararea || b->fixedarea || b->usertest) {
4719 pooldeinit(&m->badtriangles);
4720 pooldeinit(&m->flipstackers);
4721 }
4722 }
4723#endif /* not CDT_ONLY */
4724}
4725
4728/********* Memory management routines end here *********/
4729
4730/********* Constructors begin here *********/
4734/*****************************************************************************/
4735/* */
4736/* maketriangle() Create a new triangle with orientation zero. */
4737/* */
4738/*****************************************************************************/
4739
4740#ifdef ANSI_DECLARATORS
4741void maketriangle(struct mesh *m, struct behavior *b, struct otri *newotri)
4742#else /* not ANSI_DECLARATORS */
4743void maketriangle(m, b, newotri)
4744struct mesh *m;
4745struct behavior *b;
4746struct otri *newotri;
4747#endif /* not ANSI_DECLARATORS */
4748
4749{
4750 int i;
4751
4752 newotri->tri = (triangle *) poolalloc(&m->triangles);
4753 /* Initialize the three adjoining triangles to be "outer space". */
4754 newotri->tri[0] = (triangle) m->dummytri;
4755 newotri->tri[1] = (triangle) m->dummytri;
4756 newotri->tri[2] = (triangle) m->dummytri;
4757 /* Three NULL vertices. */
4758 newotri->tri[3] = (triangle) NULL;
4759 newotri->tri[4] = (triangle) NULL;
4760 newotri->tri[5] = (triangle) NULL;
4761 if (b->usesegments) {
4762 /* Initialize the three adjoining subsegments to be the omnipresent */
4763 /* subsegment. */
4764 newotri->tri[6] = (triangle) m->dummysub;
4765 newotri->tri[7] = (triangle) m->dummysub;
4766 newotri->tri[8] = (triangle) m->dummysub;
4767 }
4768 for (i = 0; i < m->eextras; i++) {
4769 setelemattribute(*newotri, i, 0.0);
4770 }
4771 if (b->vararea) {
4772 setareabound(*newotri, -1.0);
4773 }
4774
4775 newotri->orient = 0;
4776}
4777
4778/*****************************************************************************/
4779/* */
4780/* makesubseg() Create a new subsegment with orientation zero. */
4781/* */
4782/*****************************************************************************/
4783
4784#ifdef ANSI_DECLARATORS
4785void makesubseg(struct mesh *m, struct osub *newsubseg)
4786#else /* not ANSI_DECLARATORS */
4787void makesubseg(m, newsubseg)
4788struct mesh *m;
4789struct osub *newsubseg;
4790#endif /* not ANSI_DECLARATORS */
4791
4792{
4793 newsubseg->ss = (subseg *) poolalloc(&m->subsegs);
4794 /* Initialize the two adjoining subsegments to be the omnipresent */
4795 /* subsegment. */
4796 newsubseg->ss[0] = (subseg) m->dummysub;
4797 newsubseg->ss[1] = (subseg) m->dummysub;
4798 /* Four NULL vertices. */
4799 newsubseg->ss[2] = (subseg) NULL;
4800 newsubseg->ss[3] = (subseg) NULL;
4801 newsubseg->ss[4] = (subseg) NULL;
4802 newsubseg->ss[5] = (subseg) NULL;
4803 /* Initialize the two adjoining triangles to be "outer space." */
4804 newsubseg->ss[6] = (subseg) m->dummytri;
4805 newsubseg->ss[7] = (subseg) m->dummytri;
4806 /* Set the boundary marker to zero. */
4807 setmark(*newsubseg, 0);
4808
4809 newsubseg->ssorient = 0;
4810}
4811
4814/********* Constructors end here *********/
4815
4816/********* Geometric primitives begin here *********/
4820/* The adaptive exact arithmetic geometric predicates implemented herein are */
4821/* described in detail in my paper, "Adaptive Precision Floating-Point */
4822/* Arithmetic and Fast Robust Geometric Predicates." See the header for a */
4823/* full citation. */
4824
4825/* Which of the following two methods of finding the absolute values is */
4826/* fastest is compiler-dependent. A few compilers can inline and optimize */
4827/* the fabs() call; but most will incur the overhead of a function call, */
4828/* which is disastrously slow. A faster way on IEEE machines might be to */
4829/* mask the appropriate bit, but that's difficult to do in C without */
4830/* forcing the value to be stored to memory (rather than be kept in the */
4831/* register to which the optimizer assigned it). */
4832
4833#define Absolute(a) ((a) >= 0.0 ? (a) : -(a))
4834/* #define Absolute(a) fabs(a) */
4835
4836/* Many of the operations are broken up into two pieces, a main part that */
4837/* performs an approximate operation, and a "tail" that computes the */
4838/* roundoff error of that operation. */
4839/* */
4840/* The operations Fast_Two_Sum(), Fast_Two_Diff(), Two_Sum(), Two_Diff(), */
4841/* Split(), and Two_Product() are all implemented as described in the */
4842/* reference. Each of these macros requires certain variables to be */
4843/* defined in the calling routine. The variables `bvirt', `c', `abig', */
4844/* `_i', `_j', `_k', `_l', `_m', and `_n' are declared `INEXACT' because */
4845/* they store the result of an operation that may incur roundoff error. */
4846/* The input parameter `x' (or the highest numbered `x_' parameter) must */
4847/* also be declared `INEXACT'. */
4848
4849#define Fast_Two_Sum_Tail(a, b, x, y) \
4850 bvirt = x - a; \
4851 y = b - bvirt
4852
4853#define Fast_Two_Sum(a, b, x, y) \
4854 x = (REAL) (a + b); \
4855 Fast_Two_Sum_Tail(a, b, x, y)
4856
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; \
4862 y = around + bround
4863
4864#define Two_Sum(a, b, x, y) \
4865 x = (REAL) (a + b); \
4866 Two_Sum_Tail(a, b, x, y)
4867
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; \
4873 y = around + bround
4874
4875#define Two_Diff(a, b, x, y) \
4876 x = (REAL) (a - b); \
4877 Two_Diff_Tail(a, b, x, y)
4878
4879#define Split(a, ahi, alo) \
4880 c = (REAL) (splitter * a); \
4881 abig = (REAL) (c - a); \
4882 ahi = c - abig; \
4883 alo = a - ahi
4884
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
4892
4893#define Two_Product(a, b, x, y) \
4894 x = (REAL) (a * b); \
4895 Two_Product_Tail(a, b, x, y)
4896
4897/* Two_Product_Presplit() is Two_Product() where one of the inputs has */
4898/* already been split. Avoids redundant splitting. */
4899
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
4907
4908/* Square() can be done more quickly than Two_Product(). */
4909
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
4915
4916#define Square(a, x, y) \
4917 x = (REAL) (a * a); \
4918 Square_Tail(a, x, y)
4919
4920/* Macros for summing expansions of various fixed lengths. These are all */
4921/* unrolled versions of Expansion_Sum(). */
4922
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)
4926
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)
4930
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)
4934
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)
4938
4939/* Macro for multiplying a two-component expansion by a single component. */
4940
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)
4947
4948/*****************************************************************************/
4949/* */
4950/* exactinit() Initialize the variables used for exact arithmetic. */
4951/* */
4952/* `epsilon' is the largest power of two such that 1.0 + epsilon = 1.0 in */
4953/* floating-point arithmetic. `epsilon' bounds the relative roundoff */
4954/* error. It is used for floating-point error analysis. */
4955/* */
4956/* `splitter' is used to split floating-point numbers into two half- */
4957/* length significands for exact multiplication. */
4958/* */
4959/* I imagine that a highly optimizing compiler might be too smart for its */
4960/* own good, and somehow cause this routine to fail, if it pretends that */
4961/* floating-point arithmetic is too much like real arithmetic. */
4962/* */
4963/* Don't change this routine unless you fully understand it. */
4964/* */
4965/*****************************************************************************/
4966
4967void exactinit()
4968{
4969 REAL half;
4970 REAL check, lastcheck;
4971 int every_other;
4972#ifdef LINUX
4973 int cword;
4974#endif /* LINUX */
4975
4976#ifdef CPU86
4977#ifdef SINGLE
4978 _control87(_PC_24, _MCW_PC); /* Set FPU control word for single precision. */
4979#else /* not SINGLE */
4980 _control87(_PC_53, _MCW_PC); /* Set FPU control word for double precision. */
4981#endif /* not SINGLE */
4982#endif /* CPU86 */
4983#ifdef LINUX
4984#ifdef SINGLE
4985 /* cword = 4223; */
4986 cword = 4210; /* set FPU control word for single precision */
4987#else /* not SINGLE */
4988 /* cword = 4735; */
4989 cword = 4722; /* set FPU control word for double precision */
4990#endif /* not SINGLE */
4991 _FPU_SETCW(cword);
4992#endif /* LINUX */
4993
4994 every_other = 1;
4995 half = 0.5;
4996 epsilon = 1.0;
4997 splitter = 1.0;
4998 check = 1.0;
4999 /* Repeatedly divide `epsilon' by two until it is too small to add to */
5000 /* one without causing roundoff. (Also check if the sum is equal to */
5001 /* the previous sum, for machines that round up instead of using exact */
5002 /* rounding. Not that these routines will work on such machines.) */
5003 do {
5004 lastcheck = check;
5005 epsilon *= half;
5006 if (every_other) {
5007 splitter *= 2.0;
5008 }
5009 every_other = !every_other;
5010 check = 1.0 + epsilon;
5011 } while ((check != 1.0) && (check != lastcheck));
5012 splitter += 1.0;
5013 /* Error bounds for orientation and incircle tests. */
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;
5024}
5025
5026/*****************************************************************************/
5027/* */
5028/* fast_expansion_sum_zeroelim() Sum two expansions, eliminating zero */
5029/* components from the output expansion. */
5030/* */
5031/* Sets h = e + f. See my Robust Predicates paper for details. */
5032/* */
5033/* If round-to-even is used (as with IEEE 754), maintains the strongly */
5034/* nonoverlapping property. (That is, if e is strongly nonoverlapping, h */
5035/* will be also.) Does NOT maintain the nonoverlapping or nonadjacent */
5036/* properties. */
5037/* */
5038/*****************************************************************************/
5039
5040#ifdef ANSI_DECLARATORS
5041int fast_expansion_sum_zeroelim(int elen, REAL *e, int flen, REAL *f, REAL *h)
5042#else /* not ANSI_DECLARATORS */
5043int fast_expansion_sum_zeroelim(elen, e, flen, f, h) /* h cannot be e or f. */
5044int elen;
5045REAL *e;
5046int flen;
5047REAL *f;
5048REAL *h;
5049#endif /* not ANSI_DECLARATORS */
5050
5051{
5052 REAL Q;
5053 INEXACT REAL Qnew;
5054 INEXACT REAL hh;
5055 INEXACT REAL bvirt;
5056 REAL avirt, bround, around;
5057 int eindex, findex, hindex;
5058 REAL enow, fnow;
5059
5060 enow = e[0];
5061 fnow = f[0];
5062 eindex = findex = 0;
5063 if ((fnow > enow) == (fnow > -enow)) {
5064 Q = enow;
5065 enow = e[++eindex];
5066 } else {
5067 Q = fnow;
5068 fnow = f[++findex];
5069 }
5070 hindex = 0;
5071 if ((eindex < elen) && (findex < flen)) {
5072 if ((fnow > enow) == (fnow > -enow)) {
5073 Fast_Two_Sum(enow, Q, Qnew, hh);
5074 enow = e[++eindex];
5075 } else {
5076 Fast_Two_Sum(fnow, Q, Qnew, hh);
5077 fnow = f[++findex];
5078 }
5079 Q = Qnew;
5080 if (hh != 0.0) {
5081 h[hindex++] = hh;
5082 }
5083 while ((eindex < elen) && (findex < flen)) {
5084 if ((fnow > enow) == (fnow > -enow)) {
5085 Two_Sum(Q, enow, Qnew, hh);
5086 enow = e[++eindex];
5087 } else {
5088 Two_Sum(Q, fnow, Qnew, hh);
5089 fnow = f[++findex];
5090 }
5091 Q = Qnew;
5092 if (hh != 0.0) {
5093 h[hindex++] = hh;
5094 }
5095 }
5096 }
5097 while (eindex < elen) {
5098 Two_Sum(Q, enow, Qnew, hh);
5099 enow = e[++eindex];
5100 Q = Qnew;
5101 if (hh != 0.0) {
5102 h[hindex++] = hh;
5103 }
5104 }
5105 while (findex < flen) {
5106 Two_Sum(Q, fnow, Qnew, hh);
5107 fnow = f[++findex];
5108 Q = Qnew;
5109 if (hh != 0.0) {
5110 h[hindex++] = hh;
5111 }
5112 }
5113 if ((Q != 0.0) || (hindex == 0)) {
5114 h[hindex++] = Q;
5115 }
5116 return hindex;
5117}
5118
5119/*****************************************************************************/
5120/* */
5121/* scale_expansion_zeroelim() Multiply an expansion by a scalar, */
5122/* eliminating zero components from the */
5123/* output expansion. */
5124/* */
5125/* Sets h = be. See my Robust Predicates paper for details. */
5126/* */
5127/* Maintains the nonoverlapping property. If round-to-even is used (as */
5128/* with IEEE 754), maintains the strongly nonoverlapping and nonadjacent */
5129/* properties as well. (That is, if e has one of these properties, so */
5130/* will h.) */
5131/* */
5132/*****************************************************************************/
5133
5134#ifdef ANSI_DECLARATORS
5135int scale_expansion_zeroelim(int elen, REAL *e, REAL b, REAL *h)
5136#else /* not ANSI_DECLARATORS */
5137int scale_expansion_zeroelim(elen, e, b, h) /* e and h cannot be the same. */
5138int elen;
5139REAL *e;
5140REAL b;
5141REAL *h;
5142#endif /* not ANSI_DECLARATORS */
5143
5144{
5145 INEXACT REAL Q, sum;
5146 REAL hh;
5147 INEXACT REAL product1;
5148 REAL product0;
5149 int eindex, hindex;
5150 REAL enow;
5151 INEXACT REAL bvirt;
5152 REAL avirt, bround, around;
5153 INEXACT REAL c;
5154 INEXACT REAL abig;
5155 REAL ahi, alo, bhi, blo;
5156 REAL err1, err2, err3;
5157
5158 Split(b, bhi, blo);
5159 Two_Product_Presplit(e[0], b, bhi, blo, Q, hh);
5160 hindex = 0;
5161 if (hh != 0) {
5162 h[hindex++] = hh;
5163 }
5164 for (eindex = 1; eindex < elen; eindex++) {
5165 enow = e[eindex];
5166 Two_Product_Presplit(enow, b, bhi, blo, product1, product0);
5167 Two_Sum(Q, product0, sum, hh);
5168 if (hh != 0) {
5169 h[hindex++] = hh;
5170 }
5171 Fast_Two_Sum(product1, sum, Q, hh);
5172 if (hh != 0) {
5173 h[hindex++] = hh;
5174 }
5175 }
5176 if ((Q != 0.0) || (hindex == 0)) {
5177 h[hindex++] = Q;
5178 }
5179 return hindex;
5180}
5181
5182/*****************************************************************************/
5183/* */
5184/* estimate() Produce a one-word estimate of an expansion's value. */
5185/* */
5186/* See my Robust Predicates paper for details. */
5187/* */
5188/*****************************************************************************/
5189
5190#ifdef ANSI_DECLARATORS
5191REAL estimate(int elen, REAL *e)
5192#else /* not ANSI_DECLARATORS */
5193REAL estimate(elen, e)
5194int elen;
5195REAL *e;
5196#endif /* not ANSI_DECLARATORS */
5197
5198{
5199 REAL Q;
5200 int eindex;
5201
5202 Q = e[0];
5203 for (eindex = 1; eindex < elen; eindex++) {
5204 Q += e[eindex];
5205 }
5206 return Q;
5207}
5208
5209/*****************************************************************************/
5210/* */
5211/* counterclockwise() Return a positive value if the points pa, pb, and */
5212/* pc occur in counterclockwise order; a negative */
5213/* value if they occur in clockwise order; and zero */
5214/* if they are collinear. The result is also a rough */
5215/* approximation of twice the signed area of the */
5216/* triangle defined by the three points. */
5217/* */
5218/* Uses exact arithmetic if necessary to ensure a correct answer. The */
5219/* result returned is the determinant of a matrix. This determinant is */
5220/* computed adaptively, in the sense that exact arithmetic is used only to */
5221/* the degree it is needed to ensure that the returned value has the */
5222/* correct sign. Hence, this function is usually quite fast, but will run */
5223/* more slowly when the input points are collinear or nearly so. */
5224/* */
5225/* See my Robust Predicates paper for details. */
5226/* */
5227/*****************************************************************************/
5228
5229#ifdef ANSI_DECLARATORS
5230REAL counterclockwiseadapt(vertex pa, vertex pb, vertex pc, REAL detsum)
5231#else /* not ANSI_DECLARATORS */
5232REAL counterclockwiseadapt(pa, pb, pc, detsum)
5233vertex pa;
5234vertex pb;
5235vertex pc;
5236REAL detsum;
5237#endif /* not ANSI_DECLARATORS */
5238
5239{
5240 INEXACT REAL acx, acy, bcx, bcy;
5241 REAL acxtail, acytail, bcxtail, bcytail;
5242 INEXACT REAL detleft, detright;
5243 REAL detlefttail, detrighttail;
5244 REAL det, errbound;
5245 REAL B[4], C1[8], C2[12], D[16];
5246 INEXACT REAL B3;
5247 int C1length, C2length, Dlength;
5248 REAL u[4];
5249 INEXACT REAL u3;
5250 INEXACT REAL s1, t1;
5251 REAL s0, t0;
5252
5253 INEXACT REAL bvirt;
5254 REAL avirt, bround, around;
5255 INEXACT REAL c;
5256 INEXACT REAL abig;
5257 REAL ahi, alo, bhi, blo;
5258 REAL err1, err2, err3;
5259 INEXACT REAL _i, _j;
5260 REAL _0;
5261
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]);
5266
5267 Two_Product(acx, bcy, detleft, detlefttail);
5268 Two_Product(acy, bcx, detright, detrighttail);
5269
5270 Two_Two_Diff(detleft, detlefttail, detright, detrighttail,
5271 B3, B[2], B[1], B[0]);
5272 B[3] = B3;
5273
5274 det = estimate(4, B);
5275 errbound = ccwerrboundB * detsum;
5276 if ((det >= errbound) || (-det >= errbound)) {
5277 return det;
5278 }
5279
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);
5284
5285 if ((acxtail == 0.0) && (acytail == 0.0)
5286 && (bcxtail == 0.0) && (bcytail == 0.0)) {
5287 return det;
5288 }
5289
5290 errbound = ccwerrboundC * detsum + resulterrbound * Absolute(det);
5291 det += (acx * bcytail + bcy * acxtail)
5292 - (acy * bcxtail + bcx * acytail);
5293 if ((det >= errbound) || (-det >= errbound)) {
5294 return det;
5295 }
5296
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]);
5300 u[3] = u3;
5301 C1length = fast_expansion_sum_zeroelim(4, B, 4, u, C1);
5302
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]);
5306 u[3] = u3;
5307 C2length = fast_expansion_sum_zeroelim(C1length, C1, 4, u, C2);
5308
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]);
5312 u[3] = u3;
5313 Dlength = fast_expansion_sum_zeroelim(C2length, C2, 4, u, D);
5314
5315 return(D[Dlength - 1]);
5316}
5317
5318#ifdef ANSI_DECLARATORS
5319REAL counterclockwise(struct mesh *m, struct behavior *b,
5320 vertex pa, vertex pb, vertex pc)
5321#else /* not ANSI_DECLARATORS */
5322REAL counterclockwise(m, b, pa, pb, pc)
5323struct mesh *m;
5324struct behavior *b;
5325vertex pa;
5326vertex pb;
5327vertex pc;
5328#endif /* not ANSI_DECLARATORS */
5329
5330{
5331 REAL detleft, detright, det;
5332 REAL detsum, errbound;
5333
5334 m->counterclockcount++;
5335
5336 detleft = (pa[0] - pc[0]) * (pb[1] - pc[1]);
5337 detright = (pa[1] - pc[1]) * (pb[0] - pc[0]);
5338 det = detleft - detright;
5339
5340 if (b->noexact) {
5341 return det;
5342 }
5343
5344 if (detleft > 0.0) {
5345 if (detright <= 0.0) {
5346 return det;
5347 } else {
5348 detsum = detleft + detright;
5349 }
5350 } else if (detleft < 0.0) {
5351 if (detright >= 0.0) {
5352 return det;
5353 } else {
5354 detsum = -detleft - detright;
5355 }
5356 } else {
5357 return det;
5358 }
5359
5360 errbound = ccwerrboundA * detsum;
5361 if ((det >= errbound) || (-det >= errbound)) {
5362 return det;
5363 }
5364
5365 return counterclockwiseadapt(pa, pb, pc, detsum);
5366}
5367
5368/*****************************************************************************/
5369/* */
5370/* incircle() Return a positive value if the point pd lies inside the */
5371/* circle passing through pa, pb, and pc; a negative value if */
5372/* it lies outside; and zero if the four points are cocircular.*/
5373/* The points pa, pb, and pc must be in counterclockwise */
5374/* order, or the sign of the result will be reversed. */
5375/* */
5376/* Uses exact arithmetic if necessary to ensure a correct answer. The */
5377/* result returned is the determinant of a matrix. This determinant is */
5378/* computed adaptively, in the sense that exact arithmetic is used only to */
5379/* the degree it is needed to ensure that the returned value has the */
5380/* correct sign. Hence, this function is usually quite fast, but will run */
5381/* more slowly when the input points are cocircular or nearly so. */
5382/* */
5383/* See my Robust Predicates paper for details. */
5384/* */
5385/*****************************************************************************/
5386
5387#ifdef ANSI_DECLARATORS
5388REAL incircleadapt(vertex pa, vertex pb, vertex pc, vertex pd, REAL permanent)
5389#else /* not ANSI_DECLARATORS */
5390REAL incircleadapt(pa, pb, pc, pd, permanent)
5391vertex pa;
5392vertex pb;
5393vertex pc;
5394vertex pd;
5395REAL permanent;
5396#endif /* not ANSI_DECLARATORS */
5397
5398{
5399 INEXACT REAL adx, bdx, cdx, ady, bdy, cdy;
5400 REAL det, errbound;
5401
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;
5412 REAL abdet[64];
5413 int ablen;
5414 REAL fin1[1152], fin2[1152];
5415 REAL *finnow, *finother, *finswap;
5416 int finlength;
5417
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;
5424 REAL ti0, tj0;
5425 REAL u[4], v[4];
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;
5449 REAL negate;
5450
5451 INEXACT REAL bvirt;
5452 REAL avirt, bround, around;
5453 INEXACT REAL c;
5454 INEXACT REAL abig;
5455 REAL ahi, alo, bhi, blo;
5456 REAL err1, err2, err3;
5457 INEXACT REAL _i, _j;
5458 REAL _0;
5459
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]);
5466
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]);
5470 bc[3] = bc3;
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);
5476
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]);
5480 ca[3] = ca3;
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);
5486
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]);
5490 ab[3] = ab3;
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);
5496
5497 ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
5498 finlength = fast_expansion_sum_zeroelim(ablen, abdet, clen, cdet, fin1);
5499
5500 det = estimate(finlength, fin1);
5501 errbound = iccerrboundB * permanent;
5502 if ((det >= errbound) || (-det >= errbound)) {
5503 return det;
5504 }
5505
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)) {
5514 return det;
5515 }
5516
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)) {
5528 return det;
5529 }
5530
5531 finnow = fin1;
5532 finother = fin2;
5533
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]);
5539 aa[3] = aa3;
5540 }
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]);
5546 bb[3] = bb3;
5547 }
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]);
5553 cc[3] = cc3;
5554 }
5555
5556 if (adxtail != 0.0) {
5557 axtbclen = scale_expansion_zeroelim(4, bc, adxtail, axtbc);
5558 temp16alen = scale_expansion_zeroelim(axtbclen, axtbc, 2.0 * adx,
5559 temp16a);
5560
5561 axtcclen = scale_expansion_zeroelim(4, cc, adxtail, axtcc);
5562 temp16blen = scale_expansion_zeroelim(axtcclen, axtcc, bdy, temp16b);
5563
5564 axtbblen = scale_expansion_zeroelim(4, bb, adxtail, axtbb);
5565 temp16clen = scale_expansion_zeroelim(axtbblen, axtbb, -cdy, temp16c);
5566
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,
5572 temp48, finother);
5573 finswap = finnow; finnow = finother; finother = finswap;
5574 }
5575 if (adytail != 0.0) {
5576 aytbclen = scale_expansion_zeroelim(4, bc, adytail, aytbc);
5577 temp16alen = scale_expansion_zeroelim(aytbclen, aytbc, 2.0 * ady,
5578 temp16a);
5579
5580 aytbblen = scale_expansion_zeroelim(4, bb, adytail, aytbb);
5581 temp16blen = scale_expansion_zeroelim(aytbblen, aytbb, cdx, temp16b);
5582
5583 aytcclen = scale_expansion_zeroelim(4, cc, adytail, aytcc);
5584 temp16clen = scale_expansion_zeroelim(aytcclen, aytcc, -bdx, temp16c);
5585
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,
5591 temp48, finother);
5592 finswap = finnow; finnow = finother; finother = finswap;
5593 }
5594 if (bdxtail != 0.0) {
5595 bxtcalen = scale_expansion_zeroelim(4, ca, bdxtail, bxtca);
5596 temp16alen = scale_expansion_zeroelim(bxtcalen, bxtca, 2.0 * bdx,
5597 temp16a);
5598
5599 bxtaalen = scale_expansion_zeroelim(4, aa, bdxtail, bxtaa);
5600 temp16blen = scale_expansion_zeroelim(bxtaalen, bxtaa, cdy, temp16b);
5601
5602 bxtcclen = scale_expansion_zeroelim(4, cc, bdxtail, bxtcc);
5603 temp16clen = scale_expansion_zeroelim(bxtcclen, bxtcc, -ady, temp16c);
5604
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,
5610 temp48, finother);
5611 finswap = finnow; finnow = finother; finother = finswap;
5612 }
5613 if (bdytail != 0.0) {
5614 bytcalen = scale_expansion_zeroelim(4, ca, bdytail, bytca);
5615 temp16alen = scale_expansion_zeroelim(bytcalen, bytca, 2.0 * bdy,
5616 temp16a);
5617
5618 bytcclen = scale_expansion_zeroelim(4, cc, bdytail, bytcc);
5619 temp16blen = scale_expansion_zeroelim(bytcclen, bytcc, adx, temp16b);
5620
5621 bytaalen = scale_expansion_zeroelim(4, aa, bdytail, bytaa);
5622 temp16clen = scale_expansion_zeroelim(bytaalen, bytaa, -cdx, temp16c);
5623
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,
5629 temp48, finother);
5630 finswap = finnow; finnow = finother; finother = finswap;
5631 }
5632 if (cdxtail != 0.0) {
5633 cxtablen = scale_expansion_zeroelim(4, ab, cdxtail, cxtab);
5634 temp16alen = scale_expansion_zeroelim(cxtablen, cxtab, 2.0 * cdx,
5635 temp16a);
5636
5637 cxtbblen = scale_expansion_zeroelim(4, bb, cdxtail, cxtbb);
5638 temp16blen = scale_expansion_zeroelim(cxtbblen, cxtbb, ady, temp16b);
5639
5640 cxtaalen = scale_expansion_zeroelim(4, aa, cdxtail, cxtaa);
5641 temp16clen = scale_expansion_zeroelim(cxtaalen, cxtaa, -bdy, temp16c);
5642
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,
5648 temp48, finother);
5649 finswap = finnow; finnow = finother; finother = finswap;
5650 }
5651 if (cdytail != 0.0) {
5652 cytablen = scale_expansion_zeroelim(4, ab, cdytail, cytab);
5653 temp16alen = scale_expansion_zeroelim(cytablen, cytab, 2.0 * cdy,
5654 temp16a);
5655
5656 cytaalen = scale_expansion_zeroelim(4, aa, cdytail, cytaa);
5657 temp16blen = scale_expansion_zeroelim(cytaalen, cytaa, bdx, temp16b);
5658
5659 cytbblen = scale_expansion_zeroelim(4, bb, cdytail, cytbb);
5660 temp16clen = scale_expansion_zeroelim(cytbblen, cytbb, -adx, temp16c);
5661
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,
5667 temp48, finother);
5668 finswap = finnow; finnow = finother; finother = finswap;
5669 }
5670
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]);
5677 u[3] = u3;
5678 negate = -bdy;
5679 Two_Product(cdxtail, negate, ti1, ti0);
5680 negate = -bdytail;
5681 Two_Product(cdx, negate, tj1, tj0);
5682 Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]);
5683 v[3] = v3;
5684 bctlen = fast_expansion_sum_zeroelim(4, u, 4, v, bct);
5685
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]);
5689 bctt[3] = bctt3;
5690 bcttlen = 4;
5691 } else {
5692 bct[0] = 0.0;
5693 bctlen = 1;
5694 bctt[0] = 0.0;
5695 bcttlen = 1;
5696 }
5697
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,
5702 temp32a);
5703 temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
5704 temp32alen, temp32a, temp48);
5705 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
5706 temp48, finother);
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,
5711 temp16a);
5712 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
5713 temp16a, finother);
5714 finswap = finnow; finnow = finother; finother = finswap;
5715 }
5716 if (cdytail != 0.0) {
5717 temp8len = scale_expansion_zeroelim(4, bb, -adxtail, temp8);
5718 temp16alen = scale_expansion_zeroelim(temp8len, temp8, cdytail,
5719 temp16a);
5720 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
5721 temp16a, finother);
5722 finswap = finnow; finnow = finother; finother = finswap;
5723 }
5724
5725 temp32alen = scale_expansion_zeroelim(axtbctlen, axtbct, adxtail,
5726 temp32a);
5727 axtbcttlen = scale_expansion_zeroelim(bcttlen, bctt, adxtail, axtbctt);
5728 temp16alen = scale_expansion_zeroelim(axtbcttlen, axtbctt, 2.0 * adx,
5729 temp16a);
5730 temp16blen = scale_expansion_zeroelim(axtbcttlen, axtbctt, adxtail,
5731 temp16b);
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,
5737 temp64, finother);
5738 finswap = finnow; finnow = finother; finother = finswap;
5739 }
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,
5744 temp32a);
5745 temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
5746 temp32alen, temp32a, temp48);
5747 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
5748 temp48, finother);
5749 finswap = finnow; finnow = finother; finother = finswap;
5750
5751
5752 temp32alen = scale_expansion_zeroelim(aytbctlen, aytbct, adytail,
5753 temp32a);
5754 aytbcttlen = scale_expansion_zeroelim(bcttlen, bctt, adytail, aytbctt);
5755 temp16alen = scale_expansion_zeroelim(aytbcttlen, aytbctt, 2.0 * ady,
5756 temp16a);
5757 temp16blen = scale_expansion_zeroelim(aytbcttlen, aytbctt, adytail,
5758 temp16b);
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,
5764 temp64, finother);
5765 finswap = finnow; finnow = finother; finother = finswap;
5766 }
5767 }
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]);
5774 u[3] = u3;
5775 negate = -cdy;
5776 Two_Product(adxtail, negate, ti1, ti0);
5777 negate = -cdytail;
5778 Two_Product(adx, negate, tj1, tj0);
5779 Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]);
5780 v[3] = v3;
5781 catlen = fast_expansion_sum_zeroelim(4, u, 4, v, cat);
5782
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]);
5786 catt[3] = catt3;
5787 cattlen = 4;
5788 } else {
5789 cat[0] = 0.0;
5790 catlen = 1;
5791 catt[0] = 0.0;
5792 cattlen = 1;
5793 }
5794
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,
5799 temp32a);
5800 temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
5801 temp32alen, temp32a, temp48);
5802 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
5803 temp48, finother);
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,
5808 temp16a);
5809 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
5810 temp16a, finother);
5811 finswap = finnow; finnow = finother; finother = finswap;
5812 }
5813 if (adytail != 0.0) {
5814 temp8len = scale_expansion_zeroelim(4, cc, -bdxtail, temp8);
5815 temp16alen = scale_expansion_zeroelim(temp8len, temp8, adytail,
5816 temp16a);
5817 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
5818 temp16a, finother);
5819 finswap = finnow; finnow = finother; finother = finswap;
5820 }
5821
5822 temp32alen = scale_expansion_zeroelim(bxtcatlen, bxtcat, bdxtail,
5823 temp32a);
5824 bxtcattlen = scale_expansion_zeroelim(cattlen, catt, bdxtail, bxtcatt);
5825 temp16alen = scale_expansion_zeroelim(bxtcattlen, bxtcatt, 2.0 * bdx,
5826 temp16a);
5827 temp16blen = scale_expansion_zeroelim(bxtcattlen, bxtcatt, bdxtail,
5828 temp16b);
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,
5834 temp64, finother);
5835 finswap = finnow; finnow = finother; finother = finswap;
5836 }
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,
5841 temp32a);
5842 temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
5843 temp32alen, temp32a, temp48);
5844 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
5845 temp48, finother);
5846 finswap = finnow; finnow = finother; finother = finswap;
5847
5848
5849 temp32alen = scale_expansion_zeroelim(bytcatlen, bytcat, bdytail,
5850 temp32a);
5851 bytcattlen = scale_expansion_zeroelim(cattlen, catt, bdytail, bytcatt);
5852 temp16alen = scale_expansion_zeroelim(bytcattlen, bytcatt, 2.0 * bdy,
5853 temp16a);
5854 temp16blen = scale_expansion_zeroelim(bytcattlen, bytcatt, bdytail,
5855 temp16b);
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,
5861 temp64, finother);
5862 finswap = finnow; finnow = finother; finother = finswap;
5863 }
5864 }
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]);
5871 u[3] = u3;
5872 negate = -ady;
5873 Two_Product(bdxtail, negate, ti1, ti0);
5874 negate = -adytail;
5875 Two_Product(bdx, negate, tj1, tj0);
5876 Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]);
5877 v[3] = v3;
5878 abtlen = fast_expansion_sum_zeroelim(4, u, 4, v, abt);
5879
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]);
5883 abtt[3] = abtt3;
5884 abttlen = 4;
5885 } else {
5886 abt[0] = 0.0;
5887 abtlen = 1;
5888 abtt[0] = 0.0;
5889 abttlen = 1;
5890 }
5891
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,
5896 temp32a);
5897 temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
5898 temp32alen, temp32a, temp48);
5899 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
5900 temp48, finother);
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,
5905 temp16a);
5906 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
5907 temp16a, finother);
5908 finswap = finnow; finnow = finother; finother = finswap;
5909 }
5910 if (bdytail != 0.0) {
5911 temp8len = scale_expansion_zeroelim(4, aa, -cdxtail, temp8);
5912 temp16alen = scale_expansion_zeroelim(temp8len, temp8, bdytail,
5913 temp16a);
5914 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen,
5915 temp16a, finother);
5916 finswap = finnow; finnow = finother; finother = finswap;
5917 }
5918
5919 temp32alen = scale_expansion_zeroelim(cxtabtlen, cxtabt, cdxtail,
5920 temp32a);
5921 cxtabttlen = scale_expansion_zeroelim(abttlen, abtt, cdxtail, cxtabtt);
5922 temp16alen = scale_expansion_zeroelim(cxtabttlen, cxtabtt, 2.0 * cdx,
5923 temp16a);
5924 temp16blen = scale_expansion_zeroelim(cxtabttlen, cxtabtt, cdxtail,
5925 temp16b);
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,
5931 temp64, finother);
5932 finswap = finnow; finnow = finother; finother = finswap;
5933 }
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,
5938 temp32a);
5939 temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a,
5940 temp32alen, temp32a, temp48);
5941 finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len,
5942 temp48, finother);
5943 finswap = finnow; finnow = finother; finother = finswap;
5944
5945
5946 temp32alen = scale_expansion_zeroelim(cytabtlen, cytabt, cdytail,
5947 temp32a);
5948 cytabttlen = scale_expansion_zeroelim(abttlen, abtt, cdytail, cytabtt);
5949 temp16alen = scale_expansion_zeroelim(cytabttlen, cytabtt, 2.0 * cdy,
5950 temp16a);
5951 temp16blen = scale_expansion_zeroelim(cytabttlen, cytabtt, cdytail,
5952 temp16b);
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,
5958 temp64, finother);
5959 finswap = finnow; finnow = finother; finother = finswap;
5960 }
5961 }
5962
5963 return finnow[finlength - 1];
5964}
5965
5966#ifdef ANSI_DECLARATORS
5967REAL incircle(struct mesh *m, struct behavior *b,
5968 vertex pa, vertex pb, vertex pc, vertex pd)
5969#else /* not ANSI_DECLARATORS */
5970REAL incircle(m, b, pa, pb, pc, pd)
5971struct mesh *m;
5972struct behavior *b;
5973vertex pa;
5974vertex pb;
5975vertex pc;
5976vertex pd;
5977#endif /* not ANSI_DECLARATORS */
5978
5979{
5980 REAL adx, bdx, cdx, ady, bdy, cdy;
5981 REAL bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady;
5982 REAL alift, blift, clift;
5983 REAL det;
5984 REAL permanent, errbound;
5985
5986 m->incirclecount++;
5987
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];
5994
5995 bdxcdy = bdx * cdy;
5996 cdxbdy = cdx * bdy;
5997 alift = adx * adx + ady * ady;
5998
5999 cdxady = cdx * ady;
6000 adxcdy = adx * cdy;
6001 blift = bdx * bdx + bdy * bdy;
6002
6003 adxbdy = adx * bdy;
6004 bdxady = bdx * ady;
6005 clift = cdx * cdx + cdy * cdy;
6006
6007 det = alift * (bdxcdy - cdxbdy)
6008 + blift * (cdxady - adxcdy)
6009 + clift * (adxbdy - bdxady);
6010
6011 if (b->noexact) {
6012 return det;
6013 }
6014
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)) {
6020 return det;
6021 }
6022
6023 return incircleadapt(pa, pb, pc, pd, permanent);
6024}
6025
6026/*****************************************************************************/
6027/* */
6028/* orient3d() Return a positive value if the point pd lies below the */
6029/* plane passing through pa, pb, and pc; "below" is defined so */
6030/* that pa, pb, and pc appear in counterclockwise order when */
6031/* viewed from above the plane. Returns a negative value if */
6032/* pd lies above the plane. Returns zero if the points are */
6033/* coplanar. The result is also a rough approximation of six */
6034/* times the signed volume of the tetrahedron defined by the */
6035/* four points. */
6036/* */
6037/* Uses exact arithmetic if necessary to ensure a correct answer. The */
6038/* result returned is the determinant of a matrix. This determinant is */
6039/* computed adaptively, in the sense that exact arithmetic is used only to */
6040/* the degree it is needed to ensure that the returned value has the */
6041/* correct sign. Hence, this function is usually quite fast, but will run */
6042/* more slowly when the input points are coplanar or nearly so. */
6043/* */
6044/* See my Robust Predicates paper for details. */
6045/* */
6046/*****************************************************************************/
6047
6048#ifdef ANSI_DECLARATORS
6049REAL orient3dadapt(vertex pa, vertex pb, vertex pc, vertex pd,
6050 REAL aheight, REAL bheight, REAL cheight, REAL dheight,
6051 REAL permanent)
6052#else /* not ANSI_DECLARATORS */
6053REAL orient3dadapt(pa, pb, pc, pd,
6054 aheight, bheight, cheight, dheight, permanent)
6055vertex pa;
6056vertex pb;
6057vertex pc;
6058vertex pd;
6059REAL aheight;
6060REAL bheight;
6061REAL cheight;
6062REAL dheight;
6063REAL permanent;
6064#endif /* not ANSI_DECLARATORS */
6065
6066{
6067 INEXACT REAL adx, bdx, cdx, ady, bdy, cdy, adheight, bdheight, cdheight;
6068 REAL det, errbound;
6069
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;
6076 REAL abdet[16];
6077 int ablen;
6078 REAL *finnow, *finother, *finswap;
6079 REAL fin1[192], fin2[192];
6080 int finlength;
6081
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];
6105 INEXACT REAL u3;
6106 int vlength, wlength;
6107 REAL negate;
6108
6109 INEXACT REAL bvirt;
6110 REAL avirt, bround, around;
6111 INEXACT REAL c;
6112 INEXACT REAL abig;
6113 REAL ahi, alo, bhi, blo;
6114 REAL err1, err2, err3;
6115 INEXACT REAL _i, _j, _k;
6116 REAL _0;
6117
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);
6127
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]);
6131 bc[3] = bc3;
6132 alen = scale_expansion_zeroelim(4, bc, adheight, adet);
6133
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]);
6137 ca[3] = ca3;
6138 blen = scale_expansion_zeroelim(4, ca, bdheight, bdet);
6139
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]);
6143 ab[3] = ab3;
6144 clen = scale_expansion_zeroelim(4, ab, cdheight, cdet);
6145
6146 ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
6147 finlength = fast_expansion_sum_zeroelim(ablen, abdet, clen, cdet, fin1);
6148
6149 det = estimate(finlength, fin1);
6150 errbound = o3derrboundB * permanent;
6151 if ((det >= errbound) || (-det >= errbound)) {
6152 return det;
6153 }
6154
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);
6164
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)) {
6170 return det;
6171 }
6172
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)) {
6184 return det;
6185 }
6186
6187 finnow = fin1;
6188 finother = fin2;
6189
6190 if (adxtail == 0.0) {
6191 if (adytail == 0.0) {
6192 at_b[0] = 0.0;
6193 at_blen = 1;
6194 at_c[0] = 0.0;
6195 at_clen = 1;
6196 } else {
6197 negate = -adytail;
6198 Two_Product(negate, bdx, at_blarge, at_b[0]);
6199 at_b[1] = at_blarge;
6200 at_blen = 2;
6201 Two_Product(adytail, cdx, at_clarge, at_c[0]);
6202 at_c[1] = at_clarge;
6203 at_clen = 2;
6204 }
6205 } else {
6206 if (adytail == 0.0) {
6207 Two_Product(adxtail, bdy, at_blarge, at_b[0]);
6208 at_b[1] = at_blarge;
6209 at_blen = 2;
6210 negate = -adxtail;
6211 Two_Product(negate, cdy, at_clarge, at_c[0]);
6212 at_c[1] = at_clarge;
6213 at_clen = 2;
6214 } else {
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;
6220 at_blen = 4;
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;
6226 at_clen = 4;
6227 }
6228 }
6229 if (bdxtail == 0.0) {
6230 if (bdytail == 0.0) {
6231 bt_c[0] = 0.0;
6232 bt_clen = 1;
6233 bt_a[0] = 0.0;
6234 bt_alen = 1;
6235 } else {
6236 negate = -bdytail;
6237 Two_Product(negate, cdx, bt_clarge, bt_c[0]);
6238 bt_c[1] = bt_clarge;
6239 bt_clen = 2;
6240 Two_Product(bdytail, adx, bt_alarge, bt_a[0]);
6241 bt_a[1] = bt_alarge;
6242 bt_alen = 2;
6243 }
6244 } else {
6245 if (bdytail == 0.0) {
6246 Two_Product(bdxtail, cdy, bt_clarge, bt_c[0]);
6247 bt_c[1] = bt_clarge;
6248 bt_clen = 2;
6249 negate = -bdxtail;
6250 Two_Product(negate, ady, bt_alarge, bt_a[0]);
6251 bt_a[1] = bt_alarge;
6252 bt_alen = 2;
6253 } else {
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;
6259 bt_clen = 4;
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;
6265 bt_alen = 4;
6266 }
6267 }
6268 if (cdxtail == 0.0) {
6269 if (cdytail == 0.0) {
6270 ct_a[0] = 0.0;
6271 ct_alen = 1;
6272 ct_b[0] = 0.0;
6273 ct_blen = 1;
6274 } else {
6275 negate = -cdytail;
6276 Two_Product(negate, adx, ct_alarge, ct_a[0]);
6277 ct_a[1] = ct_alarge;
6278 ct_alen = 2;
6279 Two_Product(cdytail, bdx, ct_blarge, ct_b[0]);
6280 ct_b[1] = ct_blarge;
6281 ct_blen = 2;
6282 }
6283 } else {
6284 if (cdytail == 0.0) {
6285 Two_Product(cdxtail, ady, ct_alarge, ct_a[0]);
6286 ct_a[1] = ct_alarge;
6287 ct_alen = 2;
6288 negate = -cdxtail;
6289 Two_Product(negate, bdy, ct_blarge, ct_b[0]);
6290 ct_b[1] = ct_blarge;
6291 ct_blen = 2;
6292 } else {
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;
6298 ct_alen = 4;
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;
6304 ct_blen = 4;
6305 }
6306 }
6307
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,
6311 finother);
6312 finswap = finnow; finnow = finother; finother = finswap;
6313
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,
6317 finother);
6318 finswap = finnow; finnow = finother; finother = finswap;
6319
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,
6323 finother);
6324 finswap = finnow; finnow = finother; finother = finswap;
6325
6326 if (adheighttail != 0.0) {
6327 vlength = scale_expansion_zeroelim(4, bc, adheighttail, v);
6328 finlength = fast_expansion_sum_zeroelim(finlength, finnow, vlength, v,
6329 finother);
6330 finswap = finnow; finnow = finother; finother = finswap;
6331 }
6332 if (bdheighttail != 0.0) {
6333 vlength = scale_expansion_zeroelim(4, ca, bdheighttail, v);
6334 finlength = fast_expansion_sum_zeroelim(finlength, finnow, vlength, v,
6335 finother);
6336 finswap = finnow; finnow = finother; finother = finswap;
6337 }
6338 if (cdheighttail != 0.0) {
6339 vlength = scale_expansion_zeroelim(4, ab, cdheighttail, v);
6340 finlength = fast_expansion_sum_zeroelim(finlength, finnow, vlength, v,
6341 finother);
6342 finswap = finnow; finnow = finother; finother = finswap;
6343 }
6344
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]);
6349 u[3] = u3;
6350 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
6351 finother);
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]);
6356 u[3] = u3;
6357 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
6358 finother);
6359 finswap = finnow; finnow = finother; finother = finswap;
6360 }
6361 }
6362 if (cdytail != 0.0) {
6363 negate = -adxtail;
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]);
6366 u[3] = u3;
6367 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
6368 finother);
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]);
6373 u[3] = u3;
6374 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
6375 finother);
6376 finswap = finnow; finnow = finother; finother = finswap;
6377 }
6378 }
6379 }
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]);
6384 u[3] = u3;
6385 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
6386 finother);
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]);
6391 u[3] = u3;
6392 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
6393 finother);
6394 finswap = finnow; finnow = finother; finother = finswap;
6395 }
6396 }
6397 if (adytail != 0.0) {
6398 negate = -bdxtail;
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]);
6401 u[3] = u3;
6402 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
6403 finother);
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]);
6408 u[3] = u3;
6409 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
6410 finother);
6411 finswap = finnow; finnow = finother; finother = finswap;
6412 }
6413 }
6414 }
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]);
6419 u[3] = u3;
6420 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
6421 finother);
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]);
6426 u[3] = u3;
6427 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
6428 finother);
6429 finswap = finnow; finnow = finother; finother = finswap;
6430 }
6431 }
6432 if (bdytail != 0.0) {
6433 negate = -cdxtail;
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]);
6436 u[3] = u3;
6437 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
6438 finother);
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]);
6443 u[3] = u3;
6444 finlength = fast_expansion_sum_zeroelim(finlength, finnow, 4, u,
6445 finother);
6446 finswap = finnow; finnow = finother; finother = finswap;
6447 }
6448 }
6449 }
6450
6451 if (adheighttail != 0.0) {
6452 wlength = scale_expansion_zeroelim(bctlen, bct, adheighttail, w);
6453 finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
6454 finother);
6455 finswap = finnow; finnow = finother; finother = finswap;
6456 }
6457 if (bdheighttail != 0.0) {
6458 wlength = scale_expansion_zeroelim(catlen, cat, bdheighttail, w);
6459 finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
6460 finother);
6461 finswap = finnow; finnow = finother; finother = finswap;
6462 }
6463 if (cdheighttail != 0.0) {
6464 wlength = scale_expansion_zeroelim(abtlen, abt, cdheighttail, w);
6465 finlength = fast_expansion_sum_zeroelim(finlength, finnow, wlength, w,
6466 finother);
6467 finswap = finnow; finnow = finother; finother = finswap;
6468 }
6469
6470 return finnow[finlength - 1];
6471}
6472
6473#ifdef ANSI_DECLARATORS
6474REAL orient3d(struct mesh *m, struct behavior *b,
6475 vertex pa, vertex pb, vertex pc, vertex pd,
6476 REAL aheight, REAL bheight, REAL cheight, REAL dheight)
6477#else /* not ANSI_DECLARATORS */
6478REAL orient3d(m, b, pa, pb, pc, pd, aheight, bheight, cheight, dheight)
6479struct mesh *m;
6480struct behavior *b;
6481vertex pa;
6482vertex pb;
6483vertex pc;
6484vertex pd;
6485REAL aheight;
6486REAL bheight;
6487REAL cheight;
6488REAL dheight;
6489#endif /* not ANSI_DECLARATORS */
6490
6491{
6492 REAL adx, bdx, cdx, ady, bdy, cdy, adheight, bdheight, cdheight;
6493 REAL bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady;
6494 REAL det;
6495 REAL permanent, errbound;
6496
6497 m->orient3dcount++;
6498
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;
6508
6509 bdxcdy = bdx * cdy;
6510 cdxbdy = cdx * bdy;
6511
6512 cdxady = cdx * ady;
6513 adxcdy = adx * cdy;
6514
6515 adxbdy = adx * bdy;
6516 bdxady = bdx * ady;
6517
6518 det = adheight * (bdxcdy - cdxbdy)
6519 + bdheight * (cdxady - adxcdy)
6520 + cdheight * (adxbdy - bdxady);
6521
6522 if (b->noexact) {
6523 return det;
6524 }
6525
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)) {
6531 return det;
6532 }
6533
6534 return orient3dadapt(pa, pb, pc, pd, aheight, bheight, cheight, dheight,
6535 permanent);
6536}
6537
6538/*****************************************************************************/
6539/* */
6540/* nonregular() Return a positive value if the point pd is incompatible */
6541/* with the circle or plane passing through pa, pb, and pc */
6542/* (meaning that pd is inside the circle or below the */
6543/* plane); a negative value if it is compatible; and zero if */
6544/* the four points are cocircular/coplanar. The points pa, */
6545/* pb, and pc must be in counterclockwise order, or the sign */
6546/* of the result will be reversed. */
6547/* */
6548/* If the -w switch is used, the points are lifted onto the parabolic */
6549/* lifting map, then they are dropped according to their weights, then the */
6550/* 3D orientation test is applied. If the -W switch is used, the points' */
6551/* heights are already provided, so the 3D orientation test is applied */
6552/* directly. If neither switch is used, the incircle test is applied. */
6553/* */
6554/*****************************************************************************/
6555
6556#ifdef ANSI_DECLARATORS
6557REAL nonregular(struct mesh *m, struct behavior *b,
6558 vertex pa, vertex pb, vertex pc, vertex pd)
6559#else /* not ANSI_DECLARATORS */
6560REAL nonregular(m, b, pa, pb, pc, pd)
6561struct mesh *m;
6562struct behavior *b;
6563vertex pa;
6564vertex pb;
6565vertex pc;
6566vertex pd;
6567#endif /* not ANSI_DECLARATORS */
6568
6569{
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]);
6578 } else {
6579 return orient3d(m, b, pa, pb, pc, pd, pa[2], pb[2], pc[2], pd[2]);
6580 }
6581}
6582
6583/*****************************************************************************/
6584/* */
6585/* findcircumcenter() Find the circumcenter of a triangle. */
6586/* */
6587/* The result is returned both in terms of x-y coordinates and xi-eta */
6588/* (barycentric) coordinates. The xi-eta coordinate system is defined in */
6589/* terms of the triangle: the origin of the triangle is the origin of the */
6590/* coordinate system; the destination of the triangle is one unit along the */
6591/* xi axis; and the apex of the triangle is one unit along the eta axis. */
6592/* This procedure also returns the square of the length of the triangle's */
6593/* shortest edge. */
6594/* */
6595/*****************************************************************************/
6596
6597#ifdef ANSI_DECLARATORS
6598void findcircumcenter(struct mesh *m, struct behavior *b,
6599 vertex torg, vertex tdest, vertex tapex,
6600 vertex circumcenter, REAL *xi, REAL *eta, int offcenter)
6601#else /* not ANSI_DECLARATORS */
6602void findcircumcenter(m, b, torg, tdest, tapex, circumcenter, xi, eta,
6603 offcenter)
6604struct mesh *m;
6605struct behavior *b;
6606vertex torg;
6607vertex tdest;
6608vertex tapex;
6609vertex circumcenter;
6610REAL *xi;
6611REAL *eta;
6612int offcenter;
6613#endif /* not ANSI_DECLARATORS */
6614
6615{
6616 REAL xdo, ydo, xao, yao;
6617 REAL dodist, aodist, dadist;
6618 REAL denominator;
6619 REAL dx, dy, dxoff, dyoff;
6620
6621 m->circumcentercount++;
6622
6623 /* Compute the circumcenter of the triangle. */
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]);
6632 if (b->noexact) {
6633 denominator = 0.5 / (xdo * yao - xao * ydo);
6634 } else {
6635 /* Use the counterclockwise() routine to ensure a positive (and */
6636 /* reasonably accurate) result, avoiding any possibility of */
6637 /* division by zero. */
6638 denominator = 0.5 / counterclockwise(m, b, tdest, tapex, torg);
6639 /* Don't count the above as an orientation test. */
6640 m->counterclockcount--;
6641 }
6642 dx = (yao * dodist - ydo * aodist) * denominator;
6643 dy = (xdo * aodist - xao * dodist) * denominator;
6644
6645 /* Find the (squared) length of the triangle's shortest edge. This */
6646 /* serves as a conservative estimate of the insertion radius of the */
6647 /* circumcenter's parent. The estimate is used to ensure that */
6648 /* the algorithm terminates even if very small angles appear in */
6649 /* the input PSLG. */
6650 if ((dodist < aodist) && (dodist < dadist)) {
6651 if (offcenter && (b->offconstant > 0.0)) {
6652 /* Find the position of the off-center, as described by Alper Ungor. */
6653 dxoff = 0.5 * xdo - b->offconstant * ydo;
6654 dyoff = 0.5 * ydo + b->offconstant * xdo;
6655 /* If the off-center is closer to the origin than the */
6656 /* circumcenter, use the off-center instead. */
6657 if (dxoff * dxoff + dyoff * dyoff < dx * dx + dy * dy) {
6658 dx = dxoff;
6659 dy = dyoff;
6660 }
6661 }
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;
6666 /* If the off-center is closer to the origin than the */
6667 /* circumcenter, use the off-center instead. */
6668 if (dxoff * dxoff + dyoff * dyoff < dx * dx + dy * dy) {
6669 dx = dxoff;
6670 dy = dyoff;
6671 }
6672 }
6673 } else {
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]);
6679 /* If the off-center is closer to the destination than the */
6680 /* circumcenter, use the off-center instead. */
6681 if (dxoff * dxoff + dyoff * dyoff <
6682 (dx - xdo) * (dx - xdo) + (dy - ydo) * (dy - ydo)) {
6683 dx = xdo + dxoff;
6684 dy = ydo + dyoff;
6685 }
6686 }
6687 }
6688
6689 circumcenter[0] = torg[0] + dx;
6690 circumcenter[1] = torg[1] + dy;
6691
6692 /* To interpolate vertex attributes for the new vertex inserted at */
6693 /* the circumcenter, define a coordinate system with a xi-axis, */
6694 /* directed from the triangle's origin to its destination, and */
6695 /* an eta-axis, directed from its origin to its apex. */
6696 /* Calculate the xi and eta coordinates of the circumcenter. */
6697 *xi = (yao * dx - xao * dy) * (2.0 * denominator);
6698 *eta = (xdo * dy - ydo * dx) * (2.0 * denominator);
6699}
6700
6703/********* Geometric primitives end here *********/
6704
6705/*****************************************************************************/
6706/* */
6707/* triangleinit() Initialize some variables. */
6708/* */
6709/*****************************************************************************/
6710
6711#ifdef ANSI_DECLARATORS
6712void triangleinit(struct mesh *m)
6713#else /* not ANSI_DECLARATORS */
6714void triangleinit(m)
6715struct mesh *m;
6716#endif /* not ANSI_DECLARATORS */
6717
6718{
6719 poolzero(&m->vertices);
6720 poolzero(&m->triangles);
6721 poolzero(&m->subsegs);
6722 poolzero(&m->viri);
6723 poolzero(&m->badsubsegs);
6724 poolzero(&m->badtriangles);
6725 poolzero(&m->flipstackers);
6726 poolzero(&m->splaynodes);
6727
6728 m->recenttri.tri = (triangle *) NULL; /* No triangle has been visited yet. */
6729 m->undeads = 0; /* No eliminated input vertices yet. */
6730 m->samples = 1; /* Point location should take at least one sample. */
6731 m->checksegments = 0; /* There are no segments in the triangulation yet. */
6732 m->checkquality = 0; /* The quality triangulation stage has not begun. */
6733 m->incirclecount = m->counterclockcount = m->orient3dcount = 0;
6734 m->hyperbolacount = m->circletopcount = m->circumcentercount = 0;
6735 randomseed = 1;
6736
6737 exactinit(); /* Initialize exact arithmetic constants. */
6738}
6739
6740/*****************************************************************************/
6741/* */
6742/* randomnation() Generate a random number between 0 and `choices' - 1. */
6743/* */
6744/* T