- Weird shiftvalues are now hidden to the user, shiftvalues will
[xonotic/netradiant.git] / contrib / gtkgensurf / triangle.c
1 #define ANSI_DECLARATORS
2 /*****************************************************************************/
3 /*                                                                           */
4 /*      888888888        ,o,                          / 888                  */
5 /*         888    88o88o  "    o8888o  88o8888o o88888o 888  o88888o         */
6 /*         888    888    888       88b 888  888 888 888 888 d888  88b        */
7 /*         888    888    888  o88^o888 888  888 "88888" 888 8888oo888        */
8 /*         888    888    888 C888  888 888  888  /      888 q888             */
9 /*         888    888    888  "88o^888 888  888 Cb      888  "88oooo"        */
10 /*                                              "8oo8D                       */
11 /*                                                                           */
12 /*  A Two-Dimensional Quality Mesh Generator and Delaunay Triangulator.      */
13 /*  (triangle.c)                                                             */
14 /*                                                                           */
15 /*  Version 1.3                                                              */
16 /*  July 19, 1996                                                            */
17 /*                                                                           */
18 /*  Copyright 1996                                                           */
19 /*  Jonathan Richard Shewchuk                                                */
20 /*  School of Computer Science                                               */
21 /*  Carnegie Mellon University                                               */
22 /*  5000 Forbes Avenue                                                       */
23 /*  Pittsburgh, Pennsylvania  15213-3891                                     */
24 /*  jrs@cs.cmu.edu                                                           */
25 /*                                                                           */
26 /*  This program may be freely redistributed under the condition that the    */
27 /*    copyright notices (including this entire header and the copyright      */
28 /*    notice printed when the `-h' switch is selected) are not removed, and  */
29 /*    no compensation is received.  Private, research, and institutional     */
30 /*    use is free.  You may distribute modified versions of this code UNDER  */
31 /*    THE CONDITION THAT THIS CODE AND ANY MODIFICATIONS MADE TO IT IN THE   */
32 /*    SAME FILE REMAIN UNDER COPYRIGHT OF THE ORIGINAL AUTHOR, BOTH SOURCE   */
33 /*    AND OBJECT CODE ARE MADE FREELY AVAILABLE WITHOUT CHARGE, AND CLEAR    */
34 /*    NOTICE IS GIVEN OF THE MODIFICATIONS.  Distribution of this code as    */
35 /*    part of a commercial system is permissible ONLY BY DIRECT ARRANGEMENT  */
36 /*    WITH THE AUTHOR.  (If you are not directly supplying this code to a    */
37 /*    customer, and you are instead telling them how they can obtain it for  */
38 /*    free, then you are not required to make any arrangement with me.)      */
39 /*                                                                           */
40 /*  Hypertext instructions for Triangle are available on the Web at          */
41 /*                                                                           */
42 /*      http://www.cs.cmu.edu/~quake/triangle.html                           */
43 /*                                                                           */
44 /*  Some of the references listed below are marked [*].  These are available */
45 /*    for downloading from the Web page                                      */
46 /*                                                                           */
47 /*      http://www.cs.cmu.edu/~quake/triangle.research.html                  */
48 /*                                                                           */
49 /*  A paper discussing some aspects of Triangle is available.  See Jonathan  */
50 /*    Richard Shewchuk, "Triangle:  Engineering a 2D Quality Mesh Generator  */
51 /*    and Delaunay Triangulator," First Workshop on Applied Computational    */
52 /*    Geometry, ACM, May 1996.  [*]                                          */
53 /*                                                                           */
54 /*  Triangle was created as part of the Archimedes project in the School of  */
55 /*    Computer Science at Carnegie Mellon University.  Archimedes is a       */
56 /*    system for compiling parallel finite element solvers.  For further     */
57 /*    information, see Anja Feldmann, Omar Ghattas, John R. Gilbert, Gary L. */
58 /*    Miller, David R. O'Hallaron, Eric J. Schwabe, Jonathan R. Shewchuk,    */
59 /*    and Shang-Hua Teng, "Automated Parallel Solution of Unstructured PDE   */
60 /*    Problems."  To appear in Communications of the ACM, we hope.           */
61 /*                                                                           */
62 /*  The quality mesh generation algorithm is due to Jim Ruppert, "A          */
63 /*    Delaunay Refinement Algorithm for Quality 2-Dimensional Mesh           */
64 /*    Generation," Journal of Algorithms 18(3):548-585, May 1995.  [*]       */
65 /*                                                                           */
66 /*  My implementation of the divide-and-conquer and incremental Delaunay     */
67 /*    triangulation algorithms follows closely the presentation of Guibas    */
68 /*    and Stolfi, even though I use a triangle-based data structure instead  */
69 /*    of their quad-edge data structure.  (In fact, I originally implemented */
70 /*    Triangle using the quad-edge data structure, but switching to a        */
71 /*    triangle-based data structure sped Triangle by a factor of two.)  The  */
72 /*    mesh manipulation primitives and the two aforementioned Delaunay       */
73 /*    triangulation algorithms are described by Leonidas J. Guibas and Jorge */
74 /*    Stolfi, "Primitives for the Manipulation of General Subdivisions and   */
75 /*    the Computation of Voronoi Diagrams," ACM Transactions on Graphics     */
76 /*    4(2):74-123, April 1985.                                               */
77 /*                                                                           */
78 /*  Their O(n log n) divide-and-conquer algorithm is adapted from Der-Tsai   */
79 /*    Lee and Bruce J. Schachter, "Two Algorithms for Constructing the       */
80 /*    Delaunay Triangulation," International Journal of Computer and         */
81 /*    Information Science 9(3):219-242, 1980.  The idea to improve the       */
82 /*    divide-and-conquer algorithm by alternating between vertical and       */
83 /*    horizontal cuts was introduced by Rex A. Dwyer, "A Faster Divide-and-  */
84 /*    Conquer Algorithm for Constructing Delaunay Triangulations,"           */
85 /*    Algorithmica 2(2):137-151, 1987.                                       */
86 /*                                                                           */
87 /*  The incremental insertion algorithm was first proposed by C. L. Lawson,  */
88 /*    "Software for C1 Surface Interpolation," in Mathematical Software III, */
89 /*    John R. Rice, editor, Academic Press, New York, pp. 161-194, 1977.     */
90 /*    For point location, I use the algorithm of Ernst P. Mucke, Isaac       */
91 /*    Saias, and Binhai Zhu, "Fast Randomized Point Location Without         */
92 /*    Preprocessing in Two- and Three-dimensional Delaunay Triangulations,"  */
93 /*    Proceedings of the Twelfth Annual Symposium on Computational Geometry, */
94 /*    ACM, May 1996.  [*]  If I were to randomize the order of point         */
95 /*    insertion (I currently don't bother), their result combined with the   */
96 /*    result of Leonidas J. Guibas, Donald E. Knuth, and Micha Sharir,       */
97 /*    "Randomized Incremental Construction of Delaunay and Voronoi           */
98 /*    Diagrams," Algorithmica 7(4):381-413, 1992, would yield an expected    */
99 /*    O(n^{4/3}) bound on running time.                                      */
100 /*                                                                           */
101 /*  The O(n log n) sweepline Delaunay triangulation algorithm is taken from  */
102 /*    Steven Fortune, "A Sweepline Algorithm for Voronoi Diagrams",          */
103 /*    Algorithmica 2(2):153-174, 1987.  A random sample of edges on the      */
104 /*    boundary of the triangulation are maintained in a splay tree for the   */
105 /*    purpose of point location.  Splay trees are described by Daniel        */
106 /*    Dominic Sleator and Robert Endre Tarjan, "Self-Adjusting Binary Search */
107 /*    Trees," Journal of the ACM 32(3):652-686, July 1985.                   */
108 /*                                                                           */
109 /*  The algorithms for exact computation of the signs of determinants are    */
110 /*    described in Jonathan Richard Shewchuk, "Adaptive Precision Floating-  */
111 /*    Point Arithmetic and Fast Robust Geometric Predicates," Technical      */
112 /*    Report CMU-CS-96-140, School of Computer Science, Carnegie Mellon      */
113 /*    University, Pittsburgh, Pennsylvania, May 1996.  [*]  (Submitted to    */
114 /*    Discrete & Computational Geometry.)  An abbreviated version appears as */
115 /*    Jonathan Richard Shewchuk, "Robust Adaptive Floating-Point Geometric   */
116 /*    Predicates," Proceedings of the Twelfth Annual Symposium on Computa-   */
117 /*    tional Geometry, ACM, May 1996.  [*]  Many of the ideas for my exact   */
118 /*    arithmetic routines originate with Douglas M. Priest, "Algorithms for  */
119 /*    Arbitrary Precision Floating Point Arithmetic," Tenth Symposium on     */
120 /*    Computer Arithmetic, 132-143, IEEE Computer Society Press, 1991.  [*]  */
121 /*    Many of the ideas for the correct evaluation of the signs of           */
122 /*    determinants are taken from Steven Fortune and Christopher J. Van Wyk, */
123 /*    "Efficient Exact Arithmetic for Computational Geometry," Proceedings   */
124 /*    of the Ninth Annual Symposium on Computational Geometry, ACM,          */
125 /*    pp. 163-172, May 1993, and from Steven Fortune, "Numerical Stability   */
126 /*    of Algorithms for 2D Delaunay Triangulations," International Journal   */
127 /*    of Computational Geometry & Applications 5(1-2):193-213, March-June    */
128 /*    1995.                                                                  */
129 /*                                                                           */
130 /*  For definitions of and results involving Delaunay triangulations,        */
131 /*    constrained and conforming versions thereof, and other aspects of      */
132 /*    triangular mesh generation, see the excellent survey by Marshall Bern  */
133 /*    and David Eppstein, "Mesh Generation and Optimal Triangulation," in    */
134 /*    Computing and Euclidean Geometry, Ding-Zhu Du and Frank Hwang,         */
135 /*    editors, World Scientific, Singapore, pp. 23-90, 1992.                 */
136 /*                                                                           */
137 /*  The time for incrementally adding PSLG (planar straight line graph)      */
138 /*    segments to create a constrained Delaunay triangulation is probably    */
139 /*    O(n^2) per segment in the worst case and O(n) per edge in the common   */
140 /*    case, where n is the number of triangles that intersect the segment    */
141 /*    before it is inserted.  This doesn't count point location, which can   */
142 /*    be much more expensive.  (This note does not apply to conforming       */
143 /*    Delaunay triangulations, for which a different method is used to       */
144 /*    insert segments.)                                                      */
145 /*                                                                           */
146 /*  The time for adding segments to a conforming Delaunay triangulation is   */
147 /*    not clear, but does not depend upon n alone.  In some cases, very      */
148 /*    small features (like a point lying next to a segment) can cause a      */
149 /*    single segment to be split an arbitrary number of times.  Of course,   */
150 /*    floating-point precision is a practical barrier to how much this can   */
151 /*    happen.                                                                */
152 /*                                                                           */
153 /*  The time for deleting a point from a Delaunay triangulation is O(n^2) in */
154 /*    the worst case and O(n) in the common case, where n is the degree of   */
155 /*    the point being deleted.  I could improve this to expected O(n) time   */
156 /*    by "inserting" the neighboring vertices in random order, but n is      */
157 /*    usually quite small, so it's not worth the bother.  (The O(n) time     */
158 /*    for random insertion follows from L. Paul Chew, "Building Voronoi      */
159 /*    Diagrams for Convex Polygons in Linear Expected Time," Technical       */
160 /*    Report PCS-TR90-147, Department of Mathematics and Computer Science,   */
161 /*    Dartmouth College, 1990.                                               */
162 /*                                                                           */
163 /*  Ruppert's Delaunay refinement algorithm typically generates triangles    */
164 /*    at a linear rate (constant time per triangle) after the initial        */
165 /*    triangulation is formed.  There may be pathological cases where more   */
166 /*    time is required, but these never arise in practice.                   */
167 /*                                                                           */
168 /*  The segment intersection formulae are straightforward.  If you want to   */
169 /*    see them derived, see Franklin Antonio.  "Faster Line Segment          */
170 /*    Intersection."  In Graphics Gems III (David Kirk, editor), pp. 199-    */
171 /*    202.  Academic Press, Boston, 1992.                                    */
172 /*                                                                           */
173 /*  If you make any improvements to this code, please please please let me   */
174 /*    know, so that I may obtain the improvements.  Even if you don't change */
175 /*    the code, I'd still love to hear what it's being used for.             */
176 /*                                                                           */
177 /*  Disclaimer:  Neither I nor Carnegie Mellon warrant this code in any way  */
178 /*    whatsoever.  This code is provided "as-is".  Use at your own risk.     */
179 /*                                                                           */
180 /*****************************************************************************/
181
182 /* For single precision (which will save some memory and reduce paging),     */
183 /*   define the symbol SINGLE by using the -DSINGLE compiler switch or by    */
184 /*   writing "#define SINGLE" below.                                         */
185 /*                                                                           */
186 /* For double precision (which will allow you to refine meshes to a smaller  */
187 /*   edge length), leave SINGLE undefined.                                   */
188 /*                                                                           */
189 /* Double precision uses more memory, but improves the resolution of the     */
190 /*   meshes you can generate with Triangle.  It also reduces the likelihood  */
191 /*   of a floating exception due to overflow.  Finally, it is much faster    */
192 /*   than single precision on 64-bit architectures like the DEC Alpha.  I    */
193 /*   recommend double precision unless you want to generate a mesh for which */
194 /*   you do not have enough memory.                                          */
195
196 #define SINGLE
197
198 #ifdef SINGLE
199 #define REAL float
200 #else /* not SINGLE */
201 #define REAL double
202 #endif /* not SINGLE */
203
204 /* If yours is not a Unix system, define the NO_TIMER compiler switch to     */
205 /*   remove the Unix-specific timing code.                                   */
206
207 #define NO_TIMER
208
209 /* To insert lots of self-checks for internal errors, define the SELF_CHECK  */
210 /*   symbol.  This will slow down the program significantly.  It is best to  */
211 /*   define the symbol using the -DSELF_CHECK compiler switch, but you could */
212 /*   write "#define SELF_CHECK" below.  If you are modifying this code, I    */
213 /*   recommend you turn self-checks on.                                      */
214
215 /* #define SELF_CHECK */
216
217 /* To compile Triangle as a callable object library (triangle.o), define the */
218 /*   TRILIBRARY symbol.  Read the file triangle.h for details on how to call */
219 /*   the procedure triangulate() that results.                               */
220
221 #define TRILIBRARY
222
223 /* It is possible to generate a smaller version of Triangle using one or     */
224 /*   both of the following symbols.  Define the REDUCED symbol to eliminate  */
225 /*   all features that are primarily of research interest; specifically, the */
226 /*   -i, -F, -s, and -C switches.  Define the CDT_ONLY symbol to eliminate   */
227 /*   all meshing algorithms above and beyond constrained Delaunay            */
228 /*   triangulation; specifically, the -r, -q, -a, -S, and -s switches.       */
229 /*   These reductions are most likely to be useful when generating an object */
230 /*   library (triangle.o) by defining the TRILIBRARY symbol.                 */
231
232 #define REDUCED
233 #define CDT_ONLY
234
235 /* On some machines, the exact arithmetic routines might be defeated by the  */
236 /*   use of internal extended precision floating-point registers.  Sometimes */
237 /*   this problem can be fixed by defining certain values to be volatile,    */
238 /*   thus forcing them to be stored to memory and rounded off.  This isn't   */
239 /*   a great solution, though, as it slows Triangle down.                    */
240 /*                                                                           */
241 /* To try this out, write "#define INEXACT volatile" below.  Normally,       */
242 /*   however, INEXACT should be defined to be nothing.  ("#define INEXACT".) */
243
244 #define INEXACT /* Nothing */
245 /* #define INEXACT volatile */
246
247 /* Maximum number of characters in a file name (including the null).         */
248
249 #define FILENAMESIZE 512
250
251 /* Maximum number of characters in a line read from a file (including the    */
252 /*   null).                                                                  */
253
254 #define INPUTLINESIZE 512
255
256 /* For efficiency, a variety of data structures are allocated in bulk.  The  */
257 /*   following constants determine how many of each structure is allocated   */
258 /*   at once.                                                                */
259
260 #define TRIPERBLOCK 4092           /* Number of triangles allocated at once. */
261 #define SHELLEPERBLOCK 508       /* Number of shell edges allocated at once. */
262 #define POINTPERBLOCK 4092            /* Number of points allocated at once. */
263 #define VIRUSPERBLOCK 1020   /* Number of virus triangles allocated at once. */
264 /* Number of encroached segments allocated at once. */
265 #define BADSEGMENTPERBLOCK 252
266 /* Number of skinny triangles allocated at once. */
267 #define BADTRIPERBLOCK 4092
268 /* Number of splay tree nodes allocated at once. */
269 #define SPLAYNODEPERBLOCK 508
270
271 /* The point marker DEADPOINT is an arbitrary number chosen large enough to  */
272 /*   (hopefully) not conflict with user boundary markers.  Make sure that it */
273 /*   is small enough to fit into your machine's integer size.                */
274
275 #define DEADPOINT -1073741824
276
277 /* The next line is used to outsmart some very stupid compilers.  If your    */
278 /*   compiler is smarter, feel free to replace the "int" with "void".        */
279 /*   Not that it matters.                                                    */
280
281 #define VOID int
282
283 /* Two constants for algorithms based on random sampling.  Both constants    */
284 /*   have been chosen empirically to optimize their respective algorithms.   */
285
286 /* Used for the point location scheme of Mucke, Saias, and Zhu, to decide    */
287 /*   how large a random sample of triangles to inspect.                      */
288 #define SAMPLEFACTOR 11
289 /* Used in Fortune's sweepline Delaunay algorithm to determine what fraction */
290 /*   of boundary edges should be maintained in the splay tree for point      */
291 /*   location on the front.                                                  */
292 #define SAMPLERATE 10
293
294 /* A number that speaks for itself, every kissable digit.                    */
295
296 #define PI 3.141592653589793238462643383279502884197169399375105820974944592308
297
298 /* Another fave.                                                             */
299
300 #define SQUAREROOTTWO 1.4142135623730950488016887242096980785696718753769480732
301
302 /* And here's one for those of you who are intimidated by math.              */
303
304 #define ONETHIRD 0.333333333333333333333333333333333333333333333333333333333333
305
306 #include <stdio.h>
307 #include <string.h>
308 #include <math.h>
309 #ifndef NO_TIMER
310 #include <sys/time.h>
311 #endif /* NO_TIMER */
312 #ifdef TRILIBRARY
313 #include "triangle.h"
314 #endif /* TRILIBRARY */
315
316 /* The following obscenity seems to be necessary to ensure that this program */
317 /* will port to Dec Alphas running OSF/1, because their stdio.h file commits */
318 /* the unpardonable sin of including stdlib.h.  Hence, malloc(), free(), and */
319 /* exit() may or may not already be defined at this point.  I declare these  */
320 /* functions explicitly because some non-ANSI C compilers lack stdlib.h.     */
321
322 #ifndef _STDLIB_H_
323 extern void *malloc();
324 extern void free();
325 extern void exit();
326 extern double strtod();
327 extern long strtol();
328 #endif /* _STDLIB_H_ */
329
330 /* A few forward declarations.                                               */
331
332 void poolrestart();
333 #ifndef TRILIBRARY
334 char *readline();
335 char *findfield();
336 #endif /* not TRILIBRARY */
337
338 /* Labels that signify whether a record consists primarily of pointers or of */
339 /*   floating-point words.  Used to make decisions about data alignment.     */
340
341 enum wordtype {POINTER, FLOATINGPOINT};
342
343 /* Labels that signify the result of point location.  The result of a        */
344 /*   search indicates that the point falls in the interior of a triangle, on */
345 /*   an edge, on a vertex, or outside the mesh.                              */
346
347 enum locateresult {INTRIANGLE, ONEDGE, ONVERTEX, OUTSIDE};
348
349 /* Labels that signify the result of site insertion.  The result indicates   */
350 /*   that the point was inserted with complete success, was inserted but     */
351 /*   encroaches on a segment, was not inserted because it lies on a segment, */
352 /*   or was not inserted because another point occupies the same location.   */
353
354 enum insertsiteresult {SUCCESSFULPOINT, ENCROACHINGPOINT, VIOLATINGPOINT,
355                        DUPLICATEPOINT};
356
357 /* Labels that signify the result of direction finding.  The result          */
358 /*   indicates that a segment connecting the two query points falls within   */
359 /*   the direction triangle, along the left edge of the direction triangle,  */
360 /*   or along the right edge of the direction triangle.                      */
361
362 enum finddirectionresult {WITHIN, LEFTCOLLINEAR, RIGHTCOLLINEAR};
363
364 /* Labels that signify the result of the circumcenter computation routine.   */
365 /*   The return value indicates which edge of the triangle is shortest.      */
366
367 enum circumcenterresult {OPPOSITEORG, OPPOSITEDEST, OPPOSITEAPEX};
368
369 /*****************************************************************************/
370 /*                                                                           */
371 /*  The basic mesh data structures                                           */
372 /*                                                                           */
373 /*  There are three:  points, triangles, and shell edges (abbreviated        */
374 /*  `shelle').  These three data structures, linked by pointers, comprise    */
375 /*  the mesh.  A point simply represents a point in space and its properties.*/
376 /*  A triangle is a triangle.  A shell edge is a special data structure used */
377 /*  to represent impenetrable segments in the mesh (including the outer      */
378 /*  boundary, boundaries of holes, and internal boundaries separating two    */
379 /*  triangulated regions).  Shell edges represent boundaries defined by the  */
380 /*  user that triangles may not lie across.                                  */
381 /*                                                                           */
382 /*  A triangle consists of a list of three vertices, a list of three         */
383 /*  adjoining triangles, a list of three adjoining shell edges (when shell   */
384 /*  edges are used), an arbitrary number of optional user-defined floating-  */
385 /*  point attributes, and an optional area constraint.  The latter is an     */
386 /*  upper bound on the permissible area of each triangle in a region, used   */
387 /*  for mesh refinement.                                                     */
388 /*                                                                           */
389 /*  For a triangle on a boundary of the mesh, some or all of the neighboring */
390 /*  triangles may not be present.  For a triangle in the interior of the     */
391 /*  mesh, often no neighboring shell edges are present.  Such absent         */
392 /*  triangles and shell edges are never represented by NULL pointers; they   */
393 /*  are represented by two special records:  `dummytri', the triangle that   */
394 /*  fills "outer space", and `dummysh', the omnipresent shell edge.          */
395 /*  `dummytri' and `dummysh' are used for several reasons; for instance,     */
396 /*  they can be dereferenced and their contents examined without causing the */
397 /*  memory protection exception that would occur if NULL were dereferenced.  */
398 /*                                                                           */
399 /*  However, it is important to understand that a triangle includes other    */
400 /*  information as well.  The pointers to adjoining vertices, triangles, and */
401 /*  shell edges are ordered in a way that indicates their geometric relation */
402 /*  to each other.  Furthermore, each of these pointers contains orientation */
403 /*  information.  Each pointer to an adjoining triangle indicates which face */
404 /*  of that triangle is contacted.  Similarly, each pointer to an adjoining  */
405 /*  shell edge indicates which side of that shell edge is contacted, and how */
406 /*  the shell edge is oriented relative to the triangle.                     */
407 /*                                                                           */
408 /*  Shell edges are found abutting edges of triangles; either sandwiched     */
409 /*  between two triangles, or resting against one triangle on an exterior    */
410 /*  boundary or hole boundary.                                               */
411 /*                                                                           */
412 /*  A shell edge consists of a list of two vertices, a list of two           */
413 /*  adjoining shell edges, and a list of two adjoining triangles.  One of    */
414 /*  the two adjoining triangles may not be present (though there should      */
415 /*  always be one), and neighboring shell edges might not be present.        */
416 /*  Shell edges also store a user-defined integer "boundary marker".         */
417 /*  Typically, this integer is used to indicate what sort of boundary        */
418 /*  conditions are to be applied at that location in a finite element        */
419 /*  simulation.                                                              */
420 /*                                                                           */
421 /*  Like triangles, shell edges maintain information about the relative      */
422 /*  orientation of neighboring objects.                                      */
423 /*                                                                           */
424 /*  Points are relatively simple.  A point is a list of floating point       */
425 /*  numbers, starting with the x, and y coordinates, followed by an          */
426 /*  arbitrary number of optional user-defined floating-point attributes,     */
427 /*  followed by an integer boundary marker.  During the segment insertion    */
428 /*  phase, there is also a pointer from each point to a triangle that may    */
429 /*  contain it.  Each pointer is not always correct, but when one is, it     */
430 /*  speeds up segment insertion.  These pointers are assigned values once    */
431 /*  at the beginning of the segment insertion phase, and are not used or     */
432 /*  updated at any other time.  Edge swapping during segment insertion will  */
433 /*  render some of them incorrect.  Hence, don't rely upon them for          */
434 /*  anything.  For the most part, points do not have any information about   */
435 /*  what triangles or shell edges they are linked to.                        */
436 /*                                                                           */
437 /*****************************************************************************/
438
439 /*****************************************************************************/
440 /*                                                                           */
441 /*  Handles                                                                  */
442 /*                                                                           */
443 /*  The oriented triangle (`triedge') and oriented shell edge (`edge') data  */
444 /*  structures defined below do not themselves store any part of the mesh.   */
445 /*  The mesh itself is made of `triangle's, `shelle's, and `point's.         */
446 /*                                                                           */
447 /*  Oriented triangles and oriented shell edges will usually be referred to  */
448 /*  as "handles".  A handle is essentially a pointer into the mesh; it       */
449 /*  allows you to "hold" one particular part of the mesh.  Handles are used  */
450 /*  to specify the regions in which one is traversing and modifying the mesh.*/
451 /*  A single `triangle' may be held by many handles, or none at all.  (The   */
452 /*  latter case is not a memory leak, because the triangle is still          */
453 /*  connected to other triangles in the mesh.)                               */
454 /*                                                                           */
455 /*  A `triedge' is a handle that holds a triangle.  It holds a specific side */
456 /*  of the triangle.  An `edge' is a handle that holds a shell edge.  It     */
457 /*  holds either the left or right side of the edge.                         */
458 /*                                                                           */
459 /*  Navigation about the mesh is accomplished through a set of mesh          */
460 /*  manipulation primitives, further below.  Many of these primitives take   */
461 /*  a handle and produce a new handle that holds the mesh near the first     */
462 /*  handle.  Other primitives take two handles and glue the corresponding    */
463 /*  parts of the mesh together.  The exact position of the handles is        */
464 /*  important.  For instance, when two triangles are glued together by the   */
465 /*  bond() primitive, they are glued by the sides on which the handles lie.  */
466 /*                                                                           */
467 /*  Because points have no information about which triangles they are        */
468 /*  attached to, I commonly represent a point by use of a handle whose       */
469 /*  origin is the point.  A single handle can simultaneously represent a     */
470 /*  triangle, an edge, and a point.                                          */
471 /*                                                                           */
472 /*****************************************************************************/
473
474 /* The triangle data structure.  Each triangle contains three pointers to    */
475 /*   adjoining triangles, plus three pointers to vertex points, plus three   */
476 /*   pointers to shell edges (defined below; these pointers are usually      */
477 /*   `dummysh').  It may or may not also contain user-defined attributes     */
478 /*   and/or a floating-point "area constraint".  It may also contain extra   */
479 /*   pointers for nodes, when the user asks for high-order elements.         */
480 /*   Because the size and structure of a `triangle' is not decided until     */
481 /*   runtime, I haven't simply defined the type `triangle' to be a struct.   */
482
483 typedef REAL **triangle;            /* Really:  typedef triangle *triangle   */
484
485 /* An oriented triangle:  includes a pointer to a triangle and orientation.  */
486 /*   The orientation denotes an edge of the triangle.  Hence, there are      */
487 /*   three possible orientations.  By convention, each edge is always        */
488 /*   directed to point counterclockwise about the corresponding triangle.    */
489
490 struct triedge {
491   triangle *tri;
492   int orient;                                         /* Ranges from 0 to 2. */
493 };
494
495 /* The shell data structure.  Each shell edge contains two pointers to       */
496 /*   adjoining shell edges, plus two pointers to vertex points, plus two     */
497 /*   pointers to adjoining triangles, plus one shell marker.                 */
498
499 typedef REAL **shelle;                  /* Really:  typedef shelle *shelle   */
500
501 /* An oriented shell edge:  includes a pointer to a shell edge and an        */
502 /*   orientation.  The orientation denotes a side of the edge.  Hence, there */
503 /*   are two possible orientations.  By convention, the edge is always       */
504 /*   directed so that the "side" denoted is the right side of the edge.      */
505
506 struct edge {
507   shelle *sh;
508   int shorient;                                       /* Ranges from 0 to 1. */
509 };
510
511 /* The point data structure.  Each point is actually an array of REALs.      */
512 /*   The number of REALs is unknown until runtime.  An integer boundary      */
513 /*   marker, and sometimes a pointer to a triangle, is appended after the    */
514 /*   REALs.                                                                  */
515
516 typedef REAL *point;
517
518 /* A queue used to store encroached segments.  Each segment's vertices are   */
519 /*   stored so that one can check whether a segment is still the same.       */
520
521 struct badsegment {
522   struct edge encsegment;                          /* An encroached segment. */
523   point segorg, segdest;                                /* The two vertices. */
524   struct badsegment *nextsegment;     /* Pointer to next encroached segment. */
525 };
526
527 /* A queue used to store bad triangles.  The key is the square of the cosine */
528 /*   of the smallest angle of the triangle.  Each triangle's vertices are    */
529 /*   stored so that one can check whether a triangle is still the same.      */
530
531 struct badface {
532   struct triedge badfacetri;                              /* A bad triangle. */
533   REAL key;                             /* cos^2 of smallest (apical) angle. */
534   point faceorg, facedest, faceapex;                  /* The three vertices. */
535   struct badface *nextface;                 /* Pointer to next bad triangle. */
536 };
537
538 /* A node in a heap used to store events for the sweepline Delaunay          */
539 /*   algorithm.  Nodes do not point directly to their parents or children in */
540 /*   the heap.  Instead, each node knows its position in the heap, and can   */
541 /*   look up its parent and children in a separate array.  The `eventptr'    */
542 /*   points either to a `point' or to a triangle (in encoded format, so that */
543 /*   an orientation is included).  In the latter case, the origin of the     */
544 /*   oriented triangle is the apex of a "circle event" of the sweepline      */
545 /*   algorithm.  To distinguish site events from circle events, all circle   */
546 /*   events are given an invalid (smaller than `xmin') x-coordinate `xkey'.  */
547
548 struct event {
549   REAL xkey, ykey;                              /* Coordinates of the event. */
550   VOID *eventptr;       /* Can be a point or the location of a circle event. */
551   int heapposition;              /* Marks this event's position in the heap. */
552 };
553
554 /* A node in the splay tree.  Each node holds an oriented ghost triangle     */
555 /*   that represents a boundary edge of the growing triangulation.  When a   */
556 /*   circle event covers two boundary edges with a triangle, so that they    */
557 /*   are no longer boundary edges, those edges are not immediately deleted   */
558 /*   from the tree; rather, they are lazily deleted when they are next       */
559 /*   encountered.  (Since only a random sample of boundary edges are kept    */
560 /*   in the tree, lazy deletion is faster.)  `keydest' is used to verify     */
561 /*   that a triangle is still the same as when it entered the splay tree; if */
562 /*   it has been rotated (due to a circle event), it no longer represents a  */
563 /*   boundary edge and should be deleted.                                    */
564
565 struct splaynode {
566   struct triedge keyedge;                  /* Lprev of an edge on the front. */
567   point keydest;            /* Used to verify that splay node is still live. */
568   struct splaynode *lchild, *rchild;              /* Children in splay tree. */
569 };
570
571 /* A type used to allocate memory.  firstblock is the first block of items.  */
572 /*   nowblock is the block from which items are currently being allocated.   */
573 /*   nextitem points to the next slab of free memory for an item.            */
574 /*   deaditemstack is the head of a linked list (stack) of deallocated items */
575 /*   that can be recycled.  unallocateditems is the number of items that     */
576 /*   remain to be allocated from nowblock.                                   */
577 /*                                                                           */
578 /* Traversal is the process of walking through the entire list of items, and */
579 /*   is separate from allocation.  Note that a traversal will visit items on */
580 /*   the "deaditemstack" stack as well as live items.  pathblock points to   */
581 /*   the block currently being traversed.  pathitem points to the next item  */
582 /*   to be traversed.  pathitemsleft is the number of items that remain to   */
583 /*   be traversed in pathblock.                                              */
584 /*                                                                           */
585 /* itemwordtype is set to POINTER or FLOATINGPOINT, and is used to suggest   */
586 /*   what sort of word the record is primarily made up of.  alignbytes       */
587 /*   determines how new records should be aligned in memory.  itembytes and  */
588 /*   itemwords are the length of a record in bytes (after rounding up) and   */
589 /*   words.  itemsperblock is the number of items allocated at once in a     */
590 /*   single block.  items is the number of currently allocated items.        */
591 /*   maxitems is the maximum number of items that have been allocated at     */
592 /*   once; it is the current number of items plus the number of records kept */
593 /*   on deaditemstack.                                                       */
594
595 struct memorypool {
596   VOID **firstblock, **nowblock;
597   VOID *nextitem;
598   VOID *deaditemstack;
599   VOID **pathblock;
600   VOID *pathitem;
601   enum wordtype itemwordtype;
602   int alignbytes;
603   int itembytes, itemwords;
604   int itemsperblock;
605   long items, maxitems;
606   int unallocateditems;
607   int pathitemsleft;
608 };
609
610 /* Variables used to allocate memory for triangles, shell edges, points,     */
611 /*   viri (triangles being eaten), bad (encroached) segments, bad (skinny    */
612 /*   or too large) triangles, and splay tree nodes.                          */
613
614 static struct memorypool triangles;
615 static struct memorypool shelles;
616 static struct memorypool points;
617 static struct memorypool viri;
618 static struct memorypool badsegments;
619 static struct memorypool badtriangles;
620 static struct memorypool splaynodes;
621
622 /* Variables that maintain the bad triangle queues.  The tails are pointers  */
623 /*   to the pointers that have to be filled in to enqueue an item.           */
624
625 static struct badface *queuefront[64];
626 static struct badface **queuetail[64];
627
628 static REAL xmin, xmax, ymin, ymax;                              /* x and y bounds. */
629 static REAL xminextreme;        /* Nonexistent x value used as a flag in sweepline. */
630 static int inpoints;                                     /* Number of input points. */
631 static int inelements;                                /* Number of input triangles. */
632 static int insegments;                                 /* Number of input segments. */
633 static int holes;                                         /* Number of input holes. */
634 static int regions;                                     /* Number of input regions. */
635 static long edges;                                       /* Number of output edges. */
636 static int mesh_dim;                                  /* Dimension (ought to be 2). */
637 static int nextras;                              /* Number of attributes per point. */
638 static int eextras;                           /* Number of attributes per triangle. */
639 static long hullsize;                            /* Number of edges of convex hull. */
640 static int triwords;                                   /* Total words per triangle. */
641 static int shwords;                                  /* Total words per shell edge. */
642 static int pointmarkindex;             /* Index to find boundary marker of a point. */
643 static int point2triindex;         /* Index to find a triangle adjacent to a point. */
644 static int highorderindex;    /* Index to find extra nodes for high-order elements. */
645 static int elemattribindex;              /* Index to find attributes of a triangle. */
646 static int areaboundindex;               /* Index to find area bound of a triangle. */
647 static int checksegments;           /* Are there segments in the triangulation yet? */
648 static int readnodefile;                             /* Has a .node file been read? */
649 static long samples;                /* Number of random samples for point location. */
650 static unsigned long randomseed;                     /* Current random number seed. */
651
652 static REAL splitter;       /* Used to split REAL factors for exact multiplication. */
653 static REAL epsilon;                             /* Floating-point machine epsilon. */
654 static REAL resulterrbound;
655 static REAL ccwerrboundA, ccwerrboundB, ccwerrboundC;
656 static REAL iccerrboundA, iccerrboundB, iccerrboundC;
657
658 static long incirclecount;                   /* Number of incircle tests performed. */
659 static long counterclockcount;       /* Number of counterclockwise tests performed. */
660 static long hyperbolacount;        /* Number of right-of-hyperbola tests performed. */
661 static long circumcentercount;    /* Number of circumcenter calculations performed. */
662 static long circletopcount;         /* Number of circle top calculations performed. */
663
664 /* Switches for the triangulator.                                            */
665 /*   poly: -p switch.  refine: -r switch.                                    */
666 /*   quality: -q switch.                                                     */
667 /*     minangle: minimum angle bound, specified after -q switch.             */
668 /*     goodangle: cosine squared of minangle.                                */
669 /*   vararea: -a switch without number.                                      */
670 /*   fixedarea: -a switch with number.                                       */
671 /*     maxarea: maximum area bound, specified after -a switch.               */
672 /*   regionattrib: -A switch.  convex: -c switch.                            */
673 /*   firstnumber: inverse of -z switch.  All items are numbered starting     */
674 /*     from firstnumber.                                                     */
675 /*   edgesout: -e switch.  voronoi: -v switch.                               */
676 /*   neighbors: -n switch.  geomview: -g switch.                             */
677 /*   nobound: -B switch.  nopolywritten: -P switch.                          */
678 /*   nonodewritten: -N switch.  noelewritten: -E switch.                     */
679 /*   noiterationnum: -I switch.  noholes: -O switch.                         */
680 /*   noexact: -X switch.                                                     */
681 /*   order: element order, specified after -o switch.                        */
682 /*   nobisect: count of how often -Y switch is selected.                     */
683 /*   steiner: maximum number of Steiner points, specified after -S switch.   */
684 /*     steinerleft: number of Steiner points not yet used.                   */
685 /*   incremental: -i switch.  sweepline: -F switch.                          */
686 /*   dwyer: inverse of -l switch.                                            */
687 /*   splitseg: -s switch.                                                    */
688 /*   docheck: -C switch.                                                     */
689 /*   quiet: -Q switch.  verbose: count of how often -V switch is selected.   */
690 /*   useshelles: -p, -r, -q, or -c switch; determines whether shell edges    */
691 /*     are used at all.                                                      */
692 /*                                                                           */
693 /* Read the instructions to find out the meaning of these switches.          */
694
695 static int poly, refine, quality, vararea, fixedarea, regionattrib, convex;
696 static int firstnumber;
697 static int edgesout, voronoi, neighbors, geomview;
698 static int nobound, nopolywritten, nonodewritten, noelewritten, noiterationnum;
699 static int noholes, noexact;
700 static int incremental, sweepline, dwyer;
701 static int splitseg;
702 static int docheck;
703 static int quiet, verbose;
704 static int useshelles;
705 static int order;
706 static int nobisect;
707 static int steiner, steinerleft;
708 static REAL minangle, goodangle;
709 static REAL maxarea;
710
711 /* Variables for file names.                                                 */
712
713 #ifndef TRILIBRARY
714 char innodefilename[FILENAMESIZE];
715 char inelefilename[FILENAMESIZE];
716 char inpolyfilename[FILENAMESIZE];
717 char areafilename[FILENAMESIZE];
718 char outnodefilename[FILENAMESIZE];
719 char outelefilename[FILENAMESIZE];
720 char outpolyfilename[FILENAMESIZE];
721 char edgefilename[FILENAMESIZE];
722 char vnodefilename[FILENAMESIZE];
723 char vedgefilename[FILENAMESIZE];
724 char neighborfilename[FILENAMESIZE];
725 char offfilename[FILENAMESIZE];
726 #endif /* not TRILIBRARY */
727
728 /* Triangular bounding box points.                                           */
729
730 static point infpoint1, infpoint2, infpoint3;
731
732 /* Pointer to the `triangle' that occupies all of "outer space".             */
733
734 static triangle *dummytri;
735 static triangle *dummytribase;      /* Keep base address so we can free() it later. */
736
737 /* Pointer to the omnipresent shell edge.  Referenced by any triangle or     */
738 /*   shell edge that isn't really connected to a shell edge at that          */
739 /*   location.                                                               */
740
741 static shelle *dummysh;
742 static shelle *dummyshbase;         /* Keep base address so we can free() it later. */
743
744 /* Pointer to a recently visited triangle.  Improves point location if       */
745 /*   proximate points are inserted sequentially.                             */
746
747 static struct triedge recenttri;
748
749 /*****************************************************************************/
750 /*                                                                           */
751 /*  Mesh manipulation primitives.  Each triangle contains three pointers to  */
752 /*  other triangles, with orientations.  Each pointer points not to the      */
753 /*  first byte of a triangle, but to one of the first three bytes of a       */
754 /*  triangle.  It is necessary to extract both the triangle itself and the   */
755 /*  orientation.  To save memory, I keep both pieces of information in one   */
756 /*  pointer.  To make this possible, I assume that all triangles are aligned */
757 /*  to four-byte boundaries.  The `decode' routine below decodes a pointer,  */
758 /*  extracting an orientation (in the range 0 to 2) and a pointer to the     */
759 /*  beginning of a triangle.  The `encode' routine compresses a pointer to a */
760 /*  triangle and an orientation into a single pointer.  My assumptions that  */
761 /*  triangles are four-byte-aligned and that the `unsigned long' type is     */
762 /*  long enough to hold a pointer are two of the few kludges in this program.*/
763 /*                                                                           */
764 /*  Shell edges are manipulated similarly.  A pointer to a shell edge        */
765 /*  carries both an address and an orientation in the range 0 to 1.          */
766 /*                                                                           */
767 /*  The other primitives take an oriented triangle or oriented shell edge,   */
768 /*  and return an oriented triangle or oriented shell edge or point; or they */
769 /*  change the connections in the data structure.                            */
770 /*                                                                           */
771 /*****************************************************************************/
772
773 /********* Mesh manipulation primitives begin here                   *********/
774 /**                                                                         **/
775 /**                                                                         **/
776
777 /* Fast lookup arrays to speed some of the mesh manipulation primitives.     */
778
779 int plus1mod3[3] = {1, 2, 0};
780 int minus1mod3[3] = {2, 0, 1};
781
782 /********* Primitives for triangles                                  *********/
783 /*                                                                           */
784 /*                                                                           */
785
786 /* decode() converts a pointer to an oriented triangle.  The orientation is  */
787 /*   extracted from the two least significant bits of the pointer.           */
788
789 #define decode(ptr, triedge)                                                  \
790   (triedge).orient = (int) ((unsigned long) (ptr) & (unsigned long) 3l);      \
791   (triedge).tri = (triangle *)                                                \
792                   ((unsigned long) (ptr) ^ (unsigned long) (triedge).orient)
793
794 /* encode() compresses an oriented triangle into a single pointer.  It       */
795 /*   relies on the assumption that all triangles are aligned to four-byte    */
796 /*   boundaries, so the two least significant bits of (triedge).tri are zero.*/
797
798 #define encode(triedge)                                                       \
799   (triangle) ((unsigned long) (triedge).tri | (unsigned long) (triedge).orient)
800
801 /* The following edge manipulation primitives are all described by Guibas    */
802 /*   and Stolfi.  However, they use an edge-based data structure, whereas I  */
803 /*   am using a triangle-based data structure.                               */
804
805 /* sym() finds the abutting triangle, on the same edge.  Note that the       */
806 /*   edge direction is necessarily reversed, because triangle/edge handles   */
807 /*   are always directed counterclockwise around the triangle.               */
808
809 #define sym(triedge1, triedge2)                                               \
810   ptr = (triedge1).tri[(triedge1).orient];                                    \
811   decode(ptr, triedge2);
812
813 #define symself(triedge)                                                      \
814   ptr = (triedge).tri[(triedge).orient];                                      \
815   decode(ptr, triedge);
816
817 /* lnext() finds the next edge (counterclockwise) of a triangle.             */
818
819 #define lnext(triedge1, triedge2)                                             \
820   (triedge2).tri = (triedge1).tri;                                            \
821   (triedge2).orient = plus1mod3[(triedge1).orient]
822
823 #define lnextself(triedge)                                                    \
824   (triedge).orient = plus1mod3[(triedge).orient]
825
826 /* lprev() finds the previous edge (clockwise) of a triangle.                */
827
828 #define lprev(triedge1, triedge2)                                             \
829   (triedge2).tri = (triedge1).tri;                                            \
830   (triedge2).orient = minus1mod3[(triedge1).orient]
831
832 #define lprevself(triedge)                                                    \
833   (triedge).orient = minus1mod3[(triedge).orient]
834
835 /* onext() spins counterclockwise around a point; that is, it finds the next */
836 /*   edge with the same origin in the counterclockwise direction.  This edge */
837 /*   will be part of a different triangle.                                   */
838
839 #define onext(triedge1, triedge2)                                             \
840   lprev(triedge1, triedge2);                                                  \
841   symself(triedge2);
842
843 #define onextself(triedge)                                                    \
844   lprevself(triedge);                                                         \
845   symself(triedge);
846
847 /* oprev() spins clockwise around a point; that is, it finds the next edge   */
848 /*   with the same origin in the clockwise direction.  This edge will be     */
849 /*   part of a different triangle.                                           */
850
851 #define oprev(triedge1, triedge2)                                             \
852   sym(triedge1, triedge2);                                                    \
853   lnextself(triedge2);
854
855 #define oprevself(triedge)                                                    \
856   symself(triedge);                                                           \
857   lnextself(triedge);
858
859 /* dnext() spins counterclockwise around a point; that is, it finds the next */
860 /*   edge with the same destination in the counterclockwise direction.  This */
861 /*   edge will be part of a different triangle.                              */
862
863 #define dnext(triedge1, triedge2)                                             \
864   sym(triedge1, triedge2);                                                    \
865   lprevself(triedge2);
866
867 #define dnextself(triedge)                                                    \
868   symself(triedge);                                                           \
869   lprevself(triedge);
870
871 /* dprev() spins clockwise around a point; that is, it finds the next edge   */
872 /*   with the same destination in the clockwise direction.  This edge will   */
873 /*   be part of a different triangle.                                        */
874
875 #define dprev(triedge1, triedge2)                                             \
876   lnext(triedge1, triedge2);                                                  \
877   symself(triedge2);
878
879 #define dprevself(triedge)                                                    \
880   lnextself(triedge);                                                         \
881   symself(triedge);
882
883 /* rnext() moves one edge counterclockwise about the adjacent triangle.      */
884 /*   (It's best understood by reading Guibas and Stolfi.  It involves        */
885 /*   changing triangles twice.)                                              */
886
887 #define rnext(triedge1, triedge2)                                             \
888   sym(triedge1, triedge2);                                                    \
889   lnextself(triedge2);                                                        \
890   symself(triedge2);
891
892 #define rnextself(triedge)                                                    \
893   symself(triedge);                                                           \
894   lnextself(triedge);                                                         \
895   symself(triedge);
896
897 /* rnext() moves one edge clockwise about the adjacent triangle.             */
898 /*   (It's best understood by reading Guibas and Stolfi.  It involves        */
899 /*   changing triangles twice.)                                              */
900
901 #define rprev(triedge1, triedge2)                                             \
902   sym(triedge1, triedge2);                                                    \
903   lprevself(triedge2);                                                        \
904   symself(triedge2);
905
906 #define rprevself(triedge)                                                    \
907   symself(triedge);                                                           \
908   lprevself(triedge);                                                         \
909   symself(triedge);
910
911 /* These primitives determine or set the origin, destination, or apex of a   */
912 /* triangle.                                                                 */
913
914 #define org(triedge, pointptr)                                                \
915   pointptr = (point) (triedge).tri[plus1mod3[(triedge).orient] + 3]
916
917 #define dest(triedge, pointptr)                                               \
918   pointptr = (point) (triedge).tri[minus1mod3[(triedge).orient] + 3]
919
920 #define apex(triedge, pointptr)                                               \
921   pointptr = (point) (triedge).tri[(triedge).orient + 3]
922
923 #define setorg(triedge, pointptr)                                             \
924   (triedge).tri[plus1mod3[(triedge).orient] + 3] = (triangle) pointptr
925
926 #define setdest(triedge, pointptr)                                            \
927   (triedge).tri[minus1mod3[(triedge).orient] + 3] = (triangle) pointptr
928
929 #define setapex(triedge, pointptr)                                            \
930   (triedge).tri[(triedge).orient + 3] = (triangle) pointptr
931
932 #define setvertices2null(triedge)                                             \
933   (triedge).tri[3] = (triangle) NULL;                                         \
934   (triedge).tri[4] = (triangle) NULL;                                         \
935   (triedge).tri[5] = (triangle) NULL;
936
937 /* Bond two triangles together.                                              */
938
939 #define bond(triedge1, triedge2)                                              \
940   (triedge1).tri[(triedge1).orient] = encode(triedge2);                       \
941   (triedge2).tri[(triedge2).orient] = encode(triedge1)
942
943 /* Dissolve a bond (from one side).  Note that the other triangle will still */
944 /*   think it's connected to this triangle.  Usually, however, the other     */
945 /*   triangle is being deleted entirely, or bonded to another triangle, so   */
946 /*   it doesn't matter.                                                      */
947
948 #define dissolve(triedge)                                                     \
949   (triedge).tri[(triedge).orient] = (triangle) dummytri
950
951 /* Copy a triangle/edge handle.                                              */
952
953 #define triedgecopy(triedge1, triedge2)                                       \
954   (triedge2).tri = (triedge1).tri;                                            \
955   (triedge2).orient = (triedge1).orient
956
957 /* Test for equality of triangle/edge handles.                               */
958
959 #define triedgeequal(triedge1, triedge2)                                      \
960   (((triedge1).tri == (triedge2).tri) &&                                      \
961    ((triedge1).orient == (triedge2).orient))
962
963 /* Primitives to infect or cure a triangle with the virus.  These rely on    */
964 /*   the assumption that all shell edges are aligned to four-byte boundaries.*/
965
966 #define infect(triedge)                                                       \
967   (triedge).tri[6] = (triangle)                                               \
968                      ((unsigned long) (triedge).tri[6] | (unsigned long) 2l)
969
970 #define uninfect(triedge)                                                     \
971   (triedge).tri[6] = (triangle)                                               \
972                      ((unsigned long) (triedge).tri[6] & ~ (unsigned long) 2l)
973
974 /* Test a triangle for viral infection.                                      */
975
976 #define infected(triedge)                                                     \
977   (((unsigned long) (triedge).tri[6] & (unsigned long) 2l) != 0)
978
979 /* Check or set a triangle's attributes.                                     */
980
981 #define elemattribute(triedge, attnum)                                        \
982   ((REAL *) (triedge).tri)[elemattribindex + (attnum)]
983
984 #define setelemattribute(triedge, attnum, value)                              \
985   ((REAL *) (triedge).tri)[elemattribindex + (attnum)] = (REAL)value
986
987 /* Check or set a triangle's maximum area bound.                             */
988
989 #define areabound(triedge)  ((REAL *) (triedge).tri)[areaboundindex]
990
991 #define setareabound(triedge, value)                                          \
992   ((REAL *) (triedge).tri)[areaboundindex] = (REAL)value
993
994 /********* Primitives for shell edges                                *********/
995 /*                                                                           */
996 /*                                                                           */
997
998 /* sdecode() converts a pointer to an oriented shell edge.  The orientation  */
999 /*   is extracted from the least significant bit of the pointer.  The two    */
1000 /*   least significant bits (one for orientation, one for viral infection)   */
1001 /*   are masked out to produce the real pointer.                             */
1002
1003 #define sdecode(sptr, edge)                                                   \
1004   (edge).shorient = (int) ((unsigned long) (sptr) & (unsigned long) 1l);      \
1005   (edge).sh = (shelle *)                                                      \
1006               ((unsigned long) (sptr) & ~ (unsigned long) 3l)
1007
1008 /* sencode() compresses an oriented shell edge into a single pointer.  It    */
1009 /*   relies on the assumption that all shell edges are aligned to two-byte   */
1010 /*   boundaries, so the least significant bit of (edge).sh is zero.          */
1011
1012 #define sencode(edge)                                                         \
1013   (shelle) ((unsigned long) (edge).sh | (unsigned long) (edge).shorient)
1014
1015 /* ssym() toggles the orientation of a shell edge.                           */
1016
1017 #define ssym(edge1, edge2)                                                    \
1018   (edge2).sh = (edge1).sh;                                                    \
1019   (edge2).shorient = 1 - (edge1).shorient
1020
1021 #define ssymself(edge)                                                        \
1022   (edge).shorient = 1 - (edge).shorient
1023
1024 /* spivot() finds the other shell edge (from the same segment) that shares   */
1025 /*   the same origin.                                                        */
1026
1027 #define spivot(edge1, edge2)                                                  \
1028   sptr = (edge1).sh[(edge1).shorient];                                        \
1029   sdecode(sptr, edge2)
1030
1031 #define spivotself(edge)                                                      \
1032   sptr = (edge).sh[(edge).shorient];                                          \
1033   sdecode(sptr, edge)
1034
1035 /* snext() finds the next shell edge (from the same segment) in sequence;    */
1036 /*   one whose origin is the input shell edge's destination.                 */
1037
1038 #define snext(edge1, edge2)                                                   \
1039   sptr = (edge1).sh[1 - (edge1).shorient];                                    \
1040   sdecode(sptr, edge2)
1041
1042 #define snextself(edge)                                                       \
1043   sptr = (edge).sh[1 - (edge).shorient];                                      \
1044   sdecode(sptr, edge)
1045
1046 /* These primitives determine or set the origin or destination of a shell    */
1047 /*   edge.                                                                   */
1048
1049 #define sorg(edge, pointptr)                                                  \
1050   pointptr = (point) (edge).sh[2 + (edge).shorient]
1051
1052 #define sdest(edge, pointptr)                                                 \
1053   pointptr = (point) (edge).sh[3 - (edge).shorient]
1054
1055 #define setsorg(edge, pointptr)                                               \
1056   (edge).sh[2 + (edge).shorient] = (shelle) pointptr
1057
1058 #define setsdest(edge, pointptr)                                              \
1059   (edge).sh[3 - (edge).shorient] = (shelle) pointptr
1060
1061 /* These primitives read or set a shell marker.  Shell markers are used to   */
1062 /*   hold user boundary information.                                         */
1063
1064 #define mark(edge)  (* (int *) ((edge).sh + 6))
1065
1066 #define setmark(edge, value)                                                  \
1067   * (int *) ((edge).sh + 6) = value
1068
1069 /* Bond two shell edges together.                                            */
1070
1071 #define sbond(edge1, edge2)                                                   \
1072   (edge1).sh[(edge1).shorient] = sencode(edge2);                              \
1073   (edge2).sh[(edge2).shorient] = sencode(edge1)
1074
1075 /* Dissolve a shell edge bond (from one side).  Note that the other shell    */
1076 /*   edge will still think it's connected to this shell edge.                */
1077
1078 #define sdissolve(edge)                                                       \
1079   (edge).sh[(edge).shorient] = (shelle) dummysh
1080
1081 /* Copy a shell edge.                                                        */
1082
1083 #define shellecopy(edge1, edge2)                                              \
1084   (edge2).sh = (edge1).sh;                                                    \
1085   (edge2).shorient = (edge1).shorient
1086
1087 /* Test for equality of shell edges.                                         */
1088
1089 #define shelleequal(edge1, edge2)                                             \
1090   (((edge1).sh == (edge2).sh) &&                                              \
1091    ((edge1).shorient == (edge2).shorient))
1092
1093 /********* Primitives for interacting triangles and shell edges      *********/
1094 /*                                                                           */
1095 /*                                                                           */
1096
1097 /* tspivot() finds a shell edge abutting a triangle.                         */
1098
1099 #define tspivot(triedge, edge)                                                \
1100   sptr = (shelle) (triedge).tri[6 + (triedge).orient];                        \
1101   sdecode(sptr, edge)
1102
1103 /* stpivot() finds a triangle abutting a shell edge.  It requires that the   */
1104 /*   variable `ptr' of type `triangle' be defined.                           */
1105
1106 #define stpivot(edge, triedge)                                                \
1107   ptr = (triangle) (edge).sh[4 + (edge).shorient];                            \
1108   decode(ptr, triedge)
1109
1110 /* Bond a triangle to a shell edge.                                          */
1111
1112 #define tsbond(triedge, edge)                                                 \
1113   (triedge).tri[6 + (triedge).orient] = (triangle) sencode(edge);             \
1114   (edge).sh[4 + (edge).shorient] = (shelle) encode(triedge)
1115
1116 /* Dissolve a bond (from the triangle side).                                 */
1117
1118 #define tsdissolve(triedge)                                                   \
1119   (triedge).tri[6 + (triedge).orient] = (triangle) dummysh
1120
1121 /* Dissolve a bond (from the shell edge side).                               */
1122
1123 #define stdissolve(edge)                                                      \
1124   (edge).sh[4 + (edge).shorient] = (shelle) dummytri
1125
1126 /********* Primitives for points                                     *********/
1127 /*                                                                           */
1128 /*                                                                           */
1129
1130 #define pointmark(pt)  ((int *) (pt))[pointmarkindex]
1131
1132 #define setpointmark(pt, value)                                               \
1133   ((int *) (pt))[pointmarkindex] = value
1134
1135 #define point2tri(pt)  ((triangle *) (pt))[point2triindex]
1136
1137 #define setpoint2tri(pt, value)                                               \
1138   ((triangle *) (pt))[point2triindex] = value
1139
1140 /**                                                                         **/
1141 /**                                                                         **/
1142 /********* Mesh manipulation primitives end here                     *********/
1143
1144 /********* User interaction routines begin here                      *********/
1145 /**                                                                         **/
1146 /**                                                                         **/
1147
1148 /*****************************************************************************/
1149 /*                                                                           */
1150 /*  syntax()   Print list of command line switches.                          */
1151 /*                                                                           */
1152 /*****************************************************************************/
1153
1154 #ifndef TRILIBRARY
1155
1156 void syntax()
1157 {
1158 #ifdef CDT_ONLY
1159 #ifdef REDUCED
1160   printf("triangle [-pAcevngBPNEIOXzo_lQVh] input_file\n");
1161 #else /* not REDUCED */
1162   printf("triangle [-pAcevngBPNEIOXzo_iFlCQVh] input_file\n");
1163 #endif /* not REDUCED */
1164 #else /* not CDT_ONLY */
1165 #ifdef REDUCED
1166   printf("triangle [-prq__a__AcevngBPNEIOXzo_YS__lQVh] input_file\n");
1167 #else /* not REDUCED */
1168   printf("triangle [-prq__a__AcevngBPNEIOXzo_YS__iFlsCQVh] input_file\n");
1169 #endif /* not REDUCED */
1170 #endif /* not CDT_ONLY */
1171
1172   printf("    -p  Triangulates a Planar Straight Line Graph (.poly file).\n");
1173 #ifndef CDT_ONLY
1174   printf("    -r  Refines a previously generated mesh.\n");
1175   printf(
1176     "    -q  Quality mesh generation.  A minimum angle may be specified.\n");
1177   printf("    -a  Applies a maximum triangle area constraint.\n");
1178 #endif /* not CDT_ONLY */
1179   printf(
1180     "    -A  Applies attributes to identify elements in certain regions.\n");
1181   printf("    -c  Encloses the convex hull with segments.\n");
1182   printf("    -e  Generates an edge list.\n");
1183   printf("    -v  Generates a Voronoi diagram.\n");
1184   printf("    -n  Generates a list of triangle neighbors.\n");
1185   printf("    -g  Generates an .off file for Geomview.\n");
1186   printf("    -B  Suppresses output of boundary information.\n");
1187   printf("    -P  Suppresses output of .poly file.\n");
1188   printf("    -N  Suppresses output of .node file.\n");
1189   printf("    -E  Suppresses output of .ele file.\n");
1190   printf("    -I  Suppresses mesh iteration numbers.\n");
1191   printf("    -O  Ignores holes in .poly file.\n");
1192   printf("    -X  Suppresses use of exact arithmetic.\n");
1193   printf("    -z  Numbers all items starting from zero (rather than one).\n");
1194   printf("    -o2 Generates second-order subparametric elements.\n");
1195 #ifndef CDT_ONLY
1196   printf("    -Y  Suppresses boundary segment splitting.\n");
1197   printf("    -S  Specifies maximum number of added Steiner points.\n");
1198 #endif /* not CDT_ONLY */
1199 #ifndef REDUCED
1200   printf("    -i  Uses incremental method, rather than divide-and-conquer.\n");
1201   printf("    -F  Uses Fortune's sweepline algorithm, rather than d-and-c.\n");
1202 #endif /* not REDUCED */
1203   printf("    -l  Uses vertical cuts only, rather than alternating cuts.\n");
1204 #ifndef REDUCED
1205 #ifndef CDT_ONLY
1206   printf(
1207     "    -s  Force segments into mesh by splitting (instead of using CDT).\n");
1208 #endif /* not CDT_ONLY */
1209   printf("    -C  Check consistency of final mesh.\n");
1210 #endif /* not REDUCED */
1211   printf("    -Q  Quiet:  No terminal output except errors.\n");
1212   printf("    -V  Verbose:  Detailed information on what I'm doing.\n");
1213   printf("    -h  Help:  Detailed instructions for Triangle.\n");
1214   exit(0);
1215 }
1216
1217 #endif /* not TRILIBRARY */
1218
1219 /*****************************************************************************/
1220 /*                                                                           */
1221 /*  info()   Print out complete instructions.                                */
1222 /*                                                                           */
1223 /*****************************************************************************/
1224
1225 #ifndef TRILIBRARY
1226
1227 void info()
1228 {
1229   printf("Triangle\n");
1230   printf(
1231 "A Two-Dimensional Quality Mesh Generator and Delaunay Triangulator.\n");
1232   printf("Version 1.3\n\n");
1233   printf(
1234 "Copyright 1996 Jonathan Richard Shewchuk  (bugs/comments to jrs@cs.cmu.edu)\n"
1235 );
1236   printf("School of Computer Science / Carnegie Mellon University\n");
1237   printf("5000 Forbes Avenue / Pittsburgh, Pennsylvania  15213-3891\n");
1238   printf(
1239 "Created as part of the Archimedes project (tools for parallel FEM).\n");
1240   printf(
1241 "Supported in part by NSF Grant CMS-9318163 and an NSERC 1967 Scholarship.\n");
1242   printf("There is no warranty whatsoever.  Use at your own risk.\n");
1243 #ifdef SINGLE
1244   printf("This executable is compiled for single precision arithmetic.\n\n\n");
1245 #else /* not SINGLE */
1246   printf("This executable is compiled for double precision arithmetic.\n\n\n");
1247 #endif /* not SINGLE */
1248   printf(
1249 "Triangle generates exact Delaunay triangulations, constrained Delaunay\n");
1250   printf(
1251 "triangulations, and quality conforming Delaunay triangulations.  The latter\n"
1252 );
1253   printf(
1254 "can be generated with no small angles, and are thus suitable for finite\n");
1255   printf(
1256 "element analysis.  If no command line switches are specified, your .node\n");
1257   printf(
1258 "input file will be read, and the Delaunay triangulation will be returned in\n"
1259 );
1260   printf(".node and .ele output files.  The command syntax is:\n\n");
1261 #ifdef CDT_ONLY
1262 #ifdef REDUCED
1263   printf("triangle [-pAcevngBPNEIOXzo_lQVh] input_file\n\n");
1264 #else /* not REDUCED */
1265   printf("triangle [-pAcevngBPNEIOXzo_iFlCQVh] input_file\n\n");
1266 #endif /* not REDUCED */
1267 #else /* not CDT_ONLY */
1268 #ifdef REDUCED
1269   printf("triangle [-prq__a__AcevngBPNEIOXzo_YS__lQVh] input_file\n\n");
1270 #else /* not REDUCED */
1271   printf("triangle [-prq__a__AcevngBPNEIOXzo_YS__iFlsCQVh] input_file\n\n");
1272 #endif /* not REDUCED */
1273 #endif /* not CDT_ONLY */
1274   printf(
1275 "Underscores indicate that numbers may optionally follow certain switches;\n");
1276   printf(
1277 "do not leave any space between a switch and its numeric parameter.\n");
1278   printf(
1279 "input_file must be a file with extension .node, or extension .poly if the\n");
1280   printf(
1281 "-p switch is used.  If -r is used, you must supply .node and .ele files,\n");
1282   printf(
1283 "and possibly a .poly file and .area file as well.  The formats of these\n");
1284   printf("files are described below.\n\n");
1285   printf("Command Line Switches:\n\n");
1286   printf(
1287 "    -p  Reads a Planar Straight Line Graph (.poly file), which can specify\n"
1288 );
1289   printf(
1290 "        points, segments, holes, and regional attributes and area\n");
1291   printf(
1292 "        constraints.  Will generate a constrained Delaunay triangulation\n");
1293   printf(
1294 "        fitting the input; or, if -s, -q, or -a is used, a conforming\n");
1295   printf(
1296 "        Delaunay triangulation.  If -p is not used, Triangle reads a .node\n"
1297 );
1298   printf("        file by default.\n");
1299   printf(
1300 "    -r  Refines a previously generated mesh.  The mesh is read from a .node\n"
1301 );
1302   printf(
1303 "        file and an .ele file.  If -p is also used, a .poly file is read\n");
1304   printf(
1305 "        and used to constrain edges in the mesh.  Further details on\n");
1306   printf("        refinement are given below.\n");
1307   printf(
1308 "    -q  Quality mesh generation by Jim Ruppert's Delaunay refinement\n");
1309   printf(
1310 "        algorithm.  Adds points to the mesh to ensure that no angles\n");
1311   printf(
1312 "        smaller than 20 degrees occur.  An alternative minimum angle may be\n"
1313 );
1314   printf(
1315 "        specified after the `q'.  If the minimum angle is 20.7 degrees or\n");
1316   printf(
1317 "        smaller, the triangulation algorithm is theoretically guaranteed to\n"
1318 );
1319   printf(
1320 "        terminate (assuming infinite precision arithmetic - Triangle may\n");
1321   printf(
1322 "        fail to terminate if you run out of precision).  In practice, the\n");
1323   printf(
1324 "        algorithm often succeeds for minimum angles up to 33.8 degrees.\n");
1325   printf(
1326 "        For highly refined meshes, however, it may be necessary to reduce\n");
1327   printf(
1328 "        the minimum angle to well below 20 to avoid problems associated\n");
1329   printf(
1330 "        with insufficient floating-point precision.  The specified angle\n");
1331   printf("        may include a decimal point.\n");
1332   printf(
1333 "    -a  Imposes a maximum triangle area.  If a number follows the `a', no\n");
1334   printf(
1335 "        triangle will be generated whose area is larger than that number.\n");
1336   printf(
1337 "        If no number is specified, an .area file (if -r is used) or .poly\n");
1338   printf(
1339 "        file (if -r is not used) specifies a number of maximum area\n");
1340   printf(
1341 "        constraints.  An .area file contains a separate area constraint for\n"
1342 );
1343   printf(
1344 "        each triangle, and is useful for refining a finite element mesh\n");
1345   printf(
1346 "        based on a posteriori error estimates.  A .poly file can optionally\n"
1347 );
1348   printf(
1349 "        contain an area constraint for each segment-bounded region, thereby\n"
1350 );
1351   printf(
1352 "        enforcing triangle densities in a first triangulation.  You can\n");
1353   printf(
1354 "        impose both a fixed area constraint and a varying area constraint\n");
1355   printf(
1356 "        by invoking the -a switch twice, once with and once without a\n");
1357   printf(
1358 "        number following.  Each area specified may include a decimal point.\n"
1359 );
1360   printf(
1361 "    -A  Assigns an additional attribute to each triangle that identifies\n");
1362   printf(
1363 "        what segment-bounded region each triangle belongs to.  Attributes\n");
1364   printf(
1365 "        are assigned to regions by the .poly file.  If a region is not\n");
1366   printf(
1367 "        explicitly marked by the .poly file, triangles in that region are\n");
1368   printf(
1369 "        assigned an attribute of zero.  The -A switch has an effect only\n");
1370   printf("        when the -p switch is used and the -r switch is not.\n");
1371   printf(
1372 "    -c  Creates segments on the convex hull of the triangulation.  If you\n");
1373   printf(
1374 "        are triangulating a point set, this switch causes a .poly file to\n");
1375   printf(
1376 "        be written, containing all edges in the convex hull.  (By default,\n"
1377 );
1378   printf(
1379 "        a .poly file is written only if a .poly file is read.)  If you are\n"
1380 );
1381   printf(
1382 "        triangulating a PSLG, this switch specifies that the interior of\n");
1383   printf(
1384 "        the convex hull of the PSLG should be triangulated.  If you do not\n"
1385 );
1386   printf(
1387 "        use this switch when triangulating a PSLG, it is assumed that you\n");
1388   printf(
1389 "        have identified the region to be triangulated by surrounding it\n");
1390   printf(
1391 "        with segments of the input PSLG.  Beware:  if you are not careful,\n"
1392 );
1393   printf(
1394 "        this switch can cause the introduction of an extremely thin angle\n");
1395   printf(
1396 "        between a PSLG segment and a convex hull segment, which can cause\n");
1397   printf(
1398 "        overrefinement or failure if Triangle runs out of precision.  If\n");
1399   printf(
1400 "        you are refining a mesh, the -c switch works differently; it\n");
1401   printf(
1402 "        generates the set of boundary edges of the mesh, rather than the\n");
1403   printf("        convex hull.\n");
1404   printf(
1405 "    -e  Outputs (to an .edge file) a list of edges of the triangulation.\n");
1406   printf(
1407 "    -v  Outputs the Voronoi diagram associated with the triangulation.\n");
1408   printf("        Does not attempt to detect degeneracies.\n");
1409   printf(
1410 "    -n  Outputs (to a .neigh file) a list of triangles neighboring each\n");
1411   printf("        triangle.\n");
1412   printf(
1413 "    -g  Outputs the mesh to an Object File Format (.off) file, suitable for\n"
1414 );
1415   printf("        viewing with the Geometry Center's Geomview package.\n");
1416   printf(
1417 "    -B  No boundary markers in the output .node, .poly, and .edge output\n");
1418   printf(
1419 "        files.  See the detailed discussion of boundary markers below.\n");
1420   printf(
1421 "    -P  No output .poly file.  Saves disk space, but you lose the ability\n");
1422   printf(
1423 "        to impose segment constraints on later refinements of the mesh.\n");
1424   printf("    -N  No output .node file.\n");
1425   printf("    -E  No output .ele file.\n");
1426   printf(
1427 "    -I  No iteration numbers.  Suppresses the output of .node and .poly\n");
1428   printf(
1429 "        files, so your input files won't be overwritten.  (If your input is\n"
1430 );
1431   printf(
1432 "        a .poly file only, a .node file will be written.)  Cannot be used\n");
1433   printf(
1434 "        with the -r switch, because that would overwrite your input .ele\n");
1435   printf(
1436 "        file.  Shouldn't be used with the -s, -q, or -a switch if you are\n");
1437   printf(
1438 "        using a .node file for input, because no .node file will be\n");
1439   printf("        written, so there will be no record of any added points.\n");
1440   printf("    -O  No holes.  Ignores the holes in the .poly file.\n");
1441   printf(
1442 "    -X  No exact arithmetic.  Normally, Triangle uses exact floating-point\n"
1443 );
1444   printf(
1445 "        arithmetic for certain tests if it thinks the inexact tests are not\n"
1446 );
1447   printf(
1448 "        accurate enough.  Exact arithmetic ensures the robustness of the\n");
1449   printf(
1450 "        triangulation algorithms, despite floating-point roundoff error.\n");
1451   printf(
1452 "        Disabling exact arithmetic with the -X switch will cause a small\n");
1453   printf(
1454 "        improvement in speed and create the possibility (albeit small) that\n"
1455 );
1456   printf(
1457 "        Triangle will fail to produce a valid mesh.  Not recommended.\n");
1458   printf(
1459 "    -z  Numbers all items starting from zero (rather than one).  Note that\n"
1460 );
1461   printf(
1462 "        this switch is normally overrided by the value used to number the\n");
1463   printf(
1464 "        first point of the input .node or .poly file.  However, this switch\n"
1465 );
1466   printf("        is useful when calling Triangle from another program.\n");
1467   printf(
1468 "    -o2 Generates second-order subparametric elements with six nodes each.\n"
1469 );
1470   printf(
1471 "    -Y  No new points on the boundary.  This switch is useful when the mesh\n"
1472 );
1473   printf(
1474 "        boundary must be preserved so that it conforms to some adjacent\n");
1475   printf(
1476 "        mesh.  Be forewarned that you will probably sacrifice some of the\n");
1477   printf(
1478 "        quality of the mesh; Triangle will try, but the resulting mesh may\n"
1479 );
1480   printf(
1481 "        contain triangles of poor aspect ratio.  Works well if all the\n");
1482   printf(
1483 "        boundary points are closely spaced.  Specify this switch twice\n");
1484   printf(
1485 "        (`-YY') to prevent all segment splitting, including internal\n");
1486   printf("        boundaries.\n");
1487   printf(
1488 "    -S  Specifies the maximum number of Steiner points (points that are not\n"
1489 );
1490   printf(
1491 "        in the input, but are added to meet the constraints of minimum\n");
1492   printf(
1493 "        angle and maximum area).  The default is to allow an unlimited\n");
1494   printf(
1495 "        number.  If you specify this switch with no number after it,\n");
1496   printf(
1497 "        the limit is set to zero.  Triangle always adds points at segment\n");
1498   printf(
1499 "        intersections, even if it needs to use more points than the limit\n");
1500   printf(
1501 "        you set.  When Triangle inserts segments by splitting (-s), it\n");
1502   printf(
1503 "        always adds enough points to ensure that all the segments appear in\n"
1504 );
1505   printf(
1506 "        the triangulation, again ignoring the limit.  Be forewarned that\n");
1507   printf(
1508 "        the -S switch may result in a conforming triangulation that is not\n"
1509 );
1510   printf(
1511 "        truly Delaunay, because Triangle may be forced to stop adding\n");
1512   printf(
1513 "        points when the mesh is in a state where a segment is non-Delaunay\n"
1514 );
1515   printf(
1516 "        and needs to be split.  If so, Triangle will print a warning.\n");
1517   printf(
1518 "    -i  Uses an incremental rather than divide-and-conquer algorithm to\n");
1519   printf(
1520 "        form a Delaunay triangulation.  Try it if the divide-and-conquer\n");
1521   printf("        algorithm fails.\n");
1522   printf(
1523 "    -F  Uses Steven Fortune's sweepline algorithm to form a Delaunay\n");
1524   printf(
1525 "        triangulation.  Warning:  does not use exact arithmetic for all\n");
1526   printf("        calculations.  An exact result is not guaranteed.\n");
1527   printf(
1528 "    -l  Uses only vertical cuts in the divide-and-conquer algorithm.  By\n");
1529   printf(
1530 "        default, Triangle uses alternating vertical and horizontal cuts,\n");
1531   printf(
1532 "        which usually improve the speed except with point sets that are\n");
1533   printf(
1534 "        small or short and wide.  This switch is primarily of theoretical\n");
1535   printf("        interest.\n");
1536   printf(
1537 "    -s  Specifies that segments should be forced into the triangulation by\n"
1538 );
1539   printf(
1540 "        recursively splitting them at their midpoints, rather than by\n");
1541   printf(
1542 "        generating a constrained Delaunay triangulation.  Segment splitting\n"
1543 );
1544   printf(
1545 "        is true to Ruppert's original algorithm, but can create needlessly\n"
1546 );
1547   printf("        small triangles near external small features.\n");
1548   printf(
1549 "    -C  Check the consistency of the final mesh.  Uses exact arithmetic for\n"
1550 );
1551   printf(
1552 "        checking, even if the -X switch is used.  Useful if you suspect\n");
1553   printf("        Triangle is buggy.\n");
1554   printf(
1555 "    -Q  Quiet: Suppresses all explanation of what Triangle is doing, unless\n"
1556 );
1557   printf("        an error occurs.\n");
1558   printf(
1559 "    -V  Verbose: Gives detailed information about what Triangle is doing.\n");
1560   printf(
1561 "        Add more `V's for increasing amount of detail.  `-V' gives\n");
1562   printf(
1563 "        information on algorithmic progress and more detailed statistics.\n");
1564   printf(
1565 "        `-VV' gives point-by-point details, and will print so much that\n");
1566   printf(
1567 "        Triangle will run much more slowly.  `-VVV' gives information only\n"
1568 );
1569   printf("        a debugger could love.\n");
1570   printf("    -h  Help:  Displays these instructions.\n");
1571   printf("\n");
1572   printf("Definitions:\n");
1573   printf("\n");
1574   printf(
1575 "  A Delaunay triangulation of a point set is a triangulation whose vertices\n"
1576 );
1577   printf(
1578 "  are the point set, having the property that no point in the point set\n");
1579   printf(
1580 "  falls in the interior of the circumcircle (circle that passes through all\n"
1581 );
1582   printf("  three vertices) of any triangle in the triangulation.\n\n");
1583   printf(
1584 "  A Voronoi diagram of a point set is a subdivision of the plane into\n");
1585   printf(
1586 "  polygonal regions (some of which may be infinite), where each region is\n");
1587   printf(
1588 "  the set of points in the plane that are closer to some input point than\n");
1589   printf(
1590 "  to any other input point.  (The Voronoi diagram is the geometric dual of\n"
1591 );
1592   printf("  the Delaunay triangulation.)\n\n");
1593   printf(
1594 "  A Planar Straight Line Graph (PSLG) is a collection of points and\n");
1595   printf(
1596 "  segments.  Segments are simply edges, whose endpoints are points in the\n");
1597   printf(
1598 "  PSLG.  The file format for PSLGs (.poly files) is described below.\n");
1599   printf("\n");
1600   printf(
1601 "  A constrained Delaunay triangulation of a PSLG is similar to a Delaunay\n");
1602   printf(
1603 "  triangulation, but each PSLG segment is present as a single edge in the\n");
1604   printf(
1605 "  triangulation.  (A constrained Delaunay triangulation is not truly a\n");
1606   printf("  Delaunay triangulation.)\n\n");
1607   printf(
1608 "  A conforming Delaunay triangulation of a PSLG is a true Delaunay\n");
1609   printf(
1610 "  triangulation in which each PSLG segment may have been subdivided into\n");
1611   printf(
1612 "  several edges by the insertion of additional points.  These inserted\n");
1613   printf(
1614 "  points are necessary to allow the segments to exist in the mesh while\n");
1615   printf("  maintaining the Delaunay property.\n\n");
1616   printf("File Formats:\n\n");
1617   printf(
1618 "  All files may contain comments prefixed by the character '#'.  Points,\n");
1619   printf(
1620 "  triangles, edges, holes, and maximum area constraints must be numbered\n");
1621   printf(
1622 "  consecutively, starting from either 1 or 0.  Whichever you choose, all\n");
1623   printf(
1624 "  input files must be consistent; if the nodes are numbered from 1, so must\n"
1625 );
1626   printf(
1627 "  be all other objects.  Triangle automatically detects your choice while\n");
1628   printf(
1629 "  reading the .node (or .poly) file.  (When calling Triangle from another\n");
1630   printf(
1631 "  program, use the -z switch if you wish to number objects from zero.)\n");
1632   printf("  Examples of these file formats are given below.\n\n");
1633   printf("  .node files:\n");
1634   printf(
1635 "    First line:  <# of points> <dimension (must be 2)> <# of attributes>\n");
1636   printf(
1637 "                                           <# of boundary markers (0 or 1)>\n"
1638 );
1639   printf(
1640 "    Remaining lines:  <point #> <x> <y> [attributes] [boundary marker]\n");
1641   printf("\n");
1642   printf(
1643 "    The attributes, which are typically floating-point values of physical\n");
1644   printf(
1645 "    quantities (such as mass or conductivity) associated with the nodes of\n"
1646 );
1647   printf(
1648 "    a finite element mesh, are copied unchanged to the output mesh.  If -s,\n"
1649 );
1650   printf(
1651 "    -q, or -a is selected, each new Steiner point added to the mesh will\n");
1652   printf("    have attributes assigned to it by linear interpolation.\n\n");
1653   printf(
1654 "    If the fourth entry of the first line is `1', the last column of the\n");
1655   printf(
1656 "    remainder of the file is assumed to contain boundary markers.  Boundary\n"
1657 );
1658   printf(
1659 "    markers are used to identify boundary points and points resting on PSLG\n"
1660 );
1661   printf(
1662 "    segments; a complete description appears in a section below.  The .node\n"
1663 );
1664   printf(
1665 "    file produced by Triangle will contain boundary markers in the last\n");
1666   printf("    column unless they are suppressed by the -B switch.\n\n");
1667   printf("  .ele files:\n");
1668   printf(
1669 "    First line:  <# of triangles> <points per triangle> <# of attributes>\n");
1670   printf(
1671 "    Remaining lines:  <triangle #> <point> <point> <point> ... [attributes]\n"
1672 );
1673   printf("\n");
1674   printf(
1675 "    Points are indices into the corresponding .node file.  The first three\n"
1676 );
1677   printf(
1678 "    points are the corners, and are listed in counterclockwise order around\n"
1679 );
1680   printf(
1681 "    each triangle.  (The remaining points, if any, depend on the type of\n");
1682   printf(
1683 "    finite element used.)  The attributes are just like those of .node\n");
1684   printf(
1685 "    files.  Because there is no simple mapping from input to output\n");
1686   printf(
1687 "    triangles, an attempt is made to interpolate attributes, which may\n");
1688   printf(
1689 "    result in a good deal of diffusion of attributes among nearby triangles\n"
1690 );
1691   printf(
1692 "    as the triangulation is refined.  Diffusion does not occur across\n");
1693   printf(
1694 "    segments, so attributes used to identify segment-bounded regions remain\n"
1695 );
1696   printf(
1697 "    intact.  In output .ele files, all triangles have three points each\n");
1698   printf(
1699 "    unless the -o2 switch is used, in which case they have six, and the\n");
1700   printf(
1701 "    fourth, fifth, and sixth points lie on the midpoints of the edges\n");
1702   printf("    opposite the first, second, and third corners.\n\n");
1703   printf("  .poly files:\n");
1704   printf(
1705 "    First line:  <# of points> <dimension (must be 2)> <# of attributes>\n");
1706   printf(
1707 "                                           <# of boundary markers (0 or 1)>\n"
1708 );
1709   printf(
1710 "    Following lines:  <point #> <x> <y> [attributes] [boundary marker]\n");
1711   printf("    One line:  <# of segments> <# of boundary markers (0 or 1)>\n");
1712   printf(
1713 "    Following lines:  <segment #> <endpoint> <endpoint> [boundary marker]\n");
1714   printf("    One line:  <# of holes>\n");
1715   printf("    Following lines:  <hole #> <x> <y>\n");
1716   printf(
1717 "    Optional line:  <# of regional attributes and/or area constraints>\n");
1718   printf(
1719 "    Optional following lines:  <constraint #> <x> <y> <attrib> <max area>\n");
1720   printf("\n");
1721   printf(
1722 "    A .poly file represents a PSLG, as well as some additional information.\n"
1723 );
1724   printf(
1725 "    The first section lists all the points, and is identical to the format\n"
1726 );
1727   printf(
1728 "    of .node files.  <# of points> may be set to zero to indicate that the\n"
1729 );
1730   printf(
1731 "    points are listed in a separate .node file; .poly files produced by\n");
1732   printf(
1733 "    Triangle always have this format.  This has the advantage that a point\n"
1734 );
1735   printf(
1736 "    set may easily be triangulated with or without segments.  (The same\n");
1737   printf(
1738 "    effect can be achieved, albeit using more disk space, by making a copy\n"
1739 );
1740   printf(
1741 "    of the .poly file with the extension .node; all sections of the file\n");
1742   printf("    but the first are ignored.)\n\n");
1743   printf(
1744 "    The second section lists the segments.  Segments are edges whose\n");
1745   printf(
1746 "    presence in the triangulation is enforced.  Each segment is specified\n");
1747   printf(
1748 "    by listing the indices of its two endpoints.  This means that you must\n"
1749 );
1750   printf(
1751 "    include its endpoints in the point list.  If -s, -q, and -a are not\n");
1752   printf(
1753 "    selected, Triangle will produce a constrained Delaunay triangulation,\n");
1754   printf(
1755 "    in which each segment appears as a single edge in the triangulation.\n");
1756   printf(
1757 "    If -q or -a is selected, Triangle will produce a conforming Delaunay\n");
1758   printf(
1759 "    triangulation, in which segments may be subdivided into smaller edges.\n"
1760 );
1761   printf("    Each segment, like each point, may have a boundary marker.\n\n");
1762   printf(
1763 "    The third section lists holes (and concavities, if -c is selected) in\n");
1764   printf(
1765 "    the triangulation.  Holes are specified by identifying a point inside\n");
1766   printf(
1767 "    each hole.  After the triangulation is formed, Triangle creates holes\n");
1768   printf(
1769 "    by eating triangles, spreading out from each hole point until its\n");
1770   printf(
1771 "    progress is blocked by PSLG segments; you must be careful to enclose\n");
1772   printf(
1773 "    each hole in segments, or your whole triangulation may be eaten away.\n");
1774   printf(
1775 "    If the two triangles abutting a segment are eaten, the segment itself\n");
1776   printf(
1777 "    is also eaten.  Do not place a hole directly on a segment; if you do,\n");
1778   printf("    Triangle will choose one side of the segment arbitrarily.\n\n");
1779   printf(
1780 "    The optional fourth section lists regional attributes (to be assigned\n");
1781   printf(
1782 "    to all triangles in a region) and regional constraints on the maximum\n");
1783   printf(
1784 "    triangle area.  Triangle will read this section only if the -A switch\n");
1785   printf(
1786 "    is used or the -a switch is used without a number following it, and the\n"
1787 );
1788   printf(
1789 "    -r switch is not used.  Regional attributes and area constraints are\n");
1790   printf(
1791 "    propagated in the same manner as holes; you specify a point for each\n");
1792   printf(
1793 "    attribute and/or constraint, and the attribute and/or constraint will\n");
1794   printf(
1795 "    affect the whole region (bounded by segments) containing the point.  If\n"
1796 );
1797   printf(
1798 "    two values are written on a line after the x and y coordinate, the\n");
1799   printf(
1800 "    former is assumed to be a regional attribute (but will only be applied\n"
1801 );
1802   printf(
1803 "    if the -A switch is selected), and the latter is assumed to be a\n");
1804   printf(
1805 "    regional area constraint (but will only be applied if the -a switch is\n"
1806 );
1807   printf(
1808 "    selected).  You may also specify just one value after the coordinates,\n"
1809 );
1810   printf(
1811 "    which can serve as both an attribute and an area constraint, depending\n"
1812 );
1813   printf(
1814 "    on the choice of switches.  If you are using the -A and -a switches\n");
1815   printf(
1816 "    simultaneously and wish to assign an attribute to some region without\n");
1817   printf("    imposing an area constraint, use a negative maximum area.\n\n");
1818   printf(
1819 "    When a triangulation is created from a .poly file, you must either\n");
1820   printf(
1821 "    enclose the entire region to be triangulated in PSLG segments, or\n");
1822   printf(
1823 "    use the -c switch, which encloses the convex hull of the input point\n");
1824   printf(
1825 "    set.  If you do not use the -c switch, Triangle will eat all triangles\n"
1826 );
1827   printf(
1828 "    on the outer boundary that are not protected by segments; if you are\n");
1829   printf(
1830 "    not careful, your whole triangulation may be eaten away.  If you do\n");
1831   printf(
1832 "    use the -c switch, you can still produce concavities by appropriate\n");
1833   printf("    placement of holes just inside the convex hull.\n\n");
1834   printf(
1835 "    An ideal PSLG has no intersecting segments, nor any points that lie\n");
1836   printf(
1837 "    upon segments (except, of course, the endpoints of each segment.)  You\n"
1838 );
1839   printf(
1840 "    aren't required to make your .poly files ideal, but you should be aware\n"
1841 );
1842   printf(
1843 "    of what can go wrong.  Segment intersections are relatively safe -\n");
1844   printf(
1845 "    Triangle will calculate the intersection points for you and add them to\n"
1846 );
1847   printf(
1848 "    the triangulation - as long as your machine's floating-point precision\n"
1849 );
1850   printf(
1851 "    doesn't become a problem.  You are tempting the fates if you have three\n"
1852 );
1853   printf(
1854 "    segments that cross at the same location, and expect Triangle to figure\n"
1855 );
1856   printf(
1857 "    out where the intersection point is.  Thanks to floating-point roundoff\n"
1858 );
1859   printf(
1860 "    error, Triangle will probably decide that the three segments intersect\n"
1861 );
1862   printf(
1863 "    at three different points, and you will find a minuscule triangle in\n");
1864   printf(
1865 "    your output - unless Triangle tries to refine the tiny triangle, uses\n");
1866   printf(
1867 "    up the last bit of machine precision, and fails to terminate at all.\n");
1868   printf(
1869 "    You're better off putting the intersection point in the input files,\n");
1870   printf(
1871 "    and manually breaking up each segment into two.  Similarly, if you\n");
1872   printf(
1873 "    place a point at the middle of a segment, and hope that Triangle will\n");
1874   printf(
1875 "    break up the segment at that point, you might get lucky.  On the other\n"
1876 );
1877   printf(
1878 "    hand, Triangle might decide that the point doesn't lie precisely on the\n"
1879 );
1880   printf(
1881 "    line, and you'll have a needle-sharp triangle in your output - or a lot\n"
1882 );
1883   printf("    of tiny triangles if you're generating a quality mesh.\n\n");
1884   printf(
1885 "    When Triangle reads a .poly file, it also writes a .poly file, which\n");
1886   printf(
1887 "    includes all edges that are part of input segments.  If the -c switch\n");
1888   printf(
1889 "    is used, the output .poly file will also include all of the edges on\n");
1890   printf(
1891 "    the convex hull.  Hence, the output .poly file is useful for finding\n");
1892   printf(
1893 "    edges associated with input segments and setting boundary conditions in\n"
1894 );
1895   printf(
1896 "    finite element simulations.  More importantly, you will need it if you\n"
1897 );
1898   printf(
1899 "    plan to refine the output mesh, and don't want segments to be missing\n");
1900   printf("    in later triangulations.\n\n");
1901   printf("  .area files:\n");
1902   printf("    First line:  <# of triangles>\n");
1903   printf("    Following lines:  <triangle #> <maximum area>\n\n");
1904   printf(
1905 "    An .area file associates with each triangle a maximum area that is used\n"
1906 );
1907   printf(
1908 "    for mesh refinement.  As with other file formats, every triangle must\n");
1909   printf(
1910 "    be represented, and they must be numbered consecutively.  A triangle\n");
1911   printf(
1912 "    may be left unconstrained by assigning it a negative maximum area.\n");
1913   printf("\n");
1914   printf("  .edge files:\n");
1915   printf("    First line:  <# of edges> <# of boundary markers (0 or 1)>\n");
1916   printf(
1917 "    Following lines:  <edge #> <endpoint> <endpoint> [boundary marker]\n");
1918   printf("\n");
1919   printf(
1920 "    Endpoints are indices into the corresponding .node file.  Triangle can\n"
1921 );
1922   printf(
1923 "    produce .edge files (use the -e switch), but cannot read them.  The\n");
1924   printf(
1925 "    optional column of boundary markers is suppressed by the -B switch.\n");
1926   printf("\n");
1927   printf(
1928 "    In Voronoi diagrams, one also finds a special kind of edge that is an\n");
1929   printf(
1930 "    infinite ray with only one endpoint.  For these edges, a different\n");
1931   printf("    format is used:\n\n");
1932   printf("        <edge #> <endpoint> -1 <direction x> <direction y>\n\n");
1933   printf(
1934 "    The `direction' is a floating-point vector that indicates the direction\n"
1935 );
1936   printf("    of the infinite ray.\n\n");
1937   printf("  .neigh files:\n");
1938   printf(
1939 "    First line:  <# of triangles> <# of neighbors per triangle (always 3)>\n"
1940 );
1941   printf(
1942 "    Following lines:  <triangle #> <neighbor> <neighbor> <neighbor>\n");
1943   printf("\n");
1944   printf(
1945 "    Neighbors are indices into the corresponding .ele file.  An index of -1\n"
1946 );
1947   printf(
1948 "    indicates a mesh boundary, and therefore no neighbor.  Triangle can\n");
1949   printf(
1950 "    produce .neigh files (use the -n switch), but cannot read them.\n");
1951   printf("\n");
1952   printf(
1953 "    The first neighbor of triangle i is opposite the first corner of\n");
1954   printf("    triangle i, and so on.\n\n");
1955   printf("Boundary Markers:\n\n");
1956   printf(
1957 "  Boundary markers are tags used mainly to identify which output points and\n"
1958 );
1959   printf(
1960 "  edges are associated with which PSLG segment, and to identify which\n");
1961   printf(
1962 "  points and edges occur on a boundary of the triangulation.  A common use\n"
1963 );
1964   printf(
1965 "  is to determine where boundary conditions should be applied to a finite\n");
1966   printf(
1967 "  element mesh.  You can prevent boundary markers from being written into\n");
1968   printf("  files produced by Triangle by using the -B switch.\n\n");
1969   printf(
1970 "  The boundary marker associated with each segment in an output .poly file\n"
1971 );
1972   printf("  or edge in an output .edge file is chosen as follows:\n");
1973   printf(
1974 "    - If an output edge is part or all of a PSLG segment with a nonzero\n");
1975   printf(
1976 "      boundary marker, then the edge is assigned the same marker.\n");
1977   printf(
1978 "    - Otherwise, if the edge occurs on a boundary of the triangulation\n");
1979   printf(
1980 "      (including boundaries of holes), then the edge is assigned the marker\n"
1981 );
1982   printf("      one (1).\n");
1983   printf("    - Otherwise, the edge is assigned the marker zero (0).\n");
1984   printf(
1985 "  The boundary marker associated with each point in an output .node file is\n"
1986 );
1987   printf("  chosen as follows:\n");
1988   printf(
1989 "    - If a point is assigned a nonzero boundary marker in the input file,\n");
1990   printf(
1991 "      then it is assigned the same marker in the output .node file.\n");
1992   printf(
1993 "    - Otherwise, if the point lies on a PSLG segment (including the\n");
1994   printf(
1995 "      segment's endpoints) with a nonzero boundary marker, then the point\n");
1996   printf(
1997 "      is assigned the same marker.  If the point lies on several such\n");
1998   printf("      segments, one of the markers is chosen arbitrarily.\n");
1999   printf(
2000 "    - Otherwise, if the point occurs on a boundary of the triangulation,\n");
2001   printf("      then the point is assigned the marker one (1).\n");
2002   printf("    - Otherwise, the point is assigned the marker zero (0).\n");
2003   printf("\n");
2004   printf(
2005 "  If you want Triangle to determine for you which points and edges are on\n");
2006   printf(
2007 "  the boundary, assign them the boundary marker zero (or use no markers at\n"
2008 );
2009   printf(
2010 "  all) in your input files.  Alternatively, you can mark some of them and\n");
2011   printf("  leave others marked zero, allowing Triangle to label them.\n\n");
2012   printf("Triangulation Iteration Numbers:\n\n");
2013   printf(
2014 "  Because Triangle can read and refine its own triangulations, input\n");
2015   printf(
2016 "  and output files have iteration numbers.  For instance, Triangle might\n");
2017   printf(
2018 "  read the files mesh.3.node, mesh.3.ele, and mesh.3.poly, refine the\n");
2019   printf(
2020 "  triangulation, and output the files mesh.4.node, mesh.4.ele, and\n");
2021   printf("  mesh.4.poly.  Files with no iteration number are treated as if\n");
2022   printf(
2023 "  their iteration number is zero; hence, Triangle might read the file\n");
2024   printf(
2025 "  points.node, triangulate it, and produce the files points.1.node and\n");
2026   printf("  points.1.ele.\n\n");
2027   printf(
2028 "  Iteration numbers allow you to create a sequence of successively finer\n");
2029   printf(
2030 "  meshes suitable for multigrid methods.  They also allow you to produce a\n"
2031 );
2032   printf(
2033 "  sequence of meshes using error estimate-driven mesh refinement.\n");
2034   printf("\n");
2035   printf(
2036 "  If you're not using refinement or quality meshing, and you don't like\n");
2037   printf(
2038 "  iteration numbers, use the -I switch to disable them.  This switch will\n");
2039   printf(
2040 "  also disable output of .node and .poly files to prevent your input files\n"
2041 );
2042   printf(
2043 "  from being overwritten.  (If the input is a .poly file that contains its\n"
2044 );
2045   printf("  own points, a .node file will be written.)\n\n");
2046   printf("Examples of How to Use Triangle:\n\n");
2047   printf(
2048 "  `triangle dots' will read points from dots.node, and write their Delaunay\n"
2049 );
2050   printf(
2051 "  triangulation to dots.1.node and dots.1.ele.  (dots.1.node will be\n");
2052   printf(
2053 "  identical to dots.node.)  `triangle -I dots' writes the triangulation to\n"
2054 );
2055   printf(
2056 "  dots.ele instead.  (No additional .node file is needed, so none is\n");
2057   printf("  written.)\n\n");
2058   printf(
2059 "  `triangle -pe object.1' will read a PSLG from object.1.poly (and possibly\n"
2060 );
2061   printf(
2062 "  object.1.node, if the points are omitted from object.1.poly) and write\n");
2063   printf("  their constrained Delaunay triangulation to object.2.node and\n");
2064   printf(
2065 "  object.2.ele.  The segments will be copied to object.2.poly, and all\n");
2066   printf("  edges will be written to object.2.edge.\n\n");
2067   printf(
2068 "  `triangle -pq31.5a.1 object' will read a PSLG from object.poly (and\n");
2069   printf(
2070 "  possibly object.node), generate a mesh whose angles are all greater than\n"
2071 );
2072   printf(
2073 "  31.5 degrees and whose triangles all have area smaller than 0.1, and\n");
2074   printf(
2075 "  write the mesh to object.1.node and object.1.ele.  Each segment may have\n"
2076 );
2077   printf(
2078 "  been broken up into multiple edges; the resulting constrained edges are\n");
2079   printf("  written to object.1.poly.\n\n");
2080   printf(
2081 "  Here is a sample file `box.poly' describing a square with a square hole:\n"
2082 );
2083   printf("\n");
2084   printf(
2085 "    # A box with eight points in 2D, no attributes, one boundary marker.\n");
2086   printf("    8 2 0 1\n");
2087   printf("    # Outer box has these vertices:\n");
2088   printf("     1   0 0   0\n");
2089   printf("     2   0 3   0\n");
2090   printf("     3   3 0   0\n");
2091   printf("     4   3 3   33     # A special marker for this point.\n");
2092   printf("    # Inner square has these vertices:\n");
2093   printf("     5   1 1   0\n");
2094   printf("     6   1 2   0\n");
2095   printf("     7   2 1   0\n");
2096   printf("     8   2 2   0\n");
2097   printf("    # Five segments with boundary markers.\n");
2098   printf("    5 1\n");
2099   printf("     1   1 2   5      # Left side of outer box.\n");
2100   printf("     2   5 7   0      # Segments 2 through 5 enclose the hole.\n");
2101   printf("     3   7 8   0\n");
2102   printf("     4   8 6   10\n");
2103   printf("     5   6 5   0\n");
2104   printf("    # One hole in the middle of the inner square.\n");
2105   printf("    1\n");
2106   printf("     1   1.5 1.5\n\n");
2107   printf(
2108 "  Note that some segments are missing from the outer square, so one must\n");
2109   printf(
2110 "  use the `-c' switch.  After `triangle -pqc box.poly', here is the output\n"
2111 );
2112   printf(
2113 "  file `box.1.node', with twelve points.  The last four points were added\n");
2114   printf(
2115 "  to meet the angle constraint.  Points 1, 2, and 9 have markers from\n");
2116   printf(
2117 "  segment 1.  Points 6 and 8 have markers from segment 4.  All the other\n");
2118   printf(
2119 "  points but 4 have been marked to indicate that they lie on a boundary.\n");
2120   printf("\n");
2121   printf("    12  2  0  1\n");
2122   printf("       1    0   0      5\n");
2123   printf("       2    0   3      5\n");
2124   printf("       3    3   0      1\n");
2125   printf("       4    3   3     33\n");
2126   printf("       5    1   1      1\n");
2127   printf("       6    1   2     10\n");
2128   printf("       7    2   1      1\n");
2129   printf("       8    2   2     10\n");
2130   printf("       9    0   1.5    5\n");
2131   printf("      10    1.5   0    1\n");
2132   printf("      11    3   1.5    1\n");
2133   printf("      12    1.5   3    1\n");
2134   printf("    # Generated by triangle -pqc box.poly\n\n");
2135   printf("  Here is the output file `box.1.ele', with twelve triangles.\n\n");
2136   printf("    12  3  0\n");
2137   printf("       1     5   6   9\n");
2138   printf("       2    10   3   7\n");
2139   printf("       3     6   8  12\n");
2140   printf("       4     9   1   5\n");
2141   printf("       5     6   2   9\n");
2142   printf("       6     7   3  11\n");
2143   printf("       7    11   4   8\n");
2144   printf("       8     7   5  10\n");
2145   printf("       9    12   2   6\n");
2146   printf("      10     8   7  11\n");
2147   printf("      11     5   1  10\n");
2148   printf("      12     8   4  12\n");
2149   printf("    # Generated by triangle -pqc box.poly\n\n");
2150   printf(
2151 "  Here is the output file `box.1.poly'.  Note that segments have been added\n"
2152 );
2153   printf(
2154 "  to represent the convex hull, and some segments have been split by newly\n"
2155 );
2156   printf(
2157 "  added points.  Note also that <# of points> is set to zero to indicate\n");
2158   printf("  that the points should be read from the .node file.\n\n");
2159   printf("    0  2  0  1\n");
2160   printf("    12  1\n");
2161   printf("       1     1   9     5\n");
2162   printf("       2     5   7     1\n");
2163   printf("       3     8   7     1\n");
2164   printf("       4     6   8    10\n");
2165   printf("       5     5   6     1\n");
2166   printf("       6     3  10     1\n");
2167   printf("       7     4  11     1\n");
2168   printf("       8     2  12     1\n");
2169   printf("       9     9   2     5\n");
2170   printf("      10    10   1     1\n");
2171   printf("      11    11   3     1\n");
2172   printf("      12    12   4     1\n");
2173   printf("    1\n");
2174   printf("       1   1.5 1.5\n");
2175   printf("    # Generated by triangle -pqc box.poly\n\n");
2176   printf("Refinement and Area Constraints:\n\n");
2177   printf(
2178 "  The -r switch causes a mesh (.node and .ele files) to be read and\n");
2179   printf(
2180 "  refined.  If the -p switch is also used, a .poly file is read and used to\n"
2181 );
2182   printf(
2183 "  specify edges that are constrained and cannot be eliminated (although\n");
2184   printf(
2185 "  they can be divided into smaller edges) by the refinement process.\n");
2186   printf("\n");
2187   printf(
2188 "  When you refine a mesh, you generally want to impose tighter quality\n");
2189   printf(
2190 "  constraints.  One way to accomplish this is to use -q with a larger\n");
2191   printf(
2192 "  angle, or -a followed by a smaller area than you used to generate the\n");
2193   printf(
2194 "  mesh you are refining.  Another way to do this is to create an .area\n");
2195   printf(
2196 "  file, which specifies a maximum area for each triangle, and use the -a\n");
2197   printf(
2198 "  switch (without a number following).  Each triangle's area constraint is\n"
2199 );
2200   printf(
2201 "  applied to that triangle.  Area constraints tend to diffuse as the mesh\n");
2202   printf(
2203 "  is refined, so if there are large variations in area constraint between\n");
2204   printf("  adjacent triangles, you may not get the results you want.\n\n");
2205   printf(
2206 "  If you are refining a mesh composed of linear (three-node) elements, the\n"
2207 );
2208   printf(
2209 "  output mesh will contain all the nodes present in the input mesh, in the\n"
2210 );
2211   printf(
2212 "  same order, with new nodes added at the end of the .node file.  However,\n"
2213 );
2214   printf(
2215 "  there is no guarantee that each output element is contained in a single\n");
2216   printf(
2217 "  input element.  Often, output elements will overlap two input elements,\n");
2218   printf(
2219 "  and input edges are not present in the output mesh.  Hence, a sequence of\n"
2220 );
2221   printf(
2222 "  refined meshes will form a hierarchy of nodes, but not a hierarchy of\n");
2223   printf(
2224 "  elements.  If you a refining a mesh of higher-order elements, the\n");
2225   printf(
2226 "  hierarchical property applies only to the nodes at the corners of an\n");
2227   printf("  element; other nodes may not be present in the refined mesh.\n\n");
2228   printf(
2229 "  It is important to understand that maximum area constraints in .poly\n");
2230   printf(
2231 "  files are handled differently from those in .area files.  A maximum area\n"
2232 );
2233   printf(
2234 "  in a .poly file applies to the whole (segment-bounded) region in which a\n"
2235 );
2236   printf(
2237 "  point falls, whereas a maximum area in an .area file applies to only one\n"
2238 );
2239   printf(
2240 "  triangle.  Area constraints in .poly files are used only when a mesh is\n");
2241   printf(
2242 "  first generated, whereas area constraints in .area files are used only to\n"
2243 );
2244   printf(
2245 "  refine an existing mesh, and are typically based on a posteriori error\n");
2246   printf(
2247 "  estimates resulting from a finite element simulation on that mesh.\n");
2248   printf("\n");
2249   printf(
2250 "  `triangle -rq25 object.1' will read object.1.node and object.1.ele, then\n"
2251 );
2252   printf(
2253 "  refine the triangulation to enforce a 25 degree minimum angle, and then\n");
2254   printf(
2255 "  write the refined triangulation to object.2.node and object.2.ele.\n");
2256   printf("\n");
2257   printf(
2258 "  `triangle -rpaa6.2 z.3' will read z.3.node, z.3.ele, z.3.poly, and\n");
2259   printf(
2260 "  z.3.area.  After reconstructing the mesh and its segments, Triangle will\n"
2261 );
2262   printf(
2263 "  refine the mesh so that no triangle has area greater than 6.2, and\n");
2264   printf(
2265 "  furthermore the triangles satisfy the maximum area constraints in\n");
2266   printf(
2267 "  z.3.area.  The output is written to z.4.node, z.4.ele, and z.4.poly.\n");
2268   printf("\n");
2269   printf(
2270 "  The sequence `triangle -qa1 x', `triangle -rqa.3 x.1', `triangle -rqa.1\n");
2271   printf(
2272 "  x.2' creates a sequence of successively finer meshes x.1, x.2, and x.3,\n");
2273   printf("  suitable for multigrid.\n\n");
2274   printf("Convex Hulls and Mesh Boundaries:\n\n");
2275   printf(
2276 "  If the input is a point set (rather than a PSLG), Triangle produces its\n");
2277   printf(
2278 "  convex hull as a by-product in the output .poly file if you use the -c\n");
2279   printf(
2280 "  switch.  There are faster algorithms for finding a two-dimensional convex\n"
2281 );
2282   printf(
2283 "  hull than triangulation, of course, but this one comes for free.  If the\n"
2284 );
2285   printf(
2286 "  input is an unconstrained mesh (you are using the -r switch but not the\n");
2287   printf(
2288 "  -p switch), Triangle produces a list of its boundary edges (including\n");
2289   printf("  hole boundaries) as a by-product if you use the -c switch.\n\n");
2290   printf("Voronoi Diagrams:\n\n");
2291   printf(
2292 "  The -v switch produces a Voronoi diagram, in files suffixed .v.node and\n");
2293   printf(
2294 "  .v.edge.  For example, `triangle -v points' will read points.node,\n");
2295   printf(
2296 "  produce its Delaunay triangulation in points.1.node and points.1.ele,\n");
2297   printf(
2298 "  and produce its Voronoi diagram in points.1.v.node and points.1.v.edge.\n");
2299   printf(
2300 "  The .v.node file contains a list of all Voronoi vertices, and the .v.edge\n"
2301 );
2302   printf(
2303 "  file contains a list of all Voronoi edges, some of which may be infinite\n"
2304 );
2305   printf(
2306 "  rays.  (The choice of filenames makes it easy to run the set of Voronoi\n");
2307   printf("  vertices through Triangle, if so desired.)\n\n");
2308   printf(
2309 "  This implementation does not use exact arithmetic to compute the Voronoi\n"
2310 );
2311   printf(
2312 "  vertices, and does not check whether neighboring vertices are identical.\n"
2313 );
2314   printf(
2315 "  Be forewarned that if the Delaunay triangulation is degenerate or\n");
2316   printf(
2317 "  near-degenerate, the Voronoi diagram may have duplicate points, crossing\n"
2318 );
2319   printf(
2320 "  edges, or infinite rays whose direction vector is zero.  Also, if you\n");
2321   printf(
2322 "  generate a constrained (as opposed to conforming) Delaunay triangulation,\n"
2323 );
2324   printf(
2325 "  or if the triangulation has holes, the corresponding Voronoi diagram is\n");
2326   printf("  likely to have crossing edges and unlikely to make sense.\n\n");
2327   printf("Mesh Topology:\n\n");
2328   printf(
2329 "  You may wish to know which triangles are adjacent to a certain Delaunay\n");
2330   printf(
2331 "  edge in an .edge file, which Voronoi regions are adjacent to a certain\n");
2332   printf(
2333 "  Voronoi edge in a .v.edge file, or which Voronoi regions are adjacent to\n"
2334 );
2335   printf(
2336 "  each other.  All of this information can be found by cross-referencing\n");
2337   printf(
2338 "  output files with the recollection that the Delaunay triangulation and\n");
2339   printf("  the Voronoi diagrams are planar duals.\n\n");
2340   printf(
2341 "  Specifically, edge i of an .edge file is the dual of Voronoi edge i of\n");
2342   printf(
2343 "  the corresponding .v.edge file, and is rotated 90 degrees counterclock-\n");
2344   printf(
2345 "  wise from the Voronoi edge.  Triangle j of an .ele file is the dual of\n");
2346   printf(
2347 "  vertex j of the corresponding .v.node file; and Voronoi region k is the\n");
2348   printf("  dual of point k of the corresponding .node file.\n\n");
2349   printf(
2350 "  Hence, to find the triangles adjacent to a Delaunay edge, look at the\n");
2351   printf(
2352 "  vertices of the corresponding Voronoi edge; their dual triangles are on\n");
2353   printf(
2354 "  the left and right of the Delaunay edge, respectively.  To find the\n");
2355   printf(
2356 "  Voronoi regions adjacent to a Voronoi edge, look at the endpoints of the\n"
2357 );
2358   printf(
2359 "  corresponding Delaunay edge; their dual regions are on the right and left\n"
2360 );
2361   printf(
2362 "  of the Voronoi edge, respectively.  To find which Voronoi regions are\n");
2363   printf("  adjacent to each other, just read the list of Delaunay edges.\n");
2364   printf("\n");
2365   printf("Statistics:\n");
2366   printf("\n");
2367   printf(
2368 "  After generating a mesh, Triangle prints a count of the number of points,\n"
2369 );
2370   printf(
2371 "  triangles, edges, boundary edges, and segments in the output mesh.  If\n");
2372   printf(
2373 "  you've forgotten the statistics for an existing mesh, the -rNEP switches\n"
2374 );
2375   printf(
2376 "  (or -rpNEP if you've got a .poly file for the existing mesh) will\n");
2377   printf("  regenerate these statistics without writing any output.\n\n");
2378   printf(
2379 "  The -V switch produces extended statistics, including a rough estimate\n");
2380   printf(
2381 "  of memory use and a histogram of triangle aspect ratios and angles in the\n"
2382 );
2383   printf("  mesh.\n\n");
2384   printf("Exact Arithmetic:\n\n");
2385   printf(
2386 "  Triangle uses adaptive exact arithmetic to perform what computational\n");
2387   printf(
2388 "  geometers call the `orientation' and `incircle' tests.  If the floating-\n"
2389 );
2390   printf(
2391 "  point arithmetic of your machine conforms to the IEEE 754 standard (as\n");
2392   printf(
2393 "  most workstations do), and does not use extended precision internal\n");
2394   printf(
2395 "  registers, then your output is guaranteed to be an absolutely true\n");
2396   printf("  Delaunay or conforming Delaunay triangulation, roundoff error\n");
2397   printf(
2398 "  notwithstanding.  The word `adaptive' implies that these arithmetic\n");
2399   printf(
2400 "  routines compute the result only to the precision necessary to guarantee\n"
2401 );
2402   printf(
2403 "  correctness, so they are usually nearly as fast as their approximate\n");
2404   printf(
2405 "  counterparts.  The exact tests can be disabled with the -X switch.  On\n");
2406   printf(
2407 "  most inputs, this switch will reduce the computation time by about eight\n"
2408 );
2409   printf(
2410 "  percent - it's not worth the risk.  There are rare difficult inputs\n");
2411   printf(
2412 "  (having many collinear and cocircular points), however, for which the\n");
2413   printf(
2414 "  difference could be a factor of two.  These are precisely the inputs most\n"
2415 );
2416   printf("  likely to cause errors if you use the -X switch.\n\n");
2417   printf(
2418 "  Unfortunately, these routines don't solve every numerical problem.  Exact\n"
2419 );
2420   printf(
2421 "  arithmetic is not used to compute the positions of points, because the\n");
2422   printf(
2423 "  bit complexity of point coordinates would grow without bound.  Hence,\n");
2424   printf(
2425 "  segment intersections aren't computed exactly; in very unusual cases,\n");
2426   printf(
2427 "  roundoff error in computing an intersection point might actually lead to\n"
2428 );
2429   printf(
2430 "  an inverted triangle and an invalid triangulation.  (This is one reason\n");
2431   printf(
2432 "  to compute your own intersection points in your .poly files.)  Similarly,\n"
2433 );
2434   printf(
2435 "  exact arithmetic is not used to compute the vertices of the Voronoi\n");
2436   printf("  diagram.\n\n");
2437   printf(
2438 "  Underflow and overflow can also cause difficulties; the exact arithmetic\n"
2439 );
2440   printf(
2441 "  routines do not ameliorate out-of-bounds exponents, which can arise\n");
2442   printf(
2443 "  during the orientation and incircle tests.  As a rule of thumb, you\n");
2444   printf(
2445 "  should ensure that your input values are within a range such that their\n");
2446   printf(
2447 "  third powers can be taken without underflow or overflow.  Underflow can\n");
2448   printf(
2449 "  silently prevent the tests from being performed exactly, while overflow\n");
2450   printf("  will typically cause a floating exception.\n\n");
2451   printf("Calling Triangle from Another Program:\n\n");
2452   printf("  Read the file triangle.h for details.\n\n");
2453   printf("Troubleshooting:\n\n");
2454   printf("  Please read this section before mailing me bugs.\n\n");
2455   printf("  `My output mesh has no triangles!'\n\n");
2456   printf(
2457 "    If you're using a PSLG, you've probably failed to specify a proper set\n"
2458 );
2459   printf(
2460 "    of bounding segments, or forgotten to use the -c switch.  Or you may\n");
2461   printf(
2462 "    have placed a hole badly.  To test these possibilities, try again with\n"
2463 );
2464   printf(
2465 "    the -c and -O switches.  Alternatively, all your input points may be\n");
2466   printf(
2467 "    collinear, in which case you can hardly expect to triangulate them.\n");
2468   printf("\n");
2469   printf("  `Triangle doesn't terminate, or just crashes.'\n");
2470   printf("\n");
2471   printf(
2472 "    Bad things can happen when triangles get so small that the distance\n");
2473   printf(
2474 "    between their vertices isn't much larger than the precision of your\n");
2475   printf(
2476 "    machine's arithmetic.  If you've compiled Triangle for single-precision\n"
2477 );
2478   printf(
2479 "    arithmetic, you might do better by recompiling it for double-precision.\n"
2480 );
2481   printf(
2482 "    Then again, you might just have to settle for more lenient constraints\n"
2483 );
2484   printf(
2485 "    on the minimum angle and the maximum area than you had planned.\n");
2486   printf("\n");
2487   printf(
2488 "    You can minimize precision problems by ensuring that the origin lies\n");
2489   printf(
2490 "    inside your point set, or even inside the densest part of your\n");
2491   printf(
2492 "    mesh.  On the other hand, if you're triangulating an object whose x\n");
2493   printf(
2494 "    coordinates all fall between 6247133 and 6247134, you're not leaving\n");
2495   printf("    much floating-point precision for Triangle to work with.\n\n");
2496   printf(
2497 "    Precision problems can occur covertly if the input PSLG contains two\n");
2498   printf(
2499 "    segments that meet (or intersect) at a very small angle, or if such an\n"
2500 );
2501   printf(
2502 "    angle is introduced by the -c switch, which may occur if a point lies\n");
2503   printf(
2504 "    ever-so-slightly inside the convex hull, and is connected by a PSLG\n");
2505   printf(
2506 "    segment to a point on the convex hull.  If you don't realize that a\n");
2507   printf(
2508 "    small angle is being formed, you might never discover why Triangle is\n");
2509   printf(
2510 "    crashing.  To check for this possibility, use the -S switch (with an\n");
2511   printf(
2512 "    appropriate limit on the number of Steiner points, found by trial-and-\n"
2513 );
2514   printf(
2515 "    error) to stop Triangle early, and view the output .poly file with\n");
2516   printf(
2517 "    Show Me (described below).  Look carefully for small angles between\n");
2518   printf(
2519 "    segments; zoom in closely, as such segments might look like a single\n");
2520   printf("    segment from a distance.\n\n");
2521   printf(
2522 "    If some of the input values are too large, Triangle may suffer a\n");
2523   printf(
2524 "    floating exception due to overflow when attempting to perform an\n");
2525   printf(
2526 "    orientation or incircle test.  (Read the section on exact arithmetic\n");
2527   printf(
2528 "    above.)  Again, I recommend compiling Triangle for double (rather\n");
2529   printf("    than single) precision arithmetic.\n\n");
2530   printf(
2531 "  `The numbering of the output points doesn't match the input points.'\n");
2532   printf("\n");
2533   printf(
2534 "    You may have eaten some of your input points with a hole, or by placing\n"
2535 );
2536   printf("    them outside the area enclosed by segments.\n\n");
2537   printf(
2538 "  `Triangle executes without incident, but when I look at the resulting\n");
2539   printf(
2540 "  mesh, it has overlapping triangles or other geometric inconsistencies.'\n");
2541   printf("\n");
2542   printf(
2543 "    If you select the -X switch, Triangle's divide-and-conquer Delaunay\n");
2544   printf(
2545 "    triangulation algorithm occasionally makes mistakes due to floating-\n");
2546   printf(
2547 "    point roundoff error.  Although these errors are rare, don't use the -X\n"
2548 );
2549   printf("    switch.  If you still have problems, please report the bug.\n");
2550   printf("\n");
2551   printf(
2552 "  Strange things can happen if you've taken liberties with your PSLG.  Do\n");
2553   printf(
2554 "  you have a point lying in the middle of a segment?  Triangle sometimes\n");
2555   printf(
2556 "  copes poorly with that sort of thing.  Do you want to lay out a collinear\n"
2557 );
2558   printf(
2559 "  row of evenly spaced, segment-connected points?  Have you simply defined\n"
2560 );
2561   printf(
2562 "  one long segment connecting the leftmost point to the rightmost point,\n");
2563   printf(
2564 "  and a bunch of points lying along it?  This method occasionally works,\n");
2565   printf(
2566 "  especially with horizontal and vertical lines, but often it doesn't, and\n"
2567 );
2568   printf(
2569 "  you'll have to connect each adjacent pair of points with a separate\n");
2570   printf("  segment.  If you don't like it, tough.\n\n");
2571   printf(
2572 "  Furthermore, if you have segments that intersect other than at their\n");
2573   printf(
2574 "  endpoints, try not to let the intersections fall extremely close to PSLG\n"
2575 );
2576   printf("  points or each other.\n\n");
2577   printf(
2578 "  If you have problems refining a triangulation not produced by Triangle:\n");
2579   printf(
2580 "  Are you sure the triangulation is geometrically valid?  Is it formatted\n");
2581   printf(
2582 "  correctly for Triangle?  Are the triangles all listed so the first three\n"
2583 );
2584   printf("  points are their corners in counterclockwise order?\n\n");
2585   printf("Show Me:\n\n");
2586   printf(
2587 "  Triangle comes with a separate program named `Show Me', whose primary\n");
2588   printf(
2589 "  purpose is to draw meshes on your screen or in PostScript.  Its secondary\n"
2590 );
2591   printf(
2592 "  purpose is to check the validity of your input files, and do so more\n");
2593   printf(
2594 "  thoroughly than Triangle does.  Show Me requires that you have the X\n");
2595   printf(
2596 "  Windows system.  If you didn't receive Show Me with Triangle, complain to\n"
2597 );
2598   printf("  whomever you obtained Triangle from, then send me mail.\n\n");
2599   printf("Triangle on the Web:\n\n");
2600   printf(
2601 "  To see an illustrated, updated version of these instructions, check out\n");
2602   printf("\n");
2603   printf("    http://www.cs.cmu.edu/~quake/triangle.html\n");
2604   printf("\n");
2605   printf("A Brief Plea:\n");
2606   printf("\n");
2607   printf(
2608 "  If you use Triangle, and especially if you use it to accomplish real\n");
2609   printf(
2610 "  work, I would like very much to hear from you.  A short letter or email\n");
2611   printf(
2612 "  (to jrs@cs.cmu.edu) describing how you use Triangle will mean a lot to\n");
2613   printf(
2614 "  me.  The more people I know are using this program, the more easily I can\n"
2615 );
2616   printf(
2617 "  justify spending time on improvements and on the three-dimensional\n");
2618   printf(
2619 "  successor to Triangle, which in turn will benefit you.  Also, I can put\n");
2620   printf(
2621 "  you on a list to receive email whenever a new version of Triangle is\n");
2622   printf("  available.\n\n");
2623   printf(
2624 "  If you use a mesh generated by Triangle in a publication, please include\n"
2625 );
2626   printf("  an acknowledgment as well.\n\n");
2627   printf("Research credit:\n\n");
2628   printf(
2629 "  Of course, I can take credit for only a fraction of the ideas that made\n");
2630   printf(
2631 "  this mesh generator possible.  Triangle owes its existence to the efforts\n"
2632 );
2633   printf(
2634 "  of many fine computational geometers and other researchers, including\n");
2635   printf(
2636 "  Marshall Bern, L. Paul Chew, Boris Delaunay, Rex A. Dwyer, David\n");
2637   printf(
2638 "  Eppstein, Steven Fortune, Leonidas J. Guibas, Donald E. Knuth, C. L.\n");
2639   printf(
2640 "  Lawson, Der-Tsai Lee, Ernst P. Mucke, Douglas M. Priest, Jim Ruppert,\n");
2641   printf(
2642 "  Isaac Saias, Bruce J. Schachter, Micha Sharir, Jorge Stolfi, Christopher\n"
2643 );
2644   printf(
2645 "  J. Van Wyk, David F. Watson, and Binhai Zhu.  See the comments at the\n");
2646   printf("  beginning of the source code for references.\n\n");
2647   exit(0);
2648 }
2649
2650 #endif /* not TRILIBRARY */
2651
2652 /*****************************************************************************/
2653 /*                                                                           */
2654 /*  internalerror()   Ask the user to send me the defective product.  Exit.  */
2655 /*                                                                           */
2656 /*****************************************************************************/
2657
2658 void internalerror()
2659 {
2660   printf("  Please report this bug to jrs@cs.cmu.edu\n");
2661   printf("  Include the message above, your input data set, and the exact\n");
2662   printf("    command line you used to run Triangle.\n");
2663   exit(1);
2664 }
2665
2666 /*****************************************************************************/
2667 /*                                                                           */
2668 /*  parsecommandline()   Read the command line, identify switches, and set   */
2669 /*                       up options and file names.                          */
2670 /*                                                                           */
2671 /*  The effects of this routine are felt entirely through global variables.  */
2672 /*                                                                           */
2673 /*****************************************************************************/
2674
2675 void parsecommandline(argc, argv)
2676 int argc;
2677 char **argv;
2678 {
2679 #ifdef TRILIBRARY
2680 #define STARTINDEX 0
2681 #else /* not TRILIBRARY */
2682 #define STARTINDEX 1
2683   int increment;
2684   int meshnumber;
2685 #endif /* not TRILIBRARY */
2686   int i, j;
2687 #ifndef CDT_ONLY
2688   int k;
2689   char workstring[FILENAMESIZE];
2690 #endif
2691
2692   poly = refine = quality = vararea = fixedarea = regionattrib = convex = 0;
2693   firstnumber = 1;
2694   edgesout = voronoi = neighbors = geomview = 0;
2695   nobound = nopolywritten = nonodewritten = noelewritten = noiterationnum = 0;
2696   noholes = noexact = 0;
2697   incremental = sweepline = 0;
2698   dwyer = 1;
2699   splitseg = 0;
2700   docheck = 0;
2701   nobisect = 0;
2702   steiner = -1;
2703   order = 1;
2704   minangle = 0.0;
2705   maxarea = -1.0;
2706   quiet = verbose = 0;
2707 #ifndef TRILIBRARY
2708   innodefilename[0] = '\0';
2709 #endif /* not TRILIBRARY */
2710
2711   for (i = STARTINDEX; i < argc; i++) {
2712 #ifndef TRILIBRARY
2713     if (argv[i][0] == '-') {
2714 #endif /* not TRILIBRARY */
2715       for (j = STARTINDEX; argv[i][j] != '\0'; j++) {
2716         if (argv[i][j] == 'p') {
2717           poly = 1;
2718         }
2719 #ifndef CDT_ONLY
2720         if (argv[i][j] == 'r') {
2721           refine = 1;
2722         }
2723         if (argv[i][j] == 'q') {
2724           quality = 1;
2725           if (((argv[i][j + 1] >= '0') && (argv[i][j + 1] <= '9')) ||
2726               (argv[i][j + 1] == '.')) {
2727             k = 0;
2728             while (((argv[i][j + 1] >= '0') && (argv[i][j + 1] <= '9')) ||
2729                    (argv[i][j + 1] == '.')) {
2730               j++;
2731               workstring[k] = argv[i][j];
2732               k++;
2733             }
2734             workstring[k] = '\0';
2735             minangle = (REAL) strtod(workstring, (char **) NULL);
2736           } else {
2737             minangle = 20.0;
2738           }
2739         }
2740         if (argv[i][j] == 'a') {
2741           quality = 1;
2742           if (((argv[i][j + 1] >= '0') && (argv[i][j + 1] <= '9')) ||
2743               (argv[i][j + 1] == '.')) {
2744             fixedarea = 1;
2745             k = 0;
2746             while (((argv[i][j + 1] >= '0') && (argv[i][j + 1] <= '9')) ||
2747                    (argv[i][j + 1] == '.')) {
2748               j++;
2749               workstring[k] = argv[i][j];
2750               k++;
2751             }
2752             workstring[k] = '\0';
2753             maxarea = (REAL) strtod(workstring, (char **) NULL);
2754             if (maxarea <= 0.0) {
2755               printf("Error:  Maximum area must be greater than zero.\n");
2756               exit(1);
2757             }
2758           } else {
2759             vararea = 1;
2760           }
2761         }
2762 #endif /* not CDT_ONLY */
2763         if (argv[i][j] == 'A') {
2764           regionattrib = 1;
2765         }
2766         if (argv[i][j] == 'c') {
2767           convex = 1;
2768         }
2769         if (argv[i][j] == 'z') {
2770           firstnumber = 0;
2771         }
2772         if (argv[i][j] == 'e') {
2773           edgesout = 1;
2774         }
2775         if (argv[i][j] == 'v') {
2776           voronoi = 1;
2777         }
2778         if (argv[i][j] == 'n') {
2779           neighbors = 1;
2780         }
2781         if (argv[i][j] == 'g') {
2782           geomview = 1;
2783         }
2784         if (argv[i][j] == 'B') {
2785           nobound = 1;
2786         }
2787         if (argv[i][j] == 'P') {
2788           nopolywritten = 1;
2789         }
2790         if (argv[i][j] == 'N') {
2791           nonodewritten = 1;
2792         }
2793         if (argv[i][j] == 'E') {
2794           noelewritten = 1;
2795         }
2796 #ifndef TRILIBRARY
2797         if (argv[i][j] == 'I') {
2798           noiterationnum = 1;
2799         }
2800 #endif /* not TRILIBRARY */
2801         if (argv[i][j] == 'O') {
2802           noholes = 1;
2803         }
2804         if (argv[i][j] == 'X') {
2805           noexact = 1;
2806         }
2807         if (argv[i][j] == 'o') {
2808           if (argv[i][j + 1] == '2') {
2809             j++;
2810             order = 2;
2811           }
2812         }
2813 #ifndef CDT_ONLY
2814         if (argv[i][j] == 'Y') {
2815           nobisect++;
2816         }
2817         if (argv[i][j] == 'S') {
2818           steiner = 0;
2819           while ((argv[i][j + 1] >= '0') && (argv[i][j + 1] <= '9')) {
2820             j++;
2821             steiner = steiner * 10 + (int) (argv[i][j] - '0');
2822           }
2823         }
2824 #endif /* not CDT_ONLY */
2825 #ifndef REDUCED
2826         if (argv[i][j] == 'i') {
2827           incremental = 1;
2828         }
2829         if (argv[i][j] == 'F') {
2830           sweepline = 1;
2831         }
2832 #endif /* not REDUCED */
2833         if (argv[i][j] == 'l') {
2834           dwyer = 0;
2835         }
2836 #ifndef REDUCED
2837 #ifndef CDT_ONLY
2838         if (argv[i][j] == 's') {
2839           splitseg = 1;
2840         }
2841 #endif /* not CDT_ONLY */
2842         if (argv[i][j] == 'C') {
2843           docheck = 1;
2844         }
2845 #endif /* not REDUCED */
2846         if (argv[i][j] == 'Q') {
2847           quiet = 1;
2848         }
2849         if (argv[i][j] == 'V') {
2850           verbose++;
2851         }
2852 #ifndef TRILIBRARY
2853         if ((argv[i][j] == 'h') || (argv[i][j] == 'H') ||
2854             (argv[i][j] == '?')) {
2855           info();
2856         }
2857 #endif /* not TRILIBRARY */
2858       }
2859 #ifndef TRILIBRARY
2860     } else {
2861       strncpy(innodefilename, argv[i], FILENAMESIZE - 1);
2862       innodefilename[FILENAMESIZE - 1] = '\0';
2863     }
2864 #endif /* not TRILIBRARY */
2865   }
2866 #ifndef TRILIBRARY
2867   if (innodefilename[0] == '\0') {
2868     syntax();
2869   }
2870   if (!strcmp(&innodefilename[strlen(innodefilename) - 5], ".node")) {
2871     innodefilename[strlen(innodefilename) - 5] = '\0';
2872   }
2873   if (!strcmp(&innodefilename[strlen(innodefilename) - 5], ".poly")) {
2874     innodefilename[strlen(innodefilename) - 5] = '\0';
2875     poly = 1;
2876   }
2877 #ifndef CDT_ONLY
2878   if (!strcmp(&innodefilename[strlen(innodefilename) - 4], ".ele")) {
2879     innodefilename[strlen(innodefilename) - 4] = '\0';
2880     refine = 1;
2881   }
2882   if (!strcmp(&innodefilename[strlen(innodefilename) - 5], ".area")) {
2883     innodefilename[strlen(innodefilename) - 5] = '\0';
2884     refine = 1;
2885     quality = 1;
2886     vararea = 1;
2887   }
2888 #endif /* not CDT_ONLY */
2889 #endif /* not TRILIBRARY */
2890   steinerleft = steiner;
2891   useshelles = poly || refine || quality || convex;
2892   goodangle = (REAL)cos(minangle * PI / 180.0);
2893   goodangle *= goodangle;
2894   if (refine && noiterationnum) {
2895     printf(
2896       "Error:  You cannot use the -I switch when refining a triangulation.\n");
2897     exit(1);
2898   }
2899   /* Be careful not to allocate space for element area constraints that */
2900   /*   will never be assigned any value (other than the default -1.0).  */
2901   if (!refine && !poly) {
2902     vararea = 0;
2903   }
2904   /* Be careful not to add an extra attribute to each element unless the */
2905   /*   input supports it (PSLG in, but not refining a preexisting mesh). */
2906   if (refine || !poly) {
2907     regionattrib = 0;
2908   }
2909
2910 #ifndef TRILIBRARY
2911   strcpy(inpolyfilename, innodefilename);
2912   strcpy(inelefilename, innodefilename);
2913   strcpy(areafilename, innodefilename);
2914   increment = 0;
2915   strcpy(workstring, innodefilename);
2916   j = 1;
2917   while (workstring[j] != '\0') {
2918     if ((workstring[j] == '.') && (workstring[j + 1] != '\0')) {
2919       increment = j + 1;
2920     }
2921     j++;
2922   }
2923   meshnumber = 0;
2924   if (increment > 0) {
2925     j = increment;
2926     do {
2927       if ((workstring[j] >= '0') && (workstring[j] <= '9')) {
2928         meshnumber = meshnumber * 10 + (int) (workstring[j] - '0');
2929       } else {
2930         increment = 0;
2931       }
2932       j++;
2933     } while (workstring[j] != '\0');
2934   }
2935   if (noiterationnum) {
2936     strcpy(outnodefilename, innodefilename);
2937     strcpy(outelefilename, innodefilename);
2938     strcpy(edgefilename, innodefilename);
2939     strcpy(vnodefilename, innodefilename);
2940     strcpy(vedgefilename, innodefilename);
2941     strcpy(neighborfilename, innodefilename);
2942     strcpy(offfilename, innodefilename);
2943     strcat(outnodefilename, ".node");
2944     strcat(outelefilename, ".ele");
2945     strcat(edgefilename, ".edge");
2946     strcat(vnodefilename, ".v.node");
2947     strcat(vedgefilename, ".v.edge");
2948     strcat(neighborfilename, ".neigh");
2949     strcat(offfilename, ".off");
2950   } else if (increment == 0) {
2951     strcpy(outnodefilename, innodefilename);
2952     strcpy(outpolyfilename, innodefilename);
2953     strcpy(outelefilename, innodefilename);
2954     strcpy(edgefilename, innodefilename);
2955     strcpy(vnodefilename, innodefilename);
2956     strcpy(vedgefilename, innodefilename);
2957     strcpy(neighborfilename, innodefilename);
2958     strcpy(offfilename, innodefilename);
2959     strcat(outnodefilename, ".1.node");
2960     strcat(outpolyfilename, ".1.poly");
2961     strcat(outelefilename, ".1.ele");
2962     strcat(edgefilename, ".1.edge");
2963     strcat(vnodefilename, ".1.v.node");
2964     strcat(vedgefilename, ".1.v.edge");
2965     strcat(neighborfilename, ".1.neigh");
2966     strcat(offfilename, ".1.off");
2967   } else {
2968     workstring[increment] = '%';
2969     workstring[increment + 1] = 'd';
2970     workstring[increment + 2] = '\0';
2971     sprintf(outnodefilename, workstring, meshnumber + 1);
2972     strcpy(outpolyfilename, outnodefilename);
2973     strcpy(outelefilename, outnodefilename);
2974     strcpy(edgefilename, outnodefilename);
2975     strcpy(vnodefilename, outnodefilename);
2976     strcpy(vedgefilename, outnodefilename);
2977     strcpy(neighborfilename, outnodefilename);
2978     strcpy(offfilename, outnodefilename);
2979     strcat(outnodefilename, ".node");
2980     strcat(outpolyfilename, ".poly");
2981     strcat(outelefilename, ".ele");
2982     strcat(edgefilename, ".edge");
2983     strcat(vnodefilename, ".v.node");
2984     strcat(vedgefilename, ".v.edge");
2985     strcat(neighborfilename, ".neigh");
2986     strcat(offfilename, ".off");
2987   }
2988   strcat(innodefilename, ".node");
2989   strcat(inpolyfilename, ".poly");
2990   strcat(inelefilename, ".ele");
2991   strcat(areafilename, ".area");
2992 #endif /* not TRILIBRARY */
2993 }
2994
2995 /**                                                                         **/
2996 /**                                                                         **/
2997 /********* User interaction routines begin here                      *********/
2998
2999 /********* Debugging routines begin here                             *********/
3000 /**                                                                         **/
3001 /**                                                                         **/
3002
3003 /*****************************************************************************/
3004 /*                                                                           */
3005 /*  printtriangle()   Print out the details of a triangle/edge handle.       */
3006 /*                                                                           */
3007 /*  I originally wrote this procedure to simplify debugging; it can be       */
3008 /*  called directly from the debugger, and presents information about a      */
3009 /*  triangle/edge handle in digestible form.  It's also used when the        */
3010 /*  highest level of verbosity (`-VVV') is specified.                        */
3011 /*                                                                           */
3012 /*****************************************************************************/
3013
3014 void printtriangle(t)
3015 struct triedge *t;
3016 {
3017   struct triedge printtri;
3018   struct edge printsh;
3019   point printpoint;
3020
3021   printf("triangle x%lx with orientation %d:\n", (unsigned long) t->tri,
3022          t->orient);
3023   decode(t->tri[0], printtri);
3024   if (printtri.tri == dummytri) {
3025     printf("    [0] = Outer space\n");
3026   } else {
3027     printf("    [0] = x%lx  %d\n", (unsigned long) printtri.tri,
3028            printtri.orient);
3029   }
3030   decode(t->tri[1], printtri);
3031   if (printtri.tri == dummytri) {
3032     printf("    [1] = Outer space\n");
3033   } else {
3034     printf("    [1] = x%lx  %d\n", (unsigned long) printtri.tri,
3035            printtri.orient);
3036   }
3037   decode(t->tri[2], printtri);
3038   if (printtri.tri == dummytri) {
3039     printf("    [2] = Outer space\n");
3040   } else {
3041     printf("    [2] = x%lx  %d\n", (unsigned long) printtri.tri,
3042            printtri.orient);
3043   }
3044   org(*t, printpoint);
3045   if (printpoint == (point) NULL)