开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 图形图像 > 3D图形编程 > 查看源码
triangle.c
资源名称：triangle.zip　[点击查看]
上传用户：yflamps
上传日期：2010-04-01
资源大小：155k
文件大小：636k
源码类别：

3D图形编程

开发平台：

C/C++

triangle.c：源码内容
  1. /*****************************************************************************/
  2. /*                                                                           */
  3. /*      888888888        ,o,                          / 888                  */
  4. /*         888    88o88o  "    o8888o  88o8888o o88888o 888  o88888o         */
  5. /*         888    888    888       88b 888  888 888 888 888 d888  88b        */
  6. /*         888    888    888  o88^o888 888  888 "88888" 888 8888oo888        */
  7. /*         888    888    888 C888  888 888  888  /      888 q888             */
  8. /*         888    888    888  "88o^888 888  888 Cb      888  "88oooo"        */
  9. /*                                              "8oo8D                       */
  10. /*                                                                           */
  11. /*  A Two-Dimensional Quality Mesh Generator and Delaunay Triangulator.      */
  12. /*  (triangle.c)                                                             */
  13. /*                                                                           */
  14. /*  Version 1.6                                                              */
  15. /*  July 28, 2005                                                            */
  16. /*                                                                           */
  17. /*  Copyright 1993, 1995, 1997, 1998, 2002, 2005                             */
  18. /*  Jonathan Richard Shewchuk                                                */
  19. /*  2360 Woolsey #H                                                          */
  20. /*  Berkeley, California  94705-1927                                         */
  21. /*  jrs@cs.berkeley.edu                                                      */
  22. /*                                                                           */
  23. /*  This program may be freely redistributed under the condition that the    */
  24. /*    copyright notices (including this entire header and the copyright      */
  25. /*    notice printed when the `-h' switch is selected) are not removed, and  */
  26. /*    no compensation is received.  Private, research, and institutional     */
  27. /*    use is free.  You may distribute modified versions of this code UNDER  */
  28. /*    THE CONDITION THAT THIS CODE AND ANY MODIFICATIONS MADE TO IT IN THE   */
  29. /*    SAME FILE REMAIN UNDER COPYRIGHT OF THE ORIGINAL AUTHOR, BOTH SOURCE   */
  30. /*    AND OBJECT CODE ARE MADE FREELY AVAILABLE WITHOUT CHARGE, AND CLEAR    */
  31. /*    NOTICE IS GIVEN OF THE MODIFICATIONS.  Distribution of this code as    */
  32. /*    part of a commercial system is permissible ONLY BY DIRECT ARRANGEMENT  */
  33. /*    WITH THE AUTHOR.  (If you are not directly supplying this code to a    */
  34. /*    customer, and you are instead telling them how they can obtain it for  */
  35. /*    free, then you are not required to make any arrangement with me.)      */
  36. /*                                                                           */
  37. /*  Hypertext instructions for Triangle are available on the Web at          */
  38. /*                                                                           */
  39. /*      http://www.cs.cmu.edu/~quake/triangle.html                           */
  40. /*                                                                           */
  41. /*  Disclaimer:  Neither I nor Carnegie Mellon warrant this code in any way  */
  42. /*    whatsoever.  This code is provided "as-is".  Use at your own risk.     */
  43. /*                                                                           */
  44. /*  Some of the references listed below are marked with an asterisk.  [*]    */
  45. /*    These references are available for downloading from the Web page       */
  46. /*                                                                           */
  47. /*      http://www.cs.cmu.edu/~quake/triangle.research.html                  */
  48. /*                                                                           */
  49. /*  Three papers discussing aspects of Triangle are available.  A short      */
  50. /*    overview appears in "Triangle:  Engineering a 2D Quality Mesh          */
  51. /*    Generator and Delaunay Triangulator," in Applied Computational         */
  52. /*    Geometry:  Towards Geometric Engineering, Ming C. Lin and Dinesh       */
  53. /*    Manocha, editors, Lecture Notes in Computer Science volume 1148,       */
  54. /*    pages 203-222, Springer-Verlag, Berlin, May 1996 (from the First ACM   */
  55. /*    Workshop on Applied Computational Geometry).  [*]                      */
  56. /*                                                                           */
  57. /*    The algorithms are discussed in the greatest detail in "Delaunay       */
  58. /*    Refinement Algorithms for Triangular Mesh Generation," Computational   */
  59. /*    Geometry:  Theory and Applications 22(1-3):21-74, May 2002.  [*]       */
  60. /*                                                                           */
  61. /*    More detail about the data structures may be found in my dissertation: */
  62. /*    "Delaunay Refinement Mesh Generation," Ph.D. thesis, Technical Report  */
  63. /*    CMU-CS-97-137, School of Computer Science, Carnegie Mellon University, */
  64. /*    Pittsburgh, Pennsylvania, 18 May 1997.  [*]                            */
  65. /*                                                                           */
  66. /*  Triangle was created as part of the Quake Project in the School of       */
  67. /*    Computer Science at Carnegie Mellon University.  For further           */
  68. /*    information, see Hesheng Bao, Jacobo Bielak, Omar Ghattas, Loukas F.   */
  69. /*    Kallivokas, David R. O'Hallaron, Jonathan R. Shewchuk, and Jifeng Xu,  */
  70. /*    "Large-scale Simulation of Elastic Wave Propagation in Heterogeneous   */
  71. /*    Media on Parallel Computers," Computer Methods in Applied Mechanics    */
  72. /*    and Engineering 152(1-2):85-102, 22 January 1998.                      */
  73. /*                                                                           */
  74. /*  Triangle's Delaunay refinement algorithm for quality mesh generation is  */
  75. /*    a hybrid of one due to Jim Ruppert, "A Delaunay Refinement Algorithm   */
  76. /*    for Quality 2-Dimensional Mesh Generation," Journal of Algorithms      */
  77. /*    18(3):548-585, May 1995 [*], and one due to L. Paul Chew, "Guaranteed- */
  78. /*    Quality Mesh Generation for Curved Surfaces," Proceedings of the Ninth */
  79. /*    Annual Symposium on Computational Geometry (San Diego, California),    */
  80. /*    pages 274-280, Association for Computing Machinery, May 1993,          */
  81. /*    http://portal.acm.org/citation.cfm?id=161150 .                         */
  82. /*                                                                           */
  83. /*  The Delaunay refinement algorithm has been modified so that it meshes    */
  84. /*    domains with small input angles well, as described in Gary L. Miller,  */
  85. /*    Steven E. Pav, and Noel J. Walkington, "When and Why Ruppert's         */
  86. /*    Algorithm Works," Twelfth International Meshing Roundtable, pages      */
  87. /*    91-102, Sandia National Laboratories, September 2003.  [*]             */
  88. /*                                                                           */
  89. /*  My implementation of the divide-and-conquer and incremental Delaunay     */
  90. /*    triangulation algorithms follows closely the presentation of Guibas    */
  91. /*    and Stolfi, even though I use a triangle-based data structure instead  */
  92. /*    of their quad-edge data structure.  (In fact, I originally implemented */
  93. /*    Triangle using the quad-edge data structure, but the switch to a       */
  94. /*    triangle-based data structure sped Triangle by a factor of two.)  The  */
  95. /*    mesh manipulation primitives and the two aforementioned Delaunay       */
  96. /*    triangulation algorithms are described by Leonidas J. Guibas and Jorge */
  97. /*    Stolfi, "Primitives for the Manipulation of General Subdivisions and   */
  98. /*    the Computation of Voronoi Diagrams," ACM Transactions on Graphics     */
  99. /*    4(2):74-123, April 1985, http://portal.acm.org/citation.cfm?id=282923 .*/
  100. /*                                                                           */
  101. /*  Their O(n log n) divide-and-conquer algorithm is adapted from Der-Tsai   */
  102. /*    Lee and Bruce J. Schachter, "Two Algorithms for Constructing the       */
  103. /*    Delaunay Triangulation," International Journal of Computer and         */
  104. /*    Information Science 9(3):219-242, 1980.  Triangle's improvement of the */
  105. /*    divide-and-conquer algorithm by alternating between vertical and       */
  106. /*    horizontal cuts was introduced by Rex A. Dwyer, "A Faster Divide-and-  */
  107. /*    Conquer Algorithm for Constructing Delaunay Triangulations,"           */
  108. /*    Algorithmica 2(2):137-151, 1987.                                       */
  109. /*                                                                           */
  110. /*  The incremental insertion algorithm was first proposed by C. L. Lawson,  */
  111. /*    "Software for C1 Surface Interpolation," in Mathematical Software III, */
  112. /*    John R. Rice, editor, Academic Press, New York, pp. 161-194, 1977.     */
  113. /*    For point location, I use the algorithm of Ernst P. Mucke, Isaac       */
  114. /*    Saias, and Binhai Zhu, "Fast Randomized Point Location Without         */
  115. /*    Preprocessing in Two- and Three-Dimensional Delaunay Triangulations,"  */
  116. /*    Proceedings of the Twelfth Annual Symposium on Computational Geometry, */
  117. /*    ACM, May 1996.  [*]  If I were to randomize the order of vertex        */
  118. /*    insertion (I currently don't bother), their result combined with the   */
  119. /*    result of Kenneth L. Clarkson and Peter W. Shor, "Applications of      */
  120. /*    Random Sampling in Computational Geometry II," Discrete &              */
  121. /*    Computational Geometry 4(1):387-421, 1989, would yield an expected     */
  122. /*    O(n^{4/3}) bound on running time.                                      */
  123. /*                                                                           */
  124. /*  The O(n log n) sweepline Delaunay triangulation algorithm is taken from  */
  125. /*    Steven Fortune, "A Sweepline Algorithm for Voronoi Diagrams",          */
  126. /*    Algorithmica 2(2):153-174, 1987.  A random sample of edges on the      */
  127. /*    boundary of the triangulation are maintained in a splay tree for the   */
  128. /*    purpose of point location.  Splay trees are described by Daniel        */
  129. /*    Dominic Sleator and Robert Endre Tarjan, "Self-Adjusting Binary Search */
  130. /*    Trees," Journal of the ACM 32(3):652-686, July 1985,                   */
  131. /*    http://portal.acm.org/citation.cfm?id=3835 .                           */
  132. /*                                                                           */
  133. /*  The algorithms for exact computation of the signs of determinants are    */
  134. /*    described in Jonathan Richard Shewchuk, "Adaptive Precision Floating-  */
  135. /*    Point Arithmetic and Fast Robust Geometric Predicates," Discrete &     */
  136. /*    Computational Geometry 18(3):305-363, October 1997.  (Also available   */
  137. /*    as Technical Report CMU-CS-96-140, School of Computer Science,         */
  138. /*    Carnegie Mellon University, Pittsburgh, Pennsylvania, May 1996.)  [*]  */
  139. /*    An abbreviated version appears as Jonathan Richard Shewchuk, "Robust   */
  140. /*    Adaptive Floating-Point Geometric Predicates," Proceedings of the      */
  141. /*    Twelfth Annual Symposium on Computational Geometry, ACM, May 1996. [*] */
  142. /*    Many of the ideas for my exact arithmetic routines originate with      */
  143. /*    Douglas M. Priest, "Algorithms for Arbitrary Precision Floating Point  */
  144. /*    Arithmetic," Tenth Symposium on Computer Arithmetic, pp. 132-143, IEEE */
  145. /*    Computer Society Press, 1991.  [*]  Many of the ideas for the correct  */
  146. /*    evaluation of the signs of determinants are taken from Steven Fortune  */
  147. /*    and Christopher J. Van Wyk, "Efficient Exact Arithmetic for Computa-   */
  148. /*    tional Geometry," Proceedings of the Ninth Annual Symposium on         */
  149. /*    Computational Geometry, ACM, pp. 163-172, May 1993, and from Steven    */
  150. /*    Fortune, "Numerical Stability of Algorithms for 2D Delaunay Triangu-   */
  151. /*    lations," International Journal of Computational Geometry & Applica-   */
  152. /*    tions 5(1-2):193-213, March-June 1995.                                 */
  153. /*                                                                           */
  154. /*  The method of inserting new vertices off-center (not precisely at the    */
  155. /*    circumcenter of every poor-quality triangle) is from Alper Ungor,      */
  156. /*    "Off-centers:  A New Type of Steiner Points for Computing Size-Optimal */
  157. /*    Quality-Guaranteed Delaunay Triangulations," Proceedings of LATIN      */
  158. /*    2004 (Buenos Aires, Argentina), April 2004.                            */
  159. /*                                                                           */
  160. /*  For definitions of and results involving Delaunay triangulations,        */
  161. /*    constrained and conforming versions thereof, and other aspects of      */
  162. /*    triangular mesh generation, see the excellent survey by Marshall Bern  */
  163. /*    and David Eppstein, "Mesh Generation and Optimal Triangulation," in    */
  164. /*    Computing and Euclidean Geometry, Ding-Zhu Du and Frank Hwang,         */
  165. /*    editors, World Scientific, Singapore, pp. 23-90, 1992.  [*]            */
  166. /*                                                                           */
  167. /*  The time for incrementally adding PSLG (planar straight line graph)      */
  168. /*    segments to create a constrained Delaunay triangulation is probably    */
  169. /*    O(t^2) per segment in the worst case and O(t) per segment in the       */
  170. /*    common case, where t is the number of triangles that intersect the     */
  171. /*    segment before it is inserted.  This doesn't count point location,     */
  172. /*    which can be much more expensive.  I could improve this to O(d log d)  */
  173. /*    time, but d is usually quite small, so it's not worth the bother.      */
  174. /*    (This note does not apply when the -s switch is used, invoking a       */
  175. /*    different method is used to insert segments.)                          */
  176. /*                                                                           */
  177. /*  The time for deleting a vertex from a Delaunay triangulation is O(d^2)   */
  178. /*    in the worst case and O(d) in the common case, where d is the degree   */
  179. /*    of the vertex being deleted.  I could improve this to O(d log d) time, */
  180. /*    but d is usually quite small, so it's not worth the bother.            */
  181. /*                                                                           */
  182. /*  Ruppert's Delaunay refinement algorithm typically generates triangles    */
  183. /*    at a linear rate (constant time per triangle) after the initial        */
  184. /*    triangulation is formed.  There may be pathological cases where        */
  185. /*    quadratic time is required, but these never arise in practice.         */
  186. /*                                                                           */
  187. /*  The geometric predicates (circumcenter calculations, segment             */
  188. /*    intersection formulae, etc.) appear in my "Lecture Notes on Geometric  */
  189. /*    Robustness" at http://www.cs.berkeley.edu/~jrs/mesh .                  */
  190. /*                                                                           */
  191. /*  If you make any improvements to this code, please please please let me   */
  192. /*    know, so that I may obtain the improvements.  Even if you don't change */
  193. /*    the code, I'd still love to hear what it's being used for.             */
  194. /*                                                                           */
  195. /*****************************************************************************/
  197. /* For single precision (which will save some memory and reduce paging),     */
  198. /*   define the symbol SINGLE by using the -DSINGLE compiler switch or by    */
  199. /*   writing "#define SINGLE" below.                                         */
  200. /*                                                                           */
  201. /* For double precision (which will allow you to refine meshes to a smaller  */
  202. /*   edge length), leave SINGLE undefined.                                   */
  203. /*                                                                           */
  204. /* Double precision uses more memory, but improves the resolution of the     */
  205. /*   meshes you can generate with Triangle.  It also reduces the likelihood  */
  206. /*   of a floating exception due to overflow.  Finally, it is much faster    */
  207. /*   than single precision on 64-bit architectures like the DEC Alpha.  I    */
  208. /*   recommend double precision unless you want to generate a mesh for which */
  209. /*   you do not have enough memory.                                          */
  211. /* #define SINGLE */
  213. #ifdef SINGLE
  214. #define REAL float
  215. #else /* not SINGLE */
  216. #define REAL double
  217. #endif /* not SINGLE */
  219. /* If yours is not a Unix system, define the NO_TIMER compiler switch to     */
  220. /*   remove the Unix-specific timing code.                                   */
  222. /* #define NO_TIMER */
  224. /* To insert lots of self-checks for internal errors, define the SELF_CHECK  */
  225. /*   symbol.  This will slow down the program significantly.  It is best to  */
  226. /*   define the symbol using the -DSELF_CHECK compiler switch, but you could */
  227. /*   write "#define SELF_CHECK" below.  If you are modifying this code, I    */
  228. /*   recommend you turn self-checks on until your work is debugged.          */
  230. /* #define SELF_CHECK */
  232. /* To compile Triangle as a callable object library (triangle.o), define the */
  233. /*   TRILIBRARY symbol.  Read the file triangle.h for details on how to call */
  234. /*   the procedure triangulate() that results.                               */
  236. /* #define TRILIBRARY */
  238. /* It is possible to generate a smaller version of Triangle using one or     */
  239. /*   both of the following symbols.  Define the REDUCED symbol to eliminate  */
  240. /*   all features that are primarily of research interest; specifically, the */
  241. /*   -i, -F, -s, and -C switches.  Define the CDT_ONLY symbol to eliminate   */
  242. /*   all meshing algorithms above and beyond constrained Delaunay            */
  243. /*   triangulation; specifically, the -r, -q, -a, -u, -D, -S, and -s         */
  244. /*   switches.  These reductions are most likely to be useful when           */
  245. /*   generating an object library (triangle.o) by defining the TRILIBRARY    */
  246. /*   symbol.                                                                 */
  248. /* #define REDUCED */
  249. /* #define CDT_ONLY */
  251. /* On some machines, my exact arithmetic routines might be defeated by the   */
  252. /*   use of internal extended precision floating-point registers.  The best  */
  253. /*   way to solve this problem is to set the floating-point registers to use */
  254. /*   single or double precision internally.  On 80x86 processors, this may   */
  255. /*   be accomplished by setting the CPU86 symbol for the Microsoft C         */
  256. /*   compiler, or the LINUX symbol for the gcc compiler running on Linux.    */
  257. /*                                                                           */
  258. /* An inferior solution is to declare certain values as `volatile', thus     */
  259. /*   forcing them to be stored to memory and rounded off.  Unfortunately,    */
  260. /*   this solution might slow Triangle down quite a bit.  To use volatile    */
  261. /*   values, write "#define INEXACT volatile" below.  Normally, however,     */
  262. /*   INEXACT should be defined to be nothing.  ("#define INEXACT".)          */
  263. /*                                                                           */
  264. /* For more discussion, see http://www.cs.cmu.edu/~quake/robust.pc.html .    */
  265. /*   For yet more discussion, see Section 5 of my paper, "Adaptive Precision */
  266. /*   Floating-Point Arithmetic and Fast Robust Geometric Predicates" (also   */
  267. /*   available as Section 6.6 of my dissertation).                           */
  269. /* #define CPU86 */
  270. /* #define LINUX */
  272. #define INEXACT /* Nothing */
  273. /* #define INEXACT volatile */
  275. /* Maximum number of characters in a file name (including the null).         */
  277. #define FILENAMESIZE 2048
  279. /* Maximum number of characters in a line read from a file (including the    */
  280. /*   null).                                                                  */
  282. #define INPUTLINESIZE 1024
  284. /* For efficiency, a variety of data structures are allocated in bulk.  The  */
  285. /*   following constants determine how many of each structure is allocated   */
  286. /*   at once.                                                                */
  288. #define TRIPERBLOCK 4092           /* Number of triangles allocated at once. */
  289. #define SUBSEGPERBLOCK 508       /* Number of subsegments allocated at once. */
  290. #define VERTEXPERBLOCK 4092         /* Number of vertices allocated at once. */
  291. #define VIRUSPERBLOCK 1020   /* Number of virus triangles allocated at once. */
  292. /* Number of encroached subsegments allocated at once. */
  293. #define BADSUBSEGPERBLOCK 252
  294. /* Number of skinny triangles allocated at once. */
  295. #define BADTRIPERBLOCK 4092
  296. /* Number of flipped triangles allocated at once. */
  297. #define FLIPSTACKERPERBLOCK 252
  298. /* Number of splay tree nodes allocated at once. */
  299. #define SPLAYNODEPERBLOCK 508
  301. /* The vertex types.   A DEADVERTEX has been deleted entirely.  An           */
  302. /*   UNDEADVERTEX is not part of the mesh, but is written to the output      */
  303. /*   .node file and affects the node indexing in the other output files.     */
  305. #define INPUTVERTEX 0
  306. #define SEGMENTVERTEX 1
  307. #define FREEVERTEX 2
  308. #define DEADVERTEX -32768
  309. #define UNDEADVERTEX -32767
  311. /* The next line is used to outsmart some very stupid compilers.  If your    */
  312. /*   compiler is smarter, feel free to replace the "int" with "void".        */
  313. /*   Not that it matters.                                                    */
  315. #define VOID int
  317. /* Two constants for algorithms based on random sampling.  Both constants    */
  318. /*   have been chosen empirically to optimize their respective algorithms.   */
  320. /* Used for the point location scheme of Mucke, Saias, and Zhu, to decide    */
  321. /*   how large a random sample of triangles to inspect.                      */
  323. #define SAMPLEFACTOR 11
  325. /* Used in Fortune's sweepline Delaunay algorithm to determine what fraction */
  326. /*   of boundary edges should be maintained in the splay tree for point      */
  327. /*   location on the front.                                                  */
  329. #define SAMPLERATE 10
  331. /* A number that speaks for itself, every kissable digit.                    */
  333. #define PI 3.141592653589793238462643383279502884197169399375105820974944592308
  335. /* Another fave.                                                             */
  337. #define SQUAREROOTTWO 1.4142135623730950488016887242096980785696718753769480732
  339. /* And here's one for those of you who are intimidated by math.              */
  341. #define ONETHIRD 0.333333333333333333333333333333333333333333333333333333333333
  343. #include <stdio.h>
  344. #include <stdlib.h>
  345. #include <string.h>
  346. #include <math.h>
  347. #ifndef NO_TIMER
  348. #include <sys/time.h>
  349. #endif /* not NO_TIMER */
  350. #ifdef CPU86
  351. #include <float.h>
  352. #endif /* CPU86 */
  353. #ifdef LINUX
  354. #include <fpu_control.h>
  355. #endif /* LINUX */
  356. #ifdef TRILIBRARY
  357. #include "triangle.h"
  358. #endif /* TRILIBRARY */
  360. /* A few forward declarations.                                               */
  362. #ifndef TRILIBRARY
  363. char *readline();
  364. char *findfield();
  365. #endif /* not TRILIBRARY */
  367. /* Labels that signify the result of point location.  The result of a        */
  368. /*   search indicates that the point falls in the interior of a triangle, on */
  369. /*   an edge, on a vertex, or outside the mesh.                              */
  371. enum locateresult {INTRIANGLE, ONEDGE, ONVERTEX, OUTSIDE};
  373. /* Labels that signify the result of vertex insertion.  The result indicates */
  374. /*   that the vertex was inserted with complete success, was inserted but    */
  375. /*   encroaches upon a subsegment, was not inserted because it lies on a     */
  376. /*   segment, or was not inserted because another vertex occupies the same   */
  377. /*   location.                                                               */
  379. enum insertvertexresult {SUCCESSFULVERTEX, ENCROACHINGVERTEX, VIOLATINGVERTEX,
  380.                          DUPLICATEVERTEX};
  382. /* Labels that signify the result of direction finding.  The result          */
  383. /*   indicates that a segment connecting the two query points falls within   */
  384. /*   the direction triangle, along the left edge of the direction triangle,  */
  385. /*   or along the right edge of the direction triangle.                      */
  387. enum finddirectionresult {WITHIN, LEFTCOLLINEAR, RIGHTCOLLINEAR};
  389. /*****************************************************************************/
  390. /*                                                                           */
  391. /*  The basic mesh data structures                                           */
  392. /*                                                                           */
  393. /*  There are three:  vertices, triangles, and subsegments (abbreviated      */
  394. /*  `subseg').  These three data structures, linked by pointers, comprise    */
  395. /*  the mesh.  A vertex simply represents a mesh vertex and its properties.  */
  396. /*  A triangle is a triangle.  A subsegment is a special data structure used */
  397. /*  to represent an impenetrable edge of the mesh (perhaps on the outer      */
  398. /*  boundary, on the boundary of a hole, or part of an internal boundary     */
  399. /*  separating two triangulated regions).  Subsegments represent boundaries, */
  400. /*  defined by the user, that triangles may not lie across.                  */
  401. /*                                                                           */
  402. /*  A triangle consists of a list of three vertices, a list of three         */
  403. /*  adjoining triangles, a list of three adjoining subsegments (when         */
  404. /*  segments exist), an arbitrary number of optional user-defined            */
  405. /*  floating-point attributes, and an optional area constraint.  The latter  */
  406. /*  is an upper bound on the permissible area of each triangle in a region,  */
  407. /*  used for mesh refinement.                                                */
  408. /*                                                                           */
  409. /*  For a triangle on a boundary of the mesh, some or all of the neighboring */
  410. /*  triangles may not be present.  For a triangle in the interior of the     */
  411. /*  mesh, often no neighboring subsegments are present.  Such absent         */
  412. /*  triangles and subsegments are never represented by NULL pointers; they   */
  413. /*  are represented by two special records:  `dummytri', the triangle that   */
  414. /*  fills "outer space", and `dummysub', the omnipresent subsegment.         */
  415. /*  `dummytri' and `dummysub' are used for several reasons; for instance,    */
  416. /*  they can be dereferenced and their contents examined without violating   */
  417. /*  protected memory.                                                        */
  418. /*                                                                           */
  419. /*  However, it is important to understand that a triangle includes other    */
  420. /*  information as well.  The pointers to adjoining vertices, triangles, and */
  421. /*  subsegments are ordered in a way that indicates their geometric relation */
  422. /*  to each other.  Furthermore, each of these pointers contains orientation */
  423. /*  information.  Each pointer to an adjoining triangle indicates which face */
  424. /*  of that triangle is contacted.  Similarly, each pointer to an adjoining  */
  425. /*  subsegment indicates which side of that subsegment is contacted, and how */
  426. /*  the subsegment is oriented relative to the triangle.                     */
  427. /*                                                                           */
  428. /*  The data structure representing a subsegment may be thought to be        */
  429. /*  abutting the edge of one or two triangle data structures:  either        */
  430. /*  sandwiched between two triangles, or resting against one triangle on an  */
  431. /*  exterior boundary or hole boundary.                                      */
  432. /*                                                                           */
  433. /*  A subsegment consists of a list of four vertices--the vertices of the    */
  434. /*  subsegment, and the vertices of the segment it is a part of--a list of   */
  435. /*  two adjoining subsegments, and a list of two adjoining triangles.  One   */
  436. /*  of the two adjoining triangles may not be present (though there should   */
  437. /*  always be one), and neighboring subsegments might not be present.        */
  438. /*  Subsegments also store a user-defined integer "boundary marker".         */
  439. /*  Typically, this integer is used to indicate what boundary conditions are */
  440. /*  to be applied at that location in a finite element simulation.           */
  441. /*                                                                           */
  442. /*  Like triangles, subsegments maintain information about the relative      */
  443. /*  orientation of neighboring objects.                                      */
  444. /*                                                                           */
  445. /*  Vertices are relatively simple.  A vertex is a list of floating-point    */
  446. /*  numbers, starting with the x, and y coordinates, followed by an          */
  447. /*  arbitrary number of optional user-defined floating-point attributes,     */
  448. /*  followed by an integer boundary marker.  During the segment insertion    */
  449. /*  phase, there is also a pointer from each vertex to a triangle that may   */
  450. /*  contain it.  Each pointer is not always correct, but when one is, it     */
  451. /*  speeds up segment insertion.  These pointers are assigned values once    */
  452. /*  at the beginning of the segment insertion phase, and are not used or     */
  453. /*  updated except during this phase.  Edge flipping during segment          */
  454. /*  insertion will render some of them incorrect.  Hence, don't rely upon    */
  455. /*  them for anything.                                                       */
  456. /*                                                                           */
  457. /*  Other than the exception mentioned above, vertices have no information   */
  458. /*  about what triangles, subfacets, or subsegments they are linked to.      */
  459. /*                                                                           */
  460. /*****************************************************************************/
  462. /*****************************************************************************/
  463. /*                                                                           */
  464. /*  Handles                                                                  */
  465. /*                                                                           */
  466. /*  The oriented triangle (`otri') and oriented subsegment (`osub') data     */
  467. /*  structures defined below do not themselves store any part of the mesh.   */
  468. /*  The mesh itself is made of `triangle's, `subseg's, and `vertex's.        */
  469. /*                                                                           */
  470. /*  Oriented triangles and oriented subsegments will usually be referred to  */
  471. /*  as "handles."  A handle is essentially a pointer into the mesh; it       */
  472. /*  allows you to "hold" one particular part of the mesh.  Handles are used  */
  473. /*  to specify the regions in which one is traversing and modifying the mesh.*/
  474. /*  A single `triangle' may be held by many handles, or none at all.  (The   */
  475. /*  latter case is not a memory leak, because the triangle is still          */
  476. /*  connected to other triangles in the mesh.)                               */
  477. /*                                                                           */
  478. /*  An `otri' is a handle that holds a triangle.  It holds a specific edge   */
  479. /*  of the triangle.  An `osub' is a handle that holds a subsegment.  It     */
  480. /*  holds either the left or right side of the subsegment.                   */
  481. /*                                                                           */
  482. /*  Navigation about the mesh is accomplished through a set of mesh          */
  483. /*  manipulation primitives, further below.  Many of these primitives take   */
  484. /*  a handle and produce a new handle that holds the mesh near the first     */
  485. /*  handle.  Other primitives take two handles and glue the corresponding    */
  486. /*  parts of the mesh together.  The orientation of the handles is           */
  487. /*  important.  For instance, when two triangles are glued together by the   */
  488. /*  bond() primitive, they are glued at the edges on which the handles lie.  */
  489. /*                                                                           */
  490. /*  Because vertices have no information about which triangles they are      */
  491. /*  attached to, I commonly represent a vertex by use of a handle whose      */
  492. /*  origin is the vertex.  A single handle can simultaneously represent a    */
  493. /*  triangle, an edge, and a vertex.                                         */
  494. /*                                                                           */
  495. /*****************************************************************************/
  497. /* The triangle data structure.  Each triangle contains three pointers to    */
  498. /*   adjoining triangles, plus three pointers to vertices, plus three        */
  499. /*   pointers to subsegments (declared below; these pointers are usually     */
  500. /*   `dummysub').  It may or may not also contain user-defined attributes    */
  501. /*   and/or a floating-point "area constraint."  It may also contain extra   */
  502. /*   pointers for nodes, when the user asks for high-order elements.         */
  503. /*   Because the size and structure of a `triangle' is not decided until     */
  504. /*   runtime, I haven't simply declared the type `triangle' as a struct.     */
  506. typedef REAL **triangle;            /* Really:  typedef triangle *triangle   */
  508. /* An oriented triangle:  includes a pointer to a triangle and orientation.  */
  509. /*   The orientation denotes an edge of the triangle.  Hence, there are      */
  510. /*   three possible orientations.  By convention, each edge always points    */
  511. /*   counterclockwise about the corresponding triangle.                      */
  513. struct otri {
  514.   triangle *tri;
  515.   int orient;                                         /* Ranges from 0 to 2. */
  516. };
  518. /* The subsegment data structure.  Each subsegment contains two pointers to  */
  519. /*   adjoining subsegments, plus four pointers to vertices, plus two         */
  520. /*   pointers to adjoining triangles, plus one boundary marker, plus one     */
  521. /*   segment number.                                                         */
  523. typedef REAL **subseg;                  /* Really:  typedef subseg *subseg   */
  525. /* An oriented subsegment:  includes a pointer to a subsegment and an        */
  526. /*   orientation.  The orientation denotes a side of the edge.  Hence, there */
  527. /*   are two possible orientations.  By convention, the edge is always       */
  528. /*   directed so that the "side" denoted is the right side of the edge.      */
  530. struct osub {
  531.   subseg *ss;
  532.   int ssorient;                                       /* Ranges from 0 to 1. */
  533. };
  535. /* The vertex data structure.  Each vertex is actually an array of REALs.    */
  536. /*   The number of REALs is unknown until runtime.  An integer boundary      */
  537. /*   marker, and sometimes a pointer to a triangle, is appended after the    */
  538. /*   REALs.                                                                  */
  540. typedef REAL *vertex;
  542. /* A queue used to store encroached subsegments.  Each subsegment's vertices */
  543. /*   are stored so that we can check whether a subsegment is still the same. */
  545. struct badsubseg {
  546.   subseg encsubseg;                             /* An encroached subsegment. */
  547.   vertex subsegorg, subsegdest;                         /* Its two vertices. */
  548. };
  550. /* A queue used to store bad triangles.  The key is the square of the cosine */
  551. /*   of the smallest angle of the triangle.  Each triangle's vertices are    */
  552. /*   stored so that one can check whether a triangle is still the same.      */
  554. struct badtriang {
  555.   triangle poortri;                       /* A skinny or too-large triangle. */
  556.   REAL key;                             /* cos^2 of smallest (apical) angle. */
  557.   vertex triangorg, triangdest, triangapex;           /* Its three vertices. */
  558.   struct badtriang *nexttriang;             /* Pointer to next bad triangle. */
  559. };
  561. /* A stack of triangles flipped during the most recent vertex insertion.     */
  562. /*   The stack is used to undo the vertex insertion if the vertex encroaches */
  563. /*   upon a subsegment.                                                      */
  565. struct flipstacker {
  566.   triangle flippedtri;                       /* A recently flipped triangle. */
  567.   struct flipstacker *prevflip;               /* Previous flip in the stack. */
  568. };
  570. /* A node in a heap used to store events for the sweepline Delaunay          */
  571. /*   algorithm.  Nodes do not point directly to their parents or children in */
  572. /*   the heap.  Instead, each node knows its position in the heap, and can   */
  573. /*   look up its parent and children in a separate array.  The `eventptr'    */
  574. /*   points either to a `vertex' or to a triangle (in encoded format, so     */
  575. /*   that an orientation is included).  In the latter case, the origin of    */
  576. /*   the oriented triangle is the apex of a "circle event" of the sweepline  */
  577. /*   algorithm.  To distinguish site events from circle events, all circle   */
  578. /*   events are given an invalid (smaller than `xmin') x-coordinate `xkey'.  */
  580. struct event {
  581.   REAL xkey, ykey;                              /* Coordinates of the event. */
  582.   VOID *eventptr;      /* Can be a vertex or the location of a circle event. */
  583.   int heapposition;              /* Marks this event's position in the heap. */
  584. };
  586. /* A node in the splay tree.  Each node holds an oriented ghost triangle     */
  587. /*   that represents a boundary edge of the growing triangulation.  When a   */
  588. /*   circle event covers two boundary edges with a triangle, so that they    */
  589. /*   are no longer boundary edges, those edges are not immediately deleted   */
  590. /*   from the tree; rather, they are lazily deleted when they are next       */
  591. /*   encountered.  (Since only a random sample of boundary edges are kept    */
  592. /*   in the tree, lazy deletion is faster.)  `keydest' is used to verify     */
  593. /*   that a triangle is still the same as when it entered the splay tree; if */
  594. /*   it has been rotated (due to a circle event), it no longer represents a  */
  595. /*   boundary edge and should be deleted.                                    */
  597. struct splaynode {
  598.   struct otri keyedge;                     /* Lprev of an edge on the front. */
  599.   vertex keydest;           /* Used to verify that splay node is still live. */
  600.   struct splaynode *lchild, *rchild;              /* Children in splay tree. */
  601. };
  603. /* A type used to allocate memory.  firstblock is the first block of items.  */
  604. /*   nowblock is the block from which items are currently being allocated.   */
  605. /*   nextitem points to the next slab of free memory for an item.            */
  606. /*   deaditemstack is the head of a linked list (stack) of deallocated items */
  607. /*   that can be recycled.  unallocateditems is the number of items that     */
  608. /*   remain to be allocated from nowblock.                                   */
  609. /*                                                                           */
  610. /* Traversal is the process of walking through the entire list of items, and */
  611. /*   is separate from allocation.  Note that a traversal will visit items on */
  612. /*   the "deaditemstack" stack as well as live items.  pathblock points to   */
  613. /*   the block currently being traversed.  pathitem points to the next item  */
  614. /*   to be traversed.  pathitemsleft is the number of items that remain to   */
  615. /*   be traversed in pathblock.                                              */
  616. /*                                                                           */
  617. /* alignbytes determines how new records should be aligned in memory.        */
  618. /*   itembytes is the length of a record in bytes (after rounding up).       */
  619. /*   itemsperblock is the number of items allocated at once in a single      */
  620. /*   block.  itemsfirstblock is the number of items in the first block,      */
  621. /*   which can vary from the others.  items is the number of currently       */
  622. /*   allocated items.  maxitems is the maximum number of items that have     */
  623. /*   been allocated at once; it is the current number of items plus the      */
  624. /*   number of records kept on deaditemstack.                                */
  626. struct memorypool {
  627.   VOID **firstblock, **nowblock;
  628.   VOID *nextitem;
  629.   VOID *deaditemstack;
  630.   VOID **pathblock;
  631.   VOID *pathitem;
  632.   int alignbytes;
  633.   int itembytes;
  634.   int itemsperblock;
  635.   int itemsfirstblock;
  636.   long items, maxitems;
  637.   int unallocateditems;
  638.   int pathitemsleft;
  639. };
  642. /* Global constants.                                                         */
  644. REAL splitter;       /* Used to split REAL factors for exact multiplication. */
  645. REAL epsilon;                             /* Floating-point machine epsilon. */
  646. REAL resulterrbound;
  647. REAL ccwerrboundA, ccwerrboundB, ccwerrboundC;
  648. REAL iccerrboundA, iccerrboundB, iccerrboundC;
  649. REAL o3derrboundA, o3derrboundB, o3derrboundC;
  651. /* Random number seed is not constant, but I've made it global anyway.       */
  653. unsigned long randomseed;                     /* Current random number seed. */
  656. /* Mesh data structure.  Triangle operates on only one mesh, but the mesh    */
  657. /*   structure is used (instead of global variables) to allow reentrancy.    */
  659. struct mesh {
  661. /* Variables used to allocate memory for triangles, subsegments, vertices,   */
  662. /*   viri (triangles being eaten), encroached segments, bad (skinny or too   */
  663. /*   large) triangles, and splay tree nodes.                                 */
  665.   struct memorypool triangles;
  666.   struct memorypool subsegs;
  667.   struct memorypool vertices;
  668.   struct memorypool viri;
  669.   struct memorypool badsubsegs;
  670.   struct memorypool badtriangles;
  671.   struct memorypool flipstackers;
  672.   struct memorypool splaynodes;
  674. /* Variables that maintain the bad triangle queues.  The queues are          */
  675. /*   ordered from 4095 (highest priority) to 0 (lowest priority).            */
  677.   struct badtriang *queuefront[4096];
  678.   struct badtriang *queuetail[4096];
  679.   int nextnonemptyq[4096];
  680.   int firstnonemptyq;
  682. /* Variable that maintains the stack of recently flipped triangles.          */
  684.   struct flipstacker *lastflip;
  686. /* Other variables. */
  688.   REAL xmin, xmax, ymin, ymax;                            /* x and y bounds. */
  689.   REAL xminextreme;      /* Nonexistent x value used as a flag in sweepline. */
  690.   int invertices;                               /* Number of input vertices. */
  691.   int inelements;                              /* Number of input triangles. */
  692.   int insegments;                               /* Number of input segments. */
  693.   int holes;                                       /* Number of input holes. */
  694.   int regions;                                   /* Number of input regions. */
  695.   int undeads;    /* Number of input vertices that don't appear in the mesh. */
  696.   long edges;                                     /* Number of output edges. */
  697.   int mesh_dim;                                /* Dimension (ought to be 2). */
  698.   int nextras;                           /* Number of attributes per vertex. */
  699.   int eextras;                         /* Number of attributes per triangle. */
  700.   long hullsize;                          /* Number of edges in convex hull. */
  701.   int steinerleft;                 /* Number of Steiner points not yet used. */
  702.   int vertexmarkindex;         /* Index to find boundary marker of a vertex. */
  703.   int vertex2triindex;     /* Index to find a triangle adjacent to a vertex. */
  704.   int highorderindex;  /* Index to find extra nodes for high-order elements. */
  705.   int elemattribindex;            /* Index to find attributes of a triangle. */
  706.   int areaboundindex;             /* Index to find area bound of a triangle. */
  707.   int checksegments;         /* Are there segments in the triangulation yet? */
  708.   int checkquality;                  /* Has quality triangulation begun yet? */
  709.   int readnodefile;                           /* Has a .node file been read? */
  710.   long samples;              /* Number of random samples for point location. */
  712.   long incirclecount;                 /* Number of incircle tests performed. */
  713.   long counterclockcount;     /* Number of counterclockwise tests performed. */
  714.   long orient3dcount;           /* Number of 3D orientation tests performed. */
  715.   long hyperbolacount;      /* Number of right-of-hyperbola tests performed. */
  716.   long circumcentercount;  /* Number of circumcenter calculations performed. */
  717.   long circletopcount;       /* Number of circle top calculations performed. */
  719. /* Triangular bounding box vertices.                                         */
  721.   vertex infvertex1, infvertex2, infvertex3;
  723. /* Pointer to the `triangle' that occupies all of "outer space."             */
  725.   triangle *dummytri;
  726.   triangle *dummytribase;    /* Keep base address so we can free() it later. */
  728. /* Pointer to the omnipresent subsegment.  Referenced by any triangle or     */
  729. /*   subsegment that isn't really connected to a subsegment at that          */
  730. /*   location.                                                               */
  732.   subseg *dummysub;
  733.   subseg *dummysubbase;      /* Keep base address so we can free() it later. */
  735. /* Pointer to a recently visited triangle.  Improves point location if       */
  736. /*   proximate vertices are inserted sequentially.                           */
  738.   struct otri recenttri;
  740. };                                                  /* End of `struct mesh'. */
  743. /* Data structure for command line switches and file names.  This structure  */
  744. /*   is used (instead of global variables) to allow reentrancy.              */
  746. struct behavior {
  748. /* Switches for the triangulator.                                            */
  749. /*   poly: -p switch.  refine: -r switch.                                    */
  750. /*   quality: -q switch.                                                     */
  751. /*     minangle: minimum angle bound, specified after -q switch.             */
  752. /*     goodangle: cosine squared of minangle.                                */
  753. /*     offconstant: constant used to place off-center Steiner points.        */
  754. /*   vararea: -a switch without number.                                      */
  755. /*   fixedarea: -a switch with number.                                       */
  756. /*     maxarea: maximum area bound, specified after -a switch.               */
  757. /*   usertest: -u switch.                                                    */
  758. /*   regionattrib: -A switch.  convex: -c switch.                            */
  759. /*   weighted: 1 for -w switch, 2 for -W switch.  jettison: -j switch        */
  760. /*   firstnumber: inverse of -z switch.  All items are numbered starting     */
  761. /*     from `firstnumber'.                                                   */
  762. /*   edgesout: -e switch.  voronoi: -v switch.                               */
  763. /*   neighbors: -n switch.  geomview: -g switch.                             */
  764. /*   nobound: -B switch.  nopolywritten: -P switch.                          */
  765. /*   nonodewritten: -N switch.  noelewritten: -E switch.                     */
  766. /*   noiterationnum: -I switch.  noholes: -O switch.                         */
  767. /*   noexact: -X switch.                                                     */
  768. /*   order: element order, specified after -o switch.                        */
  769. /*   nobisect: count of how often -Y switch is selected.                     */
  770. /*   steiner: maximum number of Steiner points, specified after -S switch.   */
  771. /*   incremental: -i switch.  sweepline: -F switch.                          */
  772. /*   dwyer: inverse of -l switch.                                            */
  773. /*   splitseg: -s switch.                                                    */
  774. /*   conformdel: -D switch.  docheck: -C switch.                             */
  775. /*   quiet: -Q switch.  verbose: count of how often -V switch is selected.   */
  776. /*   usesegments: -p, -r, -q, or -c switch; determines whether segments are  */
  777. /*     used at all.                                                          */
  778. /*                                                                           */
  779. /* Read the instructions to find out the meaning of these switches.          */
  781.   int poly, refine, quality, vararea, fixedarea, usertest;
  782.   int regionattrib, convex, weighted, jettison;
  783.   int firstnumber;
  784.   int edgesout, voronoi, neighbors, geomview;
  785.   int nobound, nopolywritten, nonodewritten, noelewritten, noiterationnum;
  786.   int noholes, noexact, conformdel;
  787.   int incremental, sweepline, dwyer;
  788.   int splitseg;
  789.   int docheck;
  790.   int quiet, verbose;
  791.   int usesegments;
  792.   int order;
  793.   int nobisect;
  794.   int steiner;
  795.   REAL minangle, goodangle, offconstant;
  796.   REAL maxarea;
  798. /* Variables for file names.                                                 */
  800. #ifndef TRILIBRARY
  801.   char innodefilename[FILENAMESIZE];
  802.   char inelefilename[FILENAMESIZE];
  803.   char inpolyfilename[FILENAMESIZE];
  804.   char areafilename[FILENAMESIZE];
  805.   char outnodefilename[FILENAMESIZE];
  806.   char outelefilename[FILENAMESIZE];
  807.   char outpolyfilename[FILENAMESIZE];
  808.   char edgefilename[FILENAMESIZE];
  809.   char vnodefilename[FILENAMESIZE];
  810.   char vedgefilename[FILENAMESIZE];
  811.   char neighborfilename[FILENAMESIZE];
  812.   char offfilename[FILENAMESIZE];
  813. #endif /* not TRILIBRARY */
  815. };                                              /* End of `struct behavior'. */
  818. /*****************************************************************************/
  819. /*                                                                           */
  820. /*  Mesh manipulation primitives.  Each triangle contains three pointers to  */
  821. /*  other triangles, with orientations.  Each pointer points not to the      */
  822. /*  first byte of a triangle, but to one of the first three bytes of a       */
  823. /*  triangle.  It is necessary to extract both the triangle itself and the   */
  824. /*  orientation.  To save memory, I keep both pieces of information in one   */
  825. /*  pointer.  To make this possible, I assume that all triangles are aligned */
  826. /*  to four-byte boundaries.  The decode() routine below decodes a pointer,  */
  827. /*  extracting an orientation (in the range 0 to 2) and a pointer to the     */
  828. /*  beginning of a triangle.  The encode() routine compresses a pointer to a */
  829. /*  triangle and an orientation into a single pointer.  My assumptions that  */
  830. /*  triangles are four-byte-aligned and that the `unsigned long' type is     */
  831. /*  long enough to hold a pointer are two of the few kludges in this program.*/
  832. /*                                                                           */
  833. /*  Subsegments are manipulated similarly.  A pointer to a subsegment        */
  834. /*  carries both an address and an orientation in the range 0 to 1.          */
  835. /*                                                                           */
  836. /*  The other primitives take an oriented triangle or oriented subsegment,   */
  837. /*  and return an oriented triangle or oriented subsegment or vertex; or     */
  838. /*  they change the connections in the data structure.                       */
  839. /*                                                                           */
  840. /*  Below, triangles and subsegments are denoted by their vertices.  The     */
  841. /*  triangle abc has origin (org) a, destination (dest) b, and apex (apex)   */
  842. /*  c.  These vertices occur in counterclockwise order about the triangle.   */
  843. /*  The handle abc may simultaneously denote vertex a, edge ab, and triangle */
  844. /*  abc.                                                                     */
  845. /*                                                                           */
  846. /*  Similarly, the subsegment ab has origin (sorg) a and destination (sdest) */
  847. /*  b.  If ab is thought to be directed upward (with b directly above a),    */
  848. /*  then the handle ab is thought to grasp the right side of ab, and may     */
  849. /*  simultaneously denote vertex a and edge ab.                              */
  850. /*                                                                           */
  851. /*  An asterisk (*) denotes a vertex whose identity is unknown.              */
  852. /*                                                                           */
  853. /*  Given this notation, a partial list of mesh manipulation primitives      */
  854. /*  follows.                                                                 */
  855. /*                                                                           */
  856. /*                                                                           */
  857. /*  For triangles:                                                           */
  858. /*                                                                           */
  859. /*  sym:  Find the abutting triangle; same edge.                             */
  860. /*  sym(abc) -> ba*                                                          */
  861. /*                                                                           */
  862. /*  lnext:  Find the next edge (counterclockwise) of a triangle.             */
  863. /*  lnext(abc) -> bca                                                        */
  864. /*                                                                           */
  865. /*  lprev:  Find the previous edge (clockwise) of a triangle.                */
  866. /*  lprev(abc) -> cab                                                        */
  867. /*                                                                           */
  868. /*  onext:  Find the next edge counterclockwise with the same origin.        */
  869. /*  onext(abc) -> ac*                                                        */
  870. /*                                                                           */
  871. /*  oprev:  Find the next edge clockwise with the same origin.               */
  872. /*  oprev(abc) -> a*b                                                        */
  873. /*                                                                           */
  874. /*  dnext:  Find the next edge counterclockwise with the same destination.   */
  875. /*  dnext(abc) -> *ba                                                        */
  876. /*                                                                           */
  877. /*  dprev:  Find the next edge clockwise with the same destination.          */
  878. /*  dprev(abc) -> cb*                                                        */
  879. /*                                                                           */
  880. /*  rnext:  Find the next edge (counterclockwise) of the adjacent triangle.  */
  881. /*  rnext(abc) -> *a*                                                        */
  882. /*                                                                           */
  883. /*  rprev:  Find the previous edge (clockwise) of the adjacent triangle.     */
  884. /*  rprev(abc) -> b**                                                        */
  885. /*                                                                           */
  886. /*  org:  Origin          dest:  Destination          apex:  Apex            */
  887. /*  org(abc) -> a         dest(abc) -> b              apex(abc) -> c         */
  888. /*                                                                           */
  889. /*  bond:  Bond two triangles together at the resepective handles.           */
  890. /*  bond(abc, bad)                                                           */
  891. /*                                                                           */
  892. /*                                                                           */
  893. /*  For subsegments:                                                         */
  894. /*                                                                           */
  895. /*  ssym:  Reverse the orientation of a subsegment.                          */
  896. /*  ssym(ab) -> ba                                                           */
  897. /*                                                                           */
  898. /*  spivot:  Find adjoining subsegment with the same origin.                 */
  899. /*  spivot(ab) -> a*                                                         */
  900. /*                                                                           */
  901. /*  snext:  Find next subsegment in sequence.                                */
  902. /*  snext(ab) -> b*                                                          */
  903. /*                                                                           */
  904. /*  sorg:  Origin                      sdest:  Destination                   */
  905. /*  sorg(ab) -> a                      sdest(ab) -> b                        */
  906. /*                                                                           */
  907. /*  sbond:  Bond two subsegments together at the respective origins.         */
  908. /*  sbond(ab, ac)                                                            */
  909. /*                                                                           */
  910. /*                                                                           */
  911. /*  For interacting tetrahedra and subfacets:                                */
  912. /*                                                                           */
  913. /*  tspivot:  Find a subsegment abutting a triangle.                         */
  914. /*  tspivot(abc) -> ba                                                       */
  915. /*                                                                           */
  916. /*  stpivot:  Find a triangle abutting a subsegment.                         */
  917. /*  stpivot(ab) -> ba*                                                       */
  918. /*                                                                           */
  919. /*  tsbond:  Bond a triangle to a subsegment.                                */
  920. /*  tsbond(abc, ba)                                                          */
  921. /*                                                                           */
  922. /*****************************************************************************/
  924. /********* Mesh manipulation primitives begin here                   *********/
  925. /**                                                                         **/
  926. /**                                                                         **/
  928. /* Fast lookup arrays to speed some of the mesh manipulation primitives.     */
  930. int plus1mod3[3] = {1, 2, 0};
  931. int minus1mod3[3] = {2, 0, 1};
  933. /********* Primitives for triangles                                  *********/
  934. /*                                                                           */
  935. /*                                                                           */
  937. /* decode() converts a pointer to an oriented triangle.  The orientation is  */
  938. /*   extracted from the two least significant bits of the pointer.           */
  940. #define decode(ptr, otri)                                                     
  941.   (otri).orient = (int) ((unsigned long) (ptr) & (unsigned long) 3l);         
  942.   (otri).tri = (triangle *)                                                   
  943.                   ((unsigned long) (ptr) ^ (unsigned long) (otri).orient)
  945. /* encode() compresses an oriented triangle into a single pointer.  It       */
  946. /*   relies on the assumption that all triangles are aligned to four-byte    */
  947. /*   boundaries, so the two least significant bits of (otri).tri are zero.   */
  949. #define encode(otri)                                                          
  950.   (triangle) ((unsigned long) (otri).tri | (unsigned long) (otri).orient)
  952. /* The following handle manipulation primitives are all described by Guibas  */
  953. /*   and Stolfi.  However, Guibas and Stolfi use an edge-based data          */
  954. /*   structure, whereas I use a triangle-based data structure.               */
  956. /* sym() finds the abutting triangle, on the same edge.  Note that the edge  */
  957. /*   direction is necessarily reversed, because the handle specified by an   */
  958. /*   oriented triangle is directed counterclockwise around the triangle.     */
  960. #define sym(otri1, otri2)                                                     
  961.   ptr = (otri1).tri[(otri1).orient];                                          
  962.   decode(ptr, otri2);
  964. #define symself(otri)                                                         
  965.   ptr = (otri).tri[(otri).orient];                                            
  966.   decode(ptr, otri);
  968. /* lnext() finds the next edge (counterclockwise) of a triangle.             */
  970. #define lnext(otri1, otri2)                                                   
  971.   (otri2).tri = (otri1).tri;                                                  
  972.   (otri2).orient = plus1mod3[(otri1).orient]
  974. #define lnextself(otri)                                                       
  975.   (otri).orient = plus1mod3[(otri).orient]
  977. /* lprev() finds the previous edge (clockwise) of a triangle.                */
  979. #define lprev(otri1, otri2)                                                   
  980.   (otri2).tri = (otri1).tri;                                                  
  981.   (otri2).orient = minus1mod3[(otri1).orient]
  983. #define lprevself(otri)                                                       
  984.   (otri).orient = minus1mod3[(otri).orient]
  986. /* onext() spins counterclockwise around a vertex; that is, it finds the     */
  987. /*   next edge with the same origin in the counterclockwise direction.  This */
  988. /*   edge is part of a different triangle.                                   */
  990. #define onext(otri1, otri2)                                                   
  991.   lprev(otri1, otri2);                                                        
  992.   symself(otri2);
  994. #define onextself(otri)                                                       
  995.   lprevself(otri);                                                            
  996.   symself(otri);
  998. /* oprev() spins clockwise around a vertex; that is, it finds the next edge  */
  999. /*   with the same origin in the clockwise direction.  This edge is part of  */
  1000. /*   a different triangle.                                                   */
  1002. #define oprev(otri1, otri2)                                                   
  1003.   sym(otri1, otri2);                                                          
  1004.   lnextself(otri2);
  1006. #define oprevself(otri)                                                       
  1007.   symself(otri);                                                              
  1008.   lnextself(otri);
  1010. /* dnext() spins counterclockwise around a vertex; that is, it finds the     */
  1011. /*   next edge with the same destination in the counterclockwise direction.  */
  1012. /*   This edge is part of a different triangle.                              */
  1014. #define dnext(otri1, otri2)                                                   
  1015.   sym(otri1, otri2);                                                          
  1016.   lprevself(otri2);
  1018. #define dnextself(otri)                                                       
  1019.   symself(otri);                                                              
  1020.   lprevself(otri);
  1022. /* dprev() spins clockwise around a vertex; that is, it finds the next edge  */
  1023. /*   with the same destination in the clockwise direction.  This edge is     */
  1024. /*   part of a different triangle.                                           */
  1026. #define dprev(otri1, otri2)                                                   
  1027.   lnext(otri1, otri2);                                                        
  1028.   symself(otri2);
  1030. #define dprevself(otri)                                                       
  1031.   lnextself(otri);                                                            
  1032.   symself(otri);
  1034. /* rnext() moves one edge counterclockwise about the adjacent triangle.      */
  1035. /*   (It's best understood by reading Guibas and Stolfi.  It involves        */
  1036. /*   changing triangles twice.)                                              */
  1038. #define rnext(otri1, otri2)                                                   
  1039.   sym(otri1, otri2);                                                          
  1040.   lnextself(otri2);                                                           
  1041.   symself(otri2);
  1043. #define rnextself(otri)                                                       
  1044.   symself(otri);                                                              
  1045.   lnextself(otri);                                                            
  1046.   symself(otri);
  1048. /* rprev() moves one edge clockwise about the adjacent triangle.             */
  1049. /*   (It's best understood by reading Guibas and Stolfi.  It involves        */
  1050. /*   changing triangles twice.)                                              */
  1052. #define rprev(otri1, otri2)                                                   
  1053.   sym(otri1, otri2);                                                          
  1054.   lprevself(otri2);                                                           
  1055.   symself(otri2);
  1057. #define rprevself(otri)                                                       
  1058.   symself(otri);                                                              
  1059.   lprevself(otri);                                                            
  1060.   symself(otri);
  1062. /* These primitives determine or set the origin, destination, or apex of a   */
  1063. /* triangle.                                                                 */
  1065. #define org(otri, vertexptr)                                                  
  1066.   vertexptr = (vertex) (otri).tri[plus1mod3[(otri).orient] + 3]
  1068. #define dest(otri, vertexptr)                                                 
  1069.   vertexptr = (vertex) (otri).tri[minus1mod3[(otri).orient] + 3]
  1071. #define apex(otri, vertexptr)                                                 
  1072.   vertexptr = (vertex) (otri).tri[(otri).orient + 3]
  1074. #define setorg(otri, vertexptr)                                               
  1075.   (otri).tri[plus1mod3[(otri).orient] + 3] = (triangle) vertexptr
  1077. #define setdest(otri, vertexptr)                                              
  1078.   (otri).tri[minus1mod3[(otri).orient] + 3] = (triangle) vertexptr
  1080. #define setapex(otri, vertexptr)                                              
  1081.   (otri).tri[(otri).orient + 3] = (triangle) vertexptr
  1083. /* Bond two triangles together.                                              */
  1085. #define bond(otri1, otri2)                                                    
  1086.   (otri1).tri[(otri1).orient] = encode(otri2);                                
  1087.   (otri2).tri[(otri2).orient] = encode(otri1)
  1089. /* Dissolve a bond (from one side).  Note that the other triangle will still */
  1090. /*   think it's connected to this triangle.  Usually, however, the other     */
  1091. /*   triangle is being deleted entirely, or bonded to another triangle, so   */
  1092. /*   it doesn't matter.                                                      */
  1094. #define dissolve(otri)                                                        
  1095.   (otri).tri[(otri).orient] = (triangle) m->dummytri
  1097. /* Copy an oriented triangle.                                                */
  1099. #define otricopy(otri1, otri2)                                                
  1100.   (otri2).tri = (otri1).tri;                                                  
  1101.   (otri2).orient = (otri1).orient
  1103. /* Test for equality of oriented triangles.                                  */
  1105. #define otriequal(otri1, otri2)                                               
  1106.   (((otri1).tri == (otri2).tri) &&                                            
  1107.    ((otri1).orient == (otri2).orient))
  1109. /* Primitives to infect or cure a triangle with the virus.  These rely on    */
  1110. /*   the assumption that all subsegments are aligned to four-byte boundaries.*/
  1112. #define infect(otri)                                                          
  1113.   (otri).tri[6] = (triangle)                                                  
  1114.                     ((unsigned long) (otri).tri[6] | (unsigned long) 2l)
  1116. #define uninfect(otri)                                                        
  1117.   (otri).tri[6] = (triangle)                                                  
  1118.                     ((unsigned long) (otri).tri[6] & ~ (unsigned long) 2l)
  1120. /* Test a triangle for viral infection.                                      */
  1122. #define infected(otri)                                                        
  1123.   (((unsigned long) (otri).tri[6] & (unsigned long) 2l) != 0l)
  1125. /* Check or set a triangle's attributes.                                     */
  1127. #define elemattribute(otri, attnum)                                           
  1128.   ((REAL *) (otri).tri)[m->elemattribindex + (attnum)]
  1130. #define setelemattribute(otri, attnum, value)                                 
  1131.   ((REAL *) (otri).tri)[m->elemattribindex + (attnum)] = value
  1133. /* Check or set a triangle's maximum area bound.                             */
  1135. #define areabound(otri)  ((REAL *) (otri).tri)[m->areaboundindex]
  1137. #define setareabound(otri, value)                                             
  1138.   ((REAL *) (otri).tri)[m->areaboundindex] = value
  1140. /* Check or set a triangle's deallocation.  Its second pointer is set to     */
  1141. /*   NULL to indicate that it is not allocated.  (Its first pointer is used  */
  1142. /*   for the stack of dead items.)  Its fourth pointer (its first vertex)    */
  1143. /*   is set to NULL in case a `badtriang' structure points to it.            */
  1145. #define deadtri(tria)  ((tria)[1] == (triangle) NULL)
  1147. #define killtri(tria)                                                         
  1148.   (tria)[1] = (triangle) NULL;                                                
  1149.   (tria)[3] = (triangle) NULL
  1151. /********* Primitives for subsegments                                *********/
  1152. /*                                                                           */
  1153. /*                                                                           */
  1155. /* sdecode() converts a pointer to an oriented subsegment.  The orientation  */
  1156. /*   is extracted from the least significant bit of the pointer.  The two    */
  1157. /*   least significant bits (one for orientation, one for viral infection)   */
  1158. /*   are masked out to produce the real pointer.                             */
  1160. #define sdecode(sptr, osub)                                                   
  1161.   (osub).ssorient = (int) ((unsigned long) (sptr) & (unsigned long) 1l);      
  1162.   (osub).ss = (subseg *)                                                      
  1163.               ((unsigned long) (sptr) & ~ (unsigned long) 3l)
  1165. /* sencode() compresses an oriented subsegment into a single pointer.  It    */
  1166. /*   relies on the assumption that all subsegments are aligned to two-byte   */
  1167. /*   boundaries, so the least significant bit of (osub).ss is zero.          */
  1169. #define sencode(osub)                                                         
  1170.   (subseg) ((unsigned long) (osub).ss | (unsigned long) (osub).ssorient)
  1172. /* ssym() toggles the orientation of a subsegment.                           */
  1174. #define ssym(osub1, osub2)                                                    
  1175.   (osub2).ss = (osub1).ss;                                                    
  1176.   (osub2).ssorient = 1 - (osub1).ssorient
  1178. #define ssymself(osub)                                                        
  1179.   (osub).ssorient = 1 - (osub).ssorient
  1181. /* spivot() finds the other subsegment (from the same segment) that shares   */
  1182. /*   the same origin.                                                        */
  1184. #define spivot(osub1, osub2)                                                  
  1185.   sptr = (osub1).ss[(osub1).ssorient];                                        
  1186.   sdecode(sptr, osub2)
  1188. #define spivotself(osub)                                                      
  1189.   sptr = (osub).ss[(osub).ssorient];                                          
  1190.   sdecode(sptr, osub)
  1192. /* snext() finds the next subsegment (from the same segment) in sequence;    */
  1193. /*   one whose origin is the input subsegment's destination.                 */
  1195. #define snext(osub1, osub2)                                                   
  1196.   sptr = (osub1).ss[1 - (osub1).ssorient];                                    
  1197.   sdecode(sptr, osub2)
  1199. #define snextself(osub)                                                       
  1200.   sptr = (osub).ss[1 - (osub).ssorient];                                      
  1201.   sdecode(sptr, osub)
  1203. /* These primitives determine or set the origin or destination of a          */
  1204. /*   subsegment or the segment that includes it.                             */
  1206. #define sorg(osub, vertexptr)                                                 
  1207.   vertexptr = (vertex) (osub).ss[2 + (osub).ssorient]
  1209. #define sdest(osub, vertexptr)                                                
  1210.   vertexptr = (vertex) (osub).ss[3 - (osub).ssorient]
  1212. #define setsorg(osub, vertexptr)                                              
  1213.   (osub).ss[2 + (osub).ssorient] = (subseg) vertexptr
  1215. #define setsdest(osub, vertexptr)                                             
  1216.   (osub).ss[3 - (osub).ssorient] = (subseg) vertexptr
  1218. #define segorg(osub, vertexptr)                                               
  1219.   vertexptr = (vertex) (osub).ss[4 + (osub).ssorient]
  1221. #define segdest(osub, vertexptr)                                              
  1222.   vertexptr = (vertex) (osub).ss[5 - (osub).ssorient]
  1224. #define setsegorg(osub, vertexptr)                                            
  1225.   (osub).ss[4 + (osub).ssorient] = (subseg) vertexptr
  1227. #define setsegdest(osub, vertexptr)                                           
  1228.   (osub).ss[5 - (osub).ssorient] = (subseg) vertexptr
  1230. /* These primitives read or set a boundary marker.  Boundary markers are     */
  1231. /*   used to hold user-defined tags for setting boundary conditions in       */
  1232. /*   finite element solvers.                                                 */
  1234. #define mark(osub)  (* (int *) ((osub).ss + 8))
  1236. #define setmark(osub, value)                                                  
  1237.   * (int *) ((osub).ss + 8) = value
  1239. /* Bond two subsegments together.                                            */
  1241. #define sbond(osub1, osub2)                                                   
  1242.   (osub1).ss[(osub1).ssorient] = sencode(osub2);                              
  1243.   (osub2).ss[(osub2).ssorient] = sencode(osub1)
  1245. /* Dissolve a subsegment bond (from one side).  Note that the other          */
  1246. /*   subsegment will still think it's connected to this subsegment.          */
  1248. #define sdissolve(osub)                                                       
  1249.   (osub).ss[(osub).ssorient] = (subseg) m->dummysub
  1251. /* Copy a subsegment.                                                        */
  1253. #define subsegcopy(osub1, osub2)                                              
  1254.   (osub2).ss = (osub1).ss;                                                    
  1255.   (osub2).ssorient = (osub1).ssorient
  1257. /* Test for equality of subsegments.                                         */
  1259. #define subsegequal(osub1, osub2)                                             
  1260.   (((osub1).ss == (osub2).ss) &&                                              
  1261.    ((osub1).ssorient == (osub2).ssorient))
  1263. /* Check or set a subsegment's deallocation.  Its second pointer is set to   */
  1264. /*   NULL to indicate that it is not allocated.  (Its first pointer is used  */
  1265. /*   for the stack of dead items.)  Its third pointer (its first vertex)     */
  1266. /*   is set to NULL in case a `badsubseg' structure points to it.            */
  1268. #define deadsubseg(sub)  ((sub)[1] == (subseg) NULL)
  1270. #define killsubseg(sub)                                                       
  1271.   (sub)[1] = (subseg) NULL;                                                   
  1272.   (sub)[2] = (subseg) NULL
  1274. /********* Primitives for interacting triangles and subsegments      *********/
  1275. /*                                                                           */
  1276. /*                                                                           */
  1278. /* tspivot() finds a subsegment abutting a triangle.                         */
  1280. #define tspivot(otri, osub)                                                   
  1281.   sptr = (subseg) (otri).tri[6 + (otri).orient];                              
  1282.   sdecode(sptr, osub)
  1284. /* stpivot() finds a triangle abutting a subsegment.  It requires that the   */
  1285. /*   variable `ptr' of type `triangle' be defined.                           */
  1287. #define stpivot(osub, otri)                                                   
  1288.   ptr = (triangle) (osub).ss[6 + (osub).ssorient];                            
  1289.   decode(ptr, otri)
  1291. /* Bond a triangle to a subsegment.                                          */
  1293. #define tsbond(otri, osub)                                                    
  1294.   (otri).tri[6 + (otri).orient] = (triangle) sencode(osub);                   
  1295.   (osub).ss[6 + (osub).ssorient] = (subseg) encode(otri)
  1297. /* Dissolve a bond (from the triangle side).                                 */
  1299. #define tsdissolve(otri)                                                      
  1300.   (otri).tri[6 + (otri).orient] = (triangle) m->dummysub
  1302. /* Dissolve a bond (from the subsegment side).                               */
  1304. #define stdissolve(osub)                                                      
  1305.   (osub).ss[6 + (osub).ssorient] = (subseg) m->dummytri
  1307. /********* Primitives for vertices                                   *********/
  1308. /*                                                                           */
  1309. /*                                                                           */
  1311. #define vertexmark(vx)  ((int *) (vx))[m->vertexmarkindex]
  1313. #define setvertexmark(vx, value)                                              
  1314.   ((int *) (vx))[m->vertexmarkindex] = value
  1316. #define vertextype(vx)  ((int *) (vx))[m->vertexmarkindex + 1]
  1318. #define setvertextype(vx, value)                                              
  1319.   ((int *) (vx))[m->vertexmarkindex + 1] = value
  1321. #define vertex2tri(vx)  ((triangle *) (vx))[m->vertex2triindex]
  1323. #define setvertex2tri(vx, value)                                              
  1324.   ((triangle *) (vx))[m->vertex2triindex] = value
  1326. /**                                                                         **/
  1327. /**                                                                         **/
  1328. /********* Mesh manipulation primitives end here                     *********/
  1330. /********* User-defined triangle evaluation routine begins here      *********/
  1331. /**                                                                         **/
  1332. /**                                                                         **/
  1334. /*****************************************************************************/
  1335. /*                                                                           */
  1336. /*  triunsuitable()   Determine if a triangle is unsuitable, and thus must   */
  1337. /*                    be further refined.                                    */
  1338. /*                                                                           */
  1339. /*  You may write your own procedure that decides whether or not a selected  */
  1340. /*  triangle is too big (and needs to be refined).  There are two ways to do */
  1341. /*  this.                                                                    */
  1342. /*                                                                           */
  1343. /*  (1)  Modify the procedure `triunsuitable' below, then recompile          */
  1344. /*  Triangle.                                                                */
  1345. /*                                                                           */
  1346. /*  (2)  Define the symbol EXTERNAL_TEST (either by adding the definition    */
  1347. /*  to this file, or by using the appropriate compiler switch).  This way,   */
  1348. /*  you can compile triangle.c separately from your test.  Write your own    */
  1349. /*  `triunsuitable' procedure in a separate C file (using the same prototype */
  1350. /*  as below).  Compile it and link the object code with triangle.o.         */
  1351. /*                                                                           */
  1352. /*  This procedure returns 1 if the triangle is too large and should be      */
  1353. /*  refined; 0 otherwise.                                                    */
  1354. /*                                                                           */
  1355. /*****************************************************************************/
  1357. #ifdef EXTERNAL_TEST
  1359. int triunsuitable();
  1361. #else /* not EXTERNAL_TEST */
  1363. #ifdef ANSI_DECLARATORS
  1364. int triunsuitable(vertex triorg, vertex tridest, vertex triapex, REAL area)
  1365. #else /* not ANSI_DECLARATORS */
  1366. int triunsuitable(triorg, tridest, triapex, area)
  1367. vertex triorg;                              /* The triangle's origin vertex. */
  1368. vertex tridest;                        /* The triangle's destination vertex. */
  1369. vertex triapex;                               /* The triangle's apex vertex. */
  1370. REAL area;                                      /* The area of the triangle. */
  1371. #endif /* not ANSI_DECLARATORS */
  1373. {
  1374.   REAL dxoa, dxda, dxod;
  1375.   REAL dyoa, dyda, dyod;
  1376.   REAL oalen, dalen, odlen;
  1377.   REAL maxlen;
  1379.   dxoa = triorg[0] - triapex[0];
  1380.   dyoa = triorg[1] - triapex[1];
  1381.   dxda = tridest[0] - triapex[0];
  1382.   dyda = tridest[1] - triapex[1];
  1383.   dxod = triorg[0] - tridest[0];
  1384.   dyod = triorg[1] - tridest[1];
  1385.   /* Find the squares of the lengths of the triangle's three edges. */
  1386.   oalen = dxoa * dxoa + dyoa * dyoa;
  1387.   dalen = dxda * dxda + dyda * dyda;
  1388.   odlen = dxod * dxod + dyod * dyod;
  1389.   /* Find the square of the length of the longest edge. */
  1390.   maxlen = (dalen > oalen) ? dalen : oalen;
  1391.   maxlen = (odlen > maxlen) ? odlen : maxlen;
  1393.   if (maxlen > 0.05 * (triorg[0] * triorg[0] + triorg[1] * triorg[1]) + 0.02) {
  1394.     return 1;
  1395.   } else {
  1396.     return 0;
  1397.   }
  1398. }
  1400. #endif /* not EXTERNAL_TEST */
  1402. /**                                                                         **/
  1403. /**                                                                         **/
  1404. /********* User-defined triangle evaluation routine ends here        *********/
  1406. /********* Memory allocation and program exit wrappers begin here    *********/
  1407. /**                                                                         **/
  1408. /**                                                                         **/
  1410. #ifdef ANSI_DECLARATORS
  1411. void triexit(int status)
  1412. #else /* not ANSI_DECLARATORS */
  1413. void triexit(status)
  1414. int status;
  1415. #endif /* not ANSI_DECLARATORS */
  1417. {
  1418.   exit(status);
  1419. }
  1421. #ifdef ANSI_DECLARATORS
  1422. VOID *trimalloc(int size)
  1423. #else /* not ANSI_DECLARATORS */
  1424. VOID *trimalloc(size)
  1425. int size;
  1426. #endif /* not ANSI_DECLARATORS */
  1428. {
  1429.   VOID *memptr;
  1431.   memptr = (VOID *) malloc((unsigned int) size);
  1432.   if (memptr == (VOID *) NULL) {
  1433.     printf("Error:  Out of memory.n");
  1434.     triexit(1);
  1435.   }
  1436.   return(memptr);
  1437. }
  1439. #ifdef ANSI_DECLARATORS
  1440. void trifree(VOID *memptr)
  1441. #else /* not ANSI_DECLARATORS */
  1442. void trifree(memptr)
  1443. VOID *memptr;
  1444. #endif /* not ANSI_DECLARATORS */
  1446. {
  1447.   free(memptr);
  1448. }
  1450. /**                                                                         **/
  1451. /**                                                                         **/
  1452. /********* Memory allocation and program exit wrappers end here      *********/
  1454. /********* User interaction routines begin here                      *********/
  1455. /**                                                                         **/
  1456. /**                                                                         **/
  1458. /*****************************************************************************/
  1459. /*                                                                           */
  1460. /*  syntax()   Print list of command line switches.                          */
  1461. /*                                                                           */
  1462. /*****************************************************************************/
  1464. #ifndef TRILIBRARY
  1466. void syntax()
  1467. {
  1468. #ifdef CDT_ONLY
  1469. #ifdef REDUCED
  1470.   printf("triangle [-pAcjevngBPNEIOXzo_lQVh] input_filen");
  1471. #else /* not REDUCED */
  1472.   printf("triangle [-pAcjevngBPNEIOXzo_iFlCQVh] input_filen");
  1473. #endif /* not REDUCED */
  1474. #else /* not CDT_ONLY */
  1475. #ifdef REDUCED
  1476.   printf("triangle [-prq__a__uAcDjevngBPNEIOXzo_YS__lQVh] input_filen");
  1477. #else /* not REDUCED */
  1478.   printf("triangle [-prq__a__uAcDjevngBPNEIOXzo_YS__iFlsCQVh] input_filen");
  1479. #endif /* not REDUCED */
  1480. #endif /* not CDT_ONLY */
  1482.   printf("    -p  Triangulates a Planar Straight Line Graph (.poly file).n");
  1483. #ifndef CDT_ONLY
  1484.   printf("    -r  Refines a previously generated mesh.n");
  1485.   printf(
  1486.     "    -q  Quality mesh generation.  A minimum angle may be specified.n");
  1487.   printf("    -a  Applies a maximum triangle area constraint.n");
  1488.   printf("    -u  Applies a user-defined triangle constraint.n");
  1489. #endif /* not CDT_ONLY */
  1490.   printf(
  1491.     "    -A  Applies attributes to identify triangles in certain regions.n");
  1492.   printf("    -c  Encloses the convex hull with segments.n");
  1493. #ifndef CDT_ONLY
  1494.   printf("    -D  Conforming Delaunay:  all triangles are truly Delaunay.n");
  1495. #endif /* not CDT_ONLY */
  1496. /*
  1497.   printf("    -w  Weighted Delaunay triangulation.n");
  1498.   printf("    -W  Regular triangulation (lower hull of a height field).n");
  1499. */
  1500.   printf("    -j  Jettison unused vertices from output .node file.n");
  1501.   printf("    -e  Generates an edge list.n");
  1502.   printf("    -v  Generates a Voronoi diagram.n");
  1503.   printf("    -n  Generates a list of triangle neighbors.n");
  1504.   printf("    -g  Generates an .off file for Geomview.n");
  1505.   printf("    -B  Suppresses output of boundary information.n");
  1506.   printf("    -P  Suppresses output of .poly file.n");
  1507.   printf("    -N  Suppresses output of .node file.n");
  1508.   printf("    -E  Suppresses output of .ele file.n");
  1509.   printf("    -I  Suppresses mesh iteration numbers.n");
  1510.   printf("    -O  Ignores holes in .poly file.n");
  1511.   printf("    -X  Suppresses use of exact arithmetic.n");
  1512.   printf("    -z  Numbers all items starting from zero (rather than one).n");
  1513.   printf("    -o2 Generates second-order subparametric elements.n");
  1514. #ifndef CDT_ONLY
  1515.   printf("    -Y  Suppresses boundary segment splitting.n");
  1516.   printf("    -S  Specifies maximum number of added Steiner points.n");
  1517. #endif /* not CDT_ONLY */
  1518. #ifndef REDUCED
  1519.   printf("    -i  Uses incremental method, rather than divide-and-conquer.n");
  1520.   printf("    -F  Uses Fortune's sweepline algorithm, rather than d-and-c.n");
  1521. #endif /* not REDUCED */
  1522.   printf("    -l  Uses vertical cuts only, rather than alternating cuts.n");
  1523. #ifndef REDUCED
  1524. #ifndef CDT_ONLY
  1525.   printf(
  1526.     "    -s  Force segments into mesh by splitting (instead of using CDT).n");
  1527. #endif /* not CDT_ONLY */
  1528.   printf("    -C  Check consistency of final mesh.n");
  1529. #endif /* not REDUCED */
  1530.   printf("    -Q  Quiet:  No terminal output except errors.n");
  1531.   printf("    -V  Verbose:  Detailed information on what I'm doing.n");
  1532.   printf("    -h  Help:  Detailed instructions for Triangle.n");
  1533.   triexit(0);
  1534. }
  1536. #endif /* not TRILIBRARY */
  1538. /*****************************************************************************/
  1539. /*                                                                           */
  1540. /*  info()   Print out complete instructions.                                */
  1541. /*                                                                           */
  1542. /*****************************************************************************/
  1544. #ifndef TRILIBRARY
  1546. void info()
  1547. {
  1548.   printf("Trianglen");
  1549.   printf(
  1550. "A Two-Dimensional Quality Mesh Generator and Delaunay Triangulator.n");
  1551.   printf("Version 1.6nn");
  1552.   printf(
  1553. "Copyright 1993, 1995, 1997, 1998, 2002, 2005 Jonathan Richard Shewchukn");
  1554.   printf("2360 Woolsey #H / Berkeley, California 94705-1927n");
  1555.   printf("Bugs/comments to jrs@cs.berkeley.edun");
  1556.   printf(
  1557. "Created as part of the Quake project (tools for earthquake simulation).n");
  1558.   printf(
  1559. "Supported in part by NSF Grant CMS-9318163 and an NSERC 1967 Scholarship.n");
  1560.   printf("There is no warranty whatsoever.  Use at your own risk.n");
  1561. #ifdef SINGLE
  1562.   printf("This executable is compiled for single precision arithmetic.nnn");
  1563. #else /* not SINGLE */
  1564.   printf("This executable is compiled for double precision arithmetic.nnn");
  1565. #endif /* not SINGLE */
  1566.   printf(
  1567. "Triangle generates exact Delaunay triangulations, constrained Delaunayn");
  1568.   printf(
  1569. "triangulations, conforming Delaunay triangulations, Voronoi diagrams, andn");
  1570.   printf(
  1571. "high-quality triangular meshes.  The latter can be generated with no smalln"
  1572. );
  1573.   printf(
  1574. "or large angles, and are thus suitable for finite element analysis.  If non"
  1575. );
  1576.   printf(
  1577. "command line switch is specified, your .node input file is read, and then");
  1578.   printf(
  1579. "Delaunay triangulation is returned in .node and .ele output files.  Then");
  1580.   printf("command syntax is:nn");
  1581.   printf("triangle [-prq__a__uAcDjevngBPNEIOXzo_YS__iFlsCQVh] input_filenn");
  1582.   printf(
  1583. "Underscores indicate that numbers may optionally follow certain switches.n");
  1584.   printf(
  1585. "Do not leave any space between a switch and its numeric parameter.n");
  1586.   printf(
  1587. "input_file must be a file with extension .node, or extension .poly if then");
  1588.   printf(
  1589. "-p switch is used.  If -r is used, you must supply .node and .ele files,n");
  1590.   printf(
  1591. "and possibly a .poly file and an .area file as well.  The formats of thesen"
  1592. );
  1593.   printf("files are described below.nn");
  1594.   printf("Command Line Switches:nn");
  1595.   printf(
  1596. "    -p  Reads a Planar Straight Line Graph (.poly file), which can specifyn"
  1597. );
  1598.   printf(
  1599. "        vertices, segments, holes, regional attributes, and regional arean");
  1600.   printf(
  1601. "        constraints.  Generates a constrained Delaunay triangulation (CDT)n"
  1602. );
  1603.   printf(
  1604. "        fitting the input; or, if -s, -q, -a, or -u is used, a conformingn");
  1605.   printf(
  1606. "        constrained Delaunay triangulation (CCDT).  If you want a trulyn");
  1607.   printf(
  1608. "        Delaunay (not just constrained Delaunay) triangulation, use -D asn");
  1609.   printf(
  1610. "        well.  When -p is not used, Triangle reads a .node file by default.n"
  1611. );
  1612.   printf(
  1613. "    -r  Refines a previously generated mesh.  The mesh is read from a .noden"
  1614. );
  1615.   printf(
  1616. "        file and an .ele file.  If -p is also used, a .poly file is readn");
  1617.   printf(
  1618. "        and used to constrain segments in the mesh.  If -a is also usedn");
  1619.   printf(
  1620. "        (with no number following), an .area file is read and used ton");
  1621.   printf(
  1622. "        impose area constraints on the mesh.  Further details on refinementn"
  1623. );
  1624.   printf("        appear below.n");
  1625.   printf(
  1626. "    -q  Quality mesh generation by Delaunay refinement (a hybrid of Pauln");
  1627.   printf(
  1628. "        Chew's and Jim Ruppert's algorithms).  Adds vertices to the mesh ton"
  1629. );
  1630.   printf(
  1631. "        ensure that all angles are between 20 and 140 degrees.  Ann");
  1632.   printf(
  1633. "        alternative bound on the minimum angle, replacing 20 degrees, mayn");
  1634.   printf(
  1635. "        be specified after the `q'.  The specified angle may include an");
  1636.   printf(
  1637. "        decimal point, but not exponential notation.  Note that a bound ofn"
  1638. );
  1639.   printf(
  1640. "        theta degrees on the smallest angle also implies a bound ofn");
  1641.   printf(
  1642. "        (180 - 2 theta) on the largest angle.  If the minimum angle is 28.6n"
  1643. );
  1644.   printf(
  1645. "        degrees or smaller, Triangle is mathematically guaranteed ton");
  1646.   printf(
  1647. "        terminate (assuming infinite precision arithmetic--Triangle mayn");
  1648.   printf(
  1649. "        fail to terminate if you run out of precision).  In practice,n");
  1650.   printf(
  1651. "        Triangle often succeeds for minimum angles up to 34 degrees.  Forn");
  1652.   printf(
  1653. "        some meshes, however, you might need to reduce the minimum angle ton"
  1654. );
  1655.   printf(
  1656. "        avoid problems associated with insufficient floating-pointn");
  1657.   printf("        precision.n");
  1658.   printf(
  1659. "    -a  Imposes a maximum triangle area.  If a number follows the `a', non");
  1660.   printf(
  1661. "        triangle is generated whose area is larger than that number.  If non"
  1662. );
  1663.   printf(
  1664. "        number is specified, an .area file (if -r is used) or .poly filen");
  1665.   printf(
  1666. "        (if -r is not used) specifies a set of maximum area constraints.n");
  1667.   printf(
  1668. "        An .area file contains a separate area constraint for eachn");
  1669.   printf(
  1670. "        triangle, and is useful for refining a finite element mesh based onn"
  1671. );
  1672.   printf(
  1673. "        a posteriori error estimates.  A .poly file can optionally containn"
  1674. );
  1675.   printf(
  1676. "        an area constraint for each segment-bounded region, therebyn");
  1677.   printf(
  1678. "        controlling triangle densities in a first triangulation of a PSLG.n"
  1679. );
  1680.   printf(
  1681. "        You can impose both a fixed area constraint and a varying arean");
  1682.   printf(
  1683. "        constraint by invoking the -a switch twice, once with and oncen");
  1684.   printf(
  1685. "        without a number following.  Each area specified may include an");
  1686.   printf("        decimal point.n");
  1687.   printf(
  1688. "    -u  Imposes a user-defined constraint on triangle size.  There are twon"
  1689. );
  1690.   printf(
  1691. "        ways to use this feature.  One is to edit the triunsuitable()n");
  1692.   printf(
  1693. "        procedure in triangle.c to encode any constraint you like, thenn");
  1694.   printf(
  1695. "        recompile Triangle.  The other is to compile triangle.c with then");
  1696.   printf(
  1697. "        EXTERNAL_TEST symbol set (compiler switch -DEXTERNAL_TEST), thenn");
  1698.   printf(
  1699. "        link Triangle with a separate object file that implementsn");
  1700.   printf(
  1701. "        triunsuitable().  In either case, the -u switch causes the user-n");
  1702.   printf("        defined test to be applied to every triangle.n");
  1703.   printf(
  1704. "    -A  Assigns an additional floating-point attribute to each trianglen");
  1705.   printf(
  1706. "        that identifies what segment-bounded region each triangle belongsn");
  1707.   printf(
  1708. "        to.  Attributes are assigned to regions by the .poly file.  If an");
  1709.   printf(
  1710. "        region is not explicitly marked by the .poly file, triangles inn");
  1711.   printf(
  1712. "        that region are assigned an attribute of zero.  The -A switch hasn");
  1713.   printf(
  1714. "        an effect only when the -p switch is used and the -r switch is not.n"
  1715. );
  1716.   printf(
  1717. "    -c  Creates segments on the convex hull of the triangulation.  If youn");
  1718.   printf(
  1719. "        are triangulating a vertex set, this switch causes a .poly file ton"
  1720. );
  1721.   printf(
  1722. "        be written, containing all edges of the convex hull.  If you aren");
  1723.   printf(
  1724. "        triangulating a PSLG, this switch specifies that the whole convexn");
  1725.   printf(
  1726. "        hull of the PSLG should be triangulated, regardless of whatn");
  1727.   printf(
  1728. "        segments the PSLG has.  If you do not use this switch whenn");
  1729.   printf(
  1730. "        triangulating a PSLG, Triangle assumes that you have identified then"
  1731. );
  1732.   printf(
  1733. "        region to be triangulated by surrounding it with segments of then");
  1734.   printf(
  1735. "        input PSLG.  Beware:  if you are not careful, this switch can causen"
  1736. );
  1737.   printf(
  1738. "        the introduction of an extremely thin angle between a PSLG segmentn"
  1739. );
  1740.   printf(
  1741. "        and a convex hull segment, which can cause overrefinement (andn");
  1742.   printf(
  1743. "        possibly failure if Triangle runs out of precision).  If you aren");
  1744.   printf(
  1745. "        refining a mesh, the -c switch works differently:  it causes an");
  1746.   printf(
  1747. "        .poly file to be written containing the boundary edges of the meshn"
  1748. );
  1749.   printf("        (useful if no .poly file was read).n");
  1750.   printf(
  1751. "    -D  Conforming Delaunay triangulation:  use this switch if you want ton"
  1752. );
  1753.   printf(
  1754. "        ensure that all the triangles in the mesh are Delaunay, and notn");
  1755.   printf(
  1756. "        merely constrained Delaunay; or if you want to ensure that all then"
  1757. );
  1758.   printf(
  1759. "        Voronoi vertices lie within the triangulation.  (Some finite volumen"
  1760. );
  1761.   printf(
  1762. "        methods have this requirement.)  This switch invokes Ruppert'sn");
  1763.   printf(
  1764. "        original algorithm, which splits every subsegment whose diametraln");
  1765.   printf(
  1766. "        circle is encroached.  It usually increases the number of verticesn"
  1767. );
  1768.   printf("        and triangles.n");
  1769.   printf(
  1770. "    -j  Jettisons vertices that are not part of the final triangulationn");
  1771.   printf(
  1772. "        from the output .node file.  By default, Triangle copies alln");
  1773.   printf(
  1774. "        vertices in the input .node file to the output .node file, in then");
  1775.   printf(
  1776. "        same order, so their indices do not change.  The -j switch preventsn"
  1777. );
  1778.   printf(
  1779. "        duplicated input vertices, or vertices `eaten' by holes, fromn");
  1780.   printf(
  1781. "        appearing in the output .node file.  Thus, if two input verticesn");
  1782.   printf(
  1783. "        have exactly the same coordinates, only the first appears in then");
  1784.   printf(
  1785. "        output.  If any vertices are jettisoned, the vertex numbering inn");
  1786.   printf(
  1787. "        the output .node file differs from that of the input .node file.n");
  1788.   printf(
  1789. "    -e  Outputs (to an .edge file) a list of edges of the triangulation.n");
  1790.   printf(
  1791. "    -v  Outputs the Voronoi diagram associated with the triangulation.n");
  1792.   printf(
  1793. "        Does not attempt to detect degeneracies, so some Voronoi verticesn");
  1794.   printf(
  1795. "        may be duplicated.  See the discussion of Voronoi diagrams below.n");
  1796.   printf(
  1797. "    -n  Outputs (to a .neigh file) a list of triangles neighboring eachn");
  1798.   printf("        triangle.n");
  1799.   printf(
  1800. "    -g  Outputs the mesh to an Object File Format (.off) file, suitable forn"
  1801. );
  1802.   printf("        viewing with the Geometry Center's Geomview package.n");
  1803.   printf(
  1804. "    -B  No boundary markers in the output .node, .poly, and .edge outputn");
  1805.   printf(
  1806. "        files.  See the detailed discussion of boundary markers below.n");
  1807.   printf(
  1808. "    -P  No output .poly file.  Saves disk space, but you lose the abilityn");
  1809.   printf(
  1810. "        to maintain constraining segments on later refinements of the mesh.n"
  1811. );
  1812.   printf("    -N  No output .node file.n");
  1813.   printf("    -E  No output .ele file.n");
  1814.   printf(
  1815. "    -I  No iteration numbers.  Suppresses the output of .node and .polyn");
  1816.   printf(
  1817. "        files, so your input files won't be overwritten.  (If your input isn"
  1818. );
  1819.   printf(
  1820. "        a .poly file only, a .node file is written.)  Cannot be used withn");
  1821.   printf(
  1822. "        the -r switch, because that would overwrite your input .ele file.n");
  1823.   printf(
  1824. "        Shouldn't be used with the -q, -a, -u, or -s switch if you aren");
  1825.   printf(
  1826. "        using a .node file for input, because no .node file is written, son"
  1827. );
  1828.   printf("        there is no record of any added Steiner points.n");
  1829.   printf("    -O  No holes.  Ignores the holes in the .poly file.n");
  1830.   printf(
  1831. "    -X  No exact arithmetic.  Normally, Triangle uses exact floating-pointn"
  1832. );
  1833.   printf(
  1834. "        arithmetic for certain tests if it thinks the inexact tests are notn"
  1835. );
  1836.   printf(
  1837. "        accurate enough.  Exact arithmetic ensures the robustness of then");
  1838.   printf(
  1839. "        triangulation algorithms, despite floating-point roundoff error.n");
  1840.   printf(
  1841. "        Disabling exact arithmetic with the -X switch causes a smalln");
  1842.   printf(
  1843. "        improvement in speed and creates the possibility that Triangle willn"
  1844. );
  1845.   printf("        fail to produce a valid mesh.  Not recommended.n");
  1846.   printf(
  1847. "    -z  Numbers all items starting from zero (rather than one).  Note thatn"
  1848. );
  1849.   printf(
  1850. "        this switch is normally overridden by the value used to number then"
  1851. );
  1852.   printf(
  1853. "        first vertex of the input .node or .poly file.  However, thisn");
  1854.   printf(
  1855. "        switch is useful when calling Triangle from another program.n");
  1856.   printf(
  1857. "    -o2 Generates second-order subparametric elements with six nodes each.n"
  1858. );
  1859.   printf(
  1860. "    -Y  No new vertices on the boundary.  This switch is useful when then");
  1861.   printf(
  1862. "        mesh boundary must be preserved so that it conforms to somen");
  1863.   printf(
  1864. "        adjacent mesh.  Be forewarned that you will probably sacrifice muchn"
  1865. );
  1866.   printf(
  1867. "        of the quality of the mesh; Triangle will try, but the resultingn");
  1868.   printf(
  1869. "        mesh may contain poorly shaped triangles.  Works well if all then");
  1870.   printf(
  1871. "        boundary vertices are closely spaced.  Specify this switch twicen");
  1872.   printf(
  1873. "        (`-YY') to prevent all segment splitting, including internaln");
  1874.   printf("        boundaries.n");
  1875.   printf(
  1876. "    -S  Specifies the maximum number of Steiner points (vertices that aren");
  1877.   printf(
  1878. "        not in the input, but are added to meet the constraints on minimumn"
  1879. );
  1880.   printf(
  1881. "        angle and maximum area).  The default is to allow an unlimitedn");
  1882.   printf(
  1883. "        number.  If you specify this switch with no number after it,n");
  1884.   printf(
  1885. "        the limit is set to zero.  Triangle always adds vertices at segmentn"
  1886. );
  1887.   printf(
  1888. "        intersections, even if it needs to use more vertices than the limitn"
  1889. );
  1890.   printf(
  1891. "        you set.  When Triangle inserts segments by splitting (-s), itn");
  1892.   printf(
  1893. "        always adds enough vertices to ensure that all the segments of then"
  1894. );
  1895.   printf("        PLSG are recovered, ignoring the limit if necessary.n");
  1896.   printf(
  1897. "    -i  Uses an incremental rather than a divide-and-conquer algorithm ton");
  1898.   printf(
  1899. "        construct a Delaunay triangulation.  Try it if the divide-and-n");
  1900.   printf("        conquer algorithm fails.n");
  1901.   printf(
  1902. "    -F  Uses Steven Fortune's sweepline algorithm to construct a Delaunayn");
  1903.   printf(
  1904. "        triangulation.  Warning:  does not use exact arithmetic for alln");
  1905.   printf("        calculations.  An exact result is not guaranteed.n");
  1906.   printf(
  1907. "    -l  Uses only vertical cuts in the divide-and-conquer algorithm.  Byn");
  1908.   printf(
  1909. "        default, Triangle alternates between vertical and horizontal cuts,n"
  1910. );
  1911.   printf(
  1912. "        which usually improve the speed except with vertex sets that aren");
  1913.   printf(
  1914. "        small or short and wide.  This switch is primarily of theoreticaln");
  1915.   printf("        interest.n");
  1916.   printf(
  1917. "    -s  Specifies that segments should be forced into the triangulation byn"
  1918. );
  1919.   printf(
  1920. "        recursively splitting them at their midpoints, rather than byn");
  1921.   printf(
  1922. "        generating a constrained Delaunay triangulation.  Segment splittingn"
  1923. );
  1924.   printf(
  1925. "        is true to Ruppert's original algorithm, but can create needlesslyn"
  1926. );
  1927.   printf(
  1928. "        small triangles.  This switch is primarily of theoretical interest.n"
  1929. );
  1930.   printf(
  1931. "    -C  Check the consistency of the final mesh.  Uses exact arithmetic forn"
  1932. );
  1933.   printf(
  1934. "        checking, even if the -X switch is used.  Useful if you suspectn");
  1935.   printf("        Triangle is buggy.n");
  1936.   printf(
  1937. "    -Q  Quiet:  Suppresses all explanation of what Triangle is doing,n");
  1938.   printf("        unless an error occurs.n");
  1939.   printf(
  1940. "    -V  Verbose:  Gives detailed information about what Triangle is doing.n"
  1941. );
  1942.   printf(
  1943. "        Add more `V's for increasing amount of detail.  `-V' is mostn");
  1944.   printf(
  1945. "        useful; itgives information on algorithmic progress and much moren");
  1946.   printf(
  1947. "        detailed statistics.  `-VV' gives vertex-by-vertex details, andn");
  1948.   printf(
  1949. "        prints so much that Triangle runs much more slowly.  `-VVVV' givesn"
  1950. );
  1951.   printf("        information only a debugger could love.n");
  1952.   printf("    -h  Help:  Displays these instructions.n");
  1953.   printf("n");
  1954.   printf("Definitions:n");
  1955.   printf("n");
  1956.   printf(
  1957. "  A Delaunay triangulation of a vertex set is a triangulation whosen");
  1958.   printf(
  1959. "  vertices are the vertex set, that covers the convex hull of the vertexn");
  1960.   printf(
  1961. "  set.  A Delaunay triangulation has the property that no vertex liesn");
  1962.   printf(
  1963. "  inside the circumscribing circle (circle that passes through all threen");
  1964.   printf("  vertices) of any triangle in the triangulation.nn");
  1965.   printf(
  1966. "  A Voronoi diagram of a vertex set is a subdivision of the plane inton");
  1967.   printf(
  1968. "  polygonal cells (some of which may be unbounded, meaning infinitelyn");
  1969.   printf(
  1970. "  large), where each cell is the set of points in the plane that are closern"
  1971. );
  1972.   printf(
  1973. "  to some input vertex than to any other input vertex.  The Voronoi diagramn"
  1974. );
  1975.   printf("  is a geometric dual of the Delaunay triangulation.nn");
  1976.   printf(
  1977. "  A Planar Straight Line Graph (PSLG) is a set of vertices and segments.n");
  1978.   printf(
  1979. "  Segments are simply edges, whose endpoints are all vertices in the PSLG.n"
  1980. );
  1981.   printf(
  1982. "  Segments may intersect each other only at their endpoints.  The filen");
  1983.   printf("  format for PSLGs (.poly files) is described below.nn");
  1984.   printf(
  1985. "  A constrained Delaunay triangulation (CDT) of a PSLG is similar to an");
  1986.   printf(
  1987. "  Delaunay triangulation, but each PSLG segment is present as a single edgen"
  1988. );
  1989.   printf(
  1990. "  of the CDT.  (A constrained Delaunay triangulation is not truly an");
  1991.   printf(
  1992. "  Delaunay triangulation, because some of its triangles might not ben");
  1993.   printf(
  1994. "  Delaunay.)  By definition, a CDT does not have any vertices other thann");
  1995.   printf(
  1996. "  those specified in the input PSLG.  Depending on context, a CDT mightn");
  1997.   printf(
  1998. "  cover the convex hull of the PSLG, or it might cover only a segment-n");
  1999.   printf("  bounded region (e.g. a polygon).nn");
  2000.   printf(
  2001. "  A conforming Delaunay triangulation of a PSLG is a triangulation in whichn"
  2002. );
  2003.   printf(
  2004. "  each triangle is truly Delaunay, and each PSLG segment is represented byn"
  2005. );
  2006.   printf(
  2007. "  a linear contiguous sequence of edges of the triangulation.  New verticesn"
  2008. );
  2009.   printf(
  2010. "  (not part of the PSLG) may appear, and each input segment may have beenn");
  2011.   printf(
  2012. "  subdivided into shorter edges (subsegments) by these additional vertices.n"
  2013. );
  2014.   printf(
  2015. "  The new vertices are frequently necessary to maintain the Delaunayn");
  2016.   printf("  property while ensuring that every segment is represented.nn");
  2017.   printf(
  2018. "  A conforming constrained Delaunay triangulation (CCDT) of a PSLG is an");
  2019.   printf(
  2020. "  triangulation of a PSLG whose triangles are constrained Delaunay.  Newn");
  2021.   printf("  vertices may appear, and input segments may be subdivided inton");
  2022.   printf(
  2023. "  subsegments, but not to guarantee that segments are respected; rather, ton"
  2024. );
  2025.   printf(
  2026. "  improve the quality of the triangles.  The high-quality meshes producedn");
  2027.   printf(
  2028. "  by the -q switch are usually CCDTs, but can be made conforming Delaunayn");
  2029.   printf("  with the -D switch.nn");
  2030.   printf("File Formats:nn");
  2031.   printf(
  2032. "  All files may contain comments prefixed by the character '#'.  Vertices,n"
  2033. );
  2034.   printf(
  2035. "  triangles, edges, holes, and maximum area constraints must be numberedn");
  2036.   printf(
  2037. "  consecutively, starting from either 1 or 0.  Whichever you choose, alln");
  2038.   printf(
  2039. "  input files must be consistent; if the vertices are numbered from 1, son");
  2040.   printf(
  2041. "  must be all other objects.  Triangle automatically detects your choicen");
  2042.   printf(
  2043. "  while reading the .node (or .poly) file.  (When calling Triangle fromn");
  2044.   printf(
  2045. "  another program, use the -z switch if you wish to number objects fromn");
  2046.   printf("  zero.)  Examples of these file formats are given below.nn");
  2047.   printf("  .node files:n");
  2048.   printf(
  2049. "    First line:  <# of vertices> <dimension (must be 2)> <# of attributes>n"
  2050. );
  2051.   printf(
  2052. "                                           <# of boundary markers (0 or 1)>n"
  2053. );
  2054.   printf(
  2055. "    Remaining lines:  <vertex #> <x> <y> [attributes] [boundary marker]n");
  2056.   printf("n");
  2057.   printf(
  2058. "    The attributes, which are typically floating-point values of physicaln");
  2059.   printf(
  2060. "    quantities (such as mass or conductivity) associated with the nodes ofn"
  2061. );
  2062.   printf(
  2063. "    a finite element mesh, are copied unchanged to the output mesh.  If -q,n"
  2064. );
  2065.   printf(
  2066. "    -a, -u, -D, or -s is selected, each new Steiner point added to the meshn"
  2067. );
  2068.   printf("    has attributes assigned to it by linear interpolation.nn");
  2069.   printf(
  2070. "    If the fourth entry of the first line is `1', the last column of then");
  2071.   printf(
  2072. "    remainder of the file is assumed to contain boundary markers.  Boundaryn"
  2073. );
  2074.   printf(
  2075. "    markers are used to identify boundary vertices and vertices resting onn"
  2076. );
  2077.   printf(
  2078. "    PSLG segments; a complete description appears in a section below.  Then"
  2079. );
  2080.   printf(
  2081. "    .node file produced by Triangle contains boundary markers in the lastn");
  2082.   printf("    column unless they are suppressed by the -B switch.nn");
  2083.   printf("  .ele files:n");
  2084.   printf(
  2085. "    First line:  <# of triangles> <nodes per triangle> <# of attributes>n");
  2086.   printf(
  2087. "    Remaining lines:  <triangle #> <node> <node> <node> ... [attributes]n");
  2088.   printf("n");
  2089.   printf(
  2090. "    Nodes are indices into the corresponding .node file.  The first threen");
  2091.   printf(
  2092. "    nodes are the corner vertices, and are listed in counterclockwise ordern"
  2093. );
  2094.   printf(
  2095. "    around each triangle.  (The remaining nodes, if any, depend on the typen"
  2096. );
  2097.   printf("    of finite element used.)nn");
  2098.   printf(
  2099. "    The attributes are just like those of .node files.  Because there is non"
  2100. );
  2101.   printf(
  2102. "    simple mapping from input to output triangles, Triangle attempts ton");
  2103.   printf(
  2104. "    interpolate attributes, and may cause a lot of diffusion of attributesn"
  2105. );
  2106.   printf(
  2107. "    among nearby triangles as the triangulation is refined.  Attributes don"
  2108. );
  2109.   printf("    not diffuse across segments, so attributes used to identifyn");
  2110.   printf("    segment-bounded regions remain intact.nn");
  2111.   printf(
  2112. "    In .ele files produced by Triangle, each triangular element has threen");
  2113.   printf(
  2114. "    nodes (vertices) unless the -o2 switch is used, in which casen");
  2115.   printf(
  2116. "    subparametric quadratic elements with six nodes each are generated.n");
  2117.   printf(
  2118. "    The first three nodes are the corners in counterclockwise order, andn");
  2119.   printf(
  2120. "    the fourth, fifth, and sixth nodes lie on the midpoints of the edgesn");
  2121.   printf(
  2122. "    opposite the first, second, and third vertices, respectively.n");
  2123.   printf("n");
  2124.   printf("  .poly files:n");
  2125.   printf(
  2126. "    First line:  <# of vertices> <dimension (must be 2)> <# of attributes>n"
  2127. );
  2128.   printf(
  2129. "                                           <# of boundary markers (0 or 1)>n"
  2130. );
  2131.   printf(
  2132. "    Following lines:  <vertex #> <x> <y> [attributes] [boundary marker]n");
  2133.   printf("    One line:  <# of segments> <# of boundary markers (0 or 1)>n");
  2134.   printf(
  2135. "    Following lines:  <segment #> <endpoint> <endpoint> [boundary marker]n");
  2136.   printf("    One line:  <# of holes>n");
  2137.   printf("    Following lines:  <hole #> <x> <y>n");
  2138.   printf(
  2139. "    Optional line:  <# of regional attributes and/or area constraints>n");
  2140.   printf(
  2141. "    Optional following lines:  <region #> <x> <y> <attribute> <max area>n");
  2142.   printf("n");
  2143.   printf(
  2144. "    A .poly file represents a PSLG, as well as some additional information.n"
  2145. );
  2146.   printf(
  2147. "    The first section lists all the vertices, and is identical to then");
  2148.   printf(
  2149. "    format of .node files.  <# of vertices> may be set to zero to indicaten"
  2150. );
  2151.   printf(
  2152. "    that the vertices are listed in a separate .node file; .poly filesn");
  2153.   printf(
  2154. "    produced by Triangle always have this format.  A vertex set representedn"
  2155. );
  2156.   printf(
  2157. "    this way has the advantage that it may easily be triangulated with orn");
  2158.   printf(
  2159. "    without segments (depending on whether the -p switch is invoked).n");
  2160.   printf("n");
  2161.   printf(
  2162. "    The second section lists the segments.  Segments are edges whosen");
  2163.   printf(
  2164. "    presence in the triangulation is enforced.  (Depending on the choice ofn"
  2165. );
  2166.   printf(
  2167. "    switches, segment might be subdivided into smaller edges).  Eachn");
  2168.   printf(
  2169. "    segment is specified by listing the indices of its two endpoints.  Thisn"
  2170. );
  2171.   printf(
  2172. "    means that you must include its endpoints in the vertex list.  Eachn");
  2173.   printf("    segment, like each point, may have a boundary marker.nn");
  2174.   printf(
  2175. "    If -q, -a, -u, and -s are not selected, Triangle produces a constrainedn"
  2176. );
  2177.   printf(
  2178. "    Delaunay triangulation (CDT), in which each segment appears as a singlen"
  2179. );
  2180.   printf(
  2181. "    edge in the triangulation.  If -q, -a, -u, or -s is selected, Trianglen"
  2182. );
  2183.   printf(
  2184. "    produces a conforming constrained Delaunay triangulation (CCDT), inn");
  2185.   printf(
  2186. "    which segments may be subdivided into smaller edges.  If -D isn");
  2187.   printf(
  2188. "    selected, Triangle produces a conforming Delaunay triangulation, son");
  2189.   printf(
  2190. "    that every triangle is Delaunay, and not just constrained Delaunay.n");
  2191.   printf("n");
  2192.   printf(
  2193. "    The third section lists holes (and concavities, if -c is selected) inn");
  2194.   printf(
  2195. "    the triangulation.  Holes are specified by identifying a point insiden");
  2196.   printf(
  2197. "    each hole.  After the triangulation is formed, Triangle creates holesn");
  2198.   printf(
  2199. "    by eating triangles, spreading out from each hole point until itsn");
  2200.   printf(
  2201. "    progress is blocked by segments in the PSLG.  You must be careful ton");
  2202.   printf(
  2203. "    enclose each hole in segments, or your whole triangulation might ben");
  2204.   printf(
  2205. "    eaten away.  If the two triangles abutting a segment are eaten, then");
  2206.   printf(
  2207. "    segment itself is also eaten.  Do not place a hole directly on an");
  2208.   printf("    segment; if you do, Triangle chooses one side of the segmentn");
  2209.   printf("    arbitrarily.nn");
  2210.   printf(
  2211. "    The optional fourth section lists regional attributes (to be assignedn");
  2212.   printf(
  2213. "    to all triangles in a region) and regional constraints on the maximumn");
  2214.   printf(
  2215. "    triangle area.  Triangle reads this section only if the -A switch isn");
  2216.   printf(
  2217. "    used or the -a switch is used without a number following it, and the -rn"
  2218. );
  2219.   printf(
  2220. "    switch is not used.  Regional attributes and area constraints aren");
  2221.   printf(
  2222. "    propagated in the same manner as holes:  you specify a point for eachn");
  2223.   printf(
  2224. "    attribute and/or constraint, and the attribute and/or constraintn");
  2225.   printf(
  2226. "    affects the whole region (bounded by segments) containing the point.n");
  2227.   printf(
  2228. "    If two values are written on a line after the x and y coordinate, then");
  2229.   printf(
  2230. "    first such value is assumed to be a regional attribute (but is onlyn");
  2231.   printf(
  2232. "    applied if the -A switch is selected), and the second value is assumedn"
  2233. );
  2234.   printf(
  2235. "    to be a regional area constraint (but is only applied if the -a switchn"
  2236. );
  2237.   printf(
  2238. "    is selected).  You may specify just one value after the coordinates,n");
  2239.   printf(
  2240. "    which can serve as both an attribute and an area constraint, dependingn"
  2241. );
  2242.   printf(
  2243. "    on the choice of switches.  If you are using the -A and -a switchesn");
  2244.   printf(
  2245. "    simultaneously and wish to assign an attribute to some region withoutn");
  2246.   printf("    imposing an area constraint, use a negative maximum area.nn");
  2247.   printf(
  2248. "    When a triangulation is created from a .poly file, you must eithern");
  2249.   printf(
  2250. "    enclose the entire region to be triangulated in PSLG segments, orn");
  2251.   printf(
  2252. "    use the -c switch, which automatically creates extra segments thatn");
  2253.   printf(
  2254. "    enclose the convex hull of the PSLG.  If you do not use the -c switch,n"
  2255. );
  2256.   printf(
  2257. "    Triangle eats all triangles that are not enclosed by segments; if youn");
  2258.   printf(
  2259. "    are not careful, your whole triangulation may be eaten away.  If you don"
  2260. );
  2261.   printf(
  2262. "    use the -c switch, you can still produce concavities by the appropriaten"
  2263. );
  2264.   printf(
  2265. "    placement of holes just inside the boundary of the convex hull.n");
  2266.   printf("n");
  2267.   printf(
  2268. "    An ideal PSLG has no intersecting segments, nor any vertices that lien");
  2269.   printf(
  2270. "    upon segments (except, of course, the endpoints of each segment).  Youn"
  2271. );
  2272.   printf(
  2273. "    aren't required to make your .poly files ideal, but you should be awaren"
  2274. );
  2275.   printf(
  2276. "    of what can go wrong.  Segment intersections are relatively safe--n");
  2277.   printf(
  2278. "    Triangle calculates the intersection points for you and adds them ton");
  2279.   printf(
  2280. "    the triangulation--as long as your machine's floating-point precisionn");
  2281.   printf(
  2282. "    doesn't become a problem.  You are tempting the fates if you have threen"
  2283. );
  2284.   printf(
  2285. "    segments that cross at the same location, and expect Triangle to figuren"
  2286. );
  2287.   printf(
  2288. "    out where the intersection point is.  Thanks to floating-point roundoffn"
  2289. );
  2290.   printf(
  2291. "    error, Triangle will probably decide that the three segments intersectn"
  2292. );
  2293.   printf(
  2294. "    at three different points, and you will find a minuscule triangle inn");
  2295.   printf(
  2296. "    your output--unless Triangle tries to refine the tiny triangle, usesn");
  2297.   printf(
  2298. "    up the last bit of machine precision, and fails to terminate at all.n");
  2299.   printf(
  2300. "    You're better off putting the intersection point in the input files,n");
  2301.   printf(
  2302. "    and manually breaking up each segment into two.  Similarly, if youn");
  2303.   printf(
  2304. "    place a vertex at the middle of a segment, and hope that Triangle willn"
  2305. );
  2306.   printf(
  2307. "    break up the segment at that vertex, you might get lucky.  On the othern"
  2308. );
  2309.   printf(
  2310. "    hand, Triangle might decide that the vertex doesn't lie precisely onn");
  2311.   printf(
  2312. "    the segment, and you'll have a needle-sharp triangle in your output--orn"
  2313. );
  2314.   printf("    a lot of tiny triangles if you're generating a quality mesh.n");
  2315.   printf("n");
  2316.   printf(
  2317. "    When Triangle reads a .poly file, it also writes a .poly file, whichn");
  2318.   printf(
  2319. "    includes all the subsegments--the edges that are parts of inputn");
  2320.   printf(
  2321. "    segments.  If the -c switch is used, the output .poly file alson");
  2322.   printf(
  2323. "    includes all of the edges on the convex hull.  Hence, the output .polyn"
  2324. );
  2325.   printf(
  2326. "    file is useful for finding edges associated with input segments and forn"
  2327. );
  2328.   printf(
  2329. "    setting boundary conditions in finite element simulations.  Moreover,n");
  2330.   printf(
  2331. "    you will need the output .poly file if you plan to refine the outputn");
  2332.   printf(
  2333. "    mesh, and don't want segments to be missing in later triangulations.n");
  2334.   printf("n");
  2335.   printf("  .area files:n");
  2336.   printf("    First line:  <# of triangles>n");
  2337.   printf("    Following lines:  <triangle #> <maximum area>n");
  2338.   printf("n");
  2339.   printf(
  2340. "    An .area file associates with each triangle a maximum area that is usedn"
  2341. );
  2342.   printf(
  2343. "    for mesh refinement.  As with other file formats, every triangle mustn");
  2344.   printf(
  2345. "    be represented, and the triangles must be numbered consecutively.  An");
  2346.   printf(
  2347. "    triangle may be left unconstrained by assigning it a negative maximumn");
  2348.   printf("    area.nn");
  2349.   printf("  .edge files:n");
  2350.   printf("    First line:  <# of edges> <# of boundary markers (0 or 1)>n");
  2351.   printf(
  2352. "    Following lines:  <edge #> <endpoint> <endpoint> [boundary marker]n");
  2353.   printf("n");
  2354.   printf(
  2355. "    Endpoints are indices into the corresponding .node file.  Triangle cann"
  2356. );
  2357.   printf(
  2358. "    produce .edge files (use the -e switch), but cannot read them.  Then");
  2359.   printf(
  2360. "    optional column of boundary markers is suppressed by the -B switch.n");
  2361.   printf("n");
  2362.   printf(
  2363. "    In Voronoi diagrams, one also finds a special kind of edge that is ann");
  2364.   printf(
  2365. "    infinite ray with only one endpoint.  For these edges, a differentn");
  2366.   printf("    format is used:nn");
  2367.   printf("        <edge #> <endpoint> -1 <direction x> <direction y>nn");
  2368.   printf(
  2369. "    The `direction' is a floating-point vector that indicates the directionn"
  2370. );
  2371.   printf("    of the infinite ray.nn");
  2372.   printf("  .neigh files:n");
  2373.   printf(
  2374. "    First line:  <# of triangles> <# of neighbors per triangle (always 3)>n"
  2375. );
  2376.   printf(
  2377. "    Following lines:  <triangle #> <neighbor> <neighbor> <neighbor>n");
  2378.   printf("n");
  2379.   printf(
  2380. "    Neighbors are indices into the corresponding .ele file.  An index of -1n"
  2381. );
  2382.   printf(
  2383. "    indicates no neighbor (because the triangle is on an exteriorn");
  2384.   printf(
  2385. "    boundary).  The first neighbor of triangle i is opposite the firstn");
  2386.   printf("    corner of triangle i, and so on.nn");
  2387.   printf(
  2388. "    Triangle can produce .neigh files (use the -n switch), but cannot readn"
  2389. );
  2390.   printf("    them.nn");
  2391.   printf("Boundary Markers:nn");
  2392.   printf(
  2393. "  Boundary markers are tags used mainly to identify which output verticesn");
  2394.   printf(
  2395. "  and edges are associated with which PSLG segment, and to identify whichn");
  2396.   printf(
  2397. "  vertices and edges occur on a boundary of the triangulation.  A commonn");
  2398.   printf(
  2399. "  use is to determine where boundary conditions should be applied to an");
  2400.   printf(
  2401. "  finite element mesh.  You can prevent boundary markers from being writtenn"
  2402. );
  2403.   printf("  into files produced by Triangle by using the -B switch.nn");
  2404.   printf(
  2405. "  The boundary marker associated with each segment in an output .poly filen"
  2406. );
  2407.   printf("  and each edge in an output .edge file is chosen as follows:n");
  2408.   printf(
  2409. "    - If an output edge is part or all of a PSLG segment with a nonzeron");
  2410.   printf(
  2411. "      boundary marker, then the edge is assigned the same marker.n");
  2412.   printf(
  2413. "    - Otherwise, if the edge lies on a boundary of the triangulationn");
  2414.   printf(
  2415. "      (even the boundary of a hole), then the edge is assigned the markern");
  2416.   printf("      one (1).n");
  2417.   printf("    - Otherwise, the edge is assigned the marker zero (0).n");
  2418.   printf(
  2419. "  The boundary marker associated with each vertex in an output .node filen");
  2420.   printf("  is chosen as follows:n");
  2421.   printf(
  2422. "    - If a vertex is assigned a nonzero boundary marker in the input file,n"
  2423. );
  2424.   printf(
  2425. "      then it is assigned the same marker in the output .node file.n");
  2426.   printf(
  2427. "    - Otherwise, if the vertex lies on a PSLG segment (even if it is ann");
  2428.   printf(
  2429. "      endpoint of the segment) with a nonzero boundary marker, then then");
  2430.   printf(
  2431. "      vertex is assigned the same marker.  If the vertex lies on severaln");
  2432.   printf("      such segments, one of the markers is chosen arbitrarily.n");
  2433.   printf(
  2434. "    - Otherwise, if the vertex occurs on a boundary of the triangulation,n");
  2435.   printf("      then the vertex is assigned the marker one (1).n");
  2436.   printf("    - Otherwise, the vertex is assigned the marker zero (0).n");
  2437.   printf("n");
  2438.   printf(
  2439. "  If you want Triangle to determine for you which vertices and edges are onn"
  2440. );
  2441.   printf(
  2442. "  the boundary, assign them the boundary marker zero (or use no markers atn"
  2443. );
  2444.   printf(
  2445. "  all) in your input files.  In the output files, all boundary vertices,n");
  2446.   printf("  edges, and segments will be assigned the value one.nn");
  2447.   printf("Triangulation Iteration Numbers:nn");
  2448.   printf(
  2449. "  Because Triangle can read and refine its own triangulations, inputn");
  2450.   printf(
  2451. "  and output files have iteration numbers.  For instance, Triangle mightn");
  2452.   printf(
  2453. "  read the files mesh.3.node, mesh.3.ele, and mesh.3.poly, refine then");
  2454.   printf(
  2455. "  triangulation, and output the files mesh.4.node, mesh.4.ele, andn");
  2456.   printf("  mesh.4.poly.  Files with no iteration number are treated as ifn");
  2457.   printf(
  2458. "  their iteration number is zero; hence, Triangle might read the filen");
  2459.   printf(
  2460. "  points.node, triangulate it, and produce the files points.1.node andn");
  2461.   printf("  points.1.ele.nn");
  2462.   printf(
  2463. "  Iteration numbers allow you to create a sequence of successively finern");
  2464.   printf(
  2465. "  meshes suitable for multigrid methods.  They also allow you to produce an"
  2466. );
  2467.   printf(
  2468. "  sequence of meshes using error estimate-driven mesh refinement.n");
  2469.   printf("n");
  2470.   printf(
  2471. "  If you're not using refinement or quality meshing, and you don't liken");
  2472.   printf(
  2473. "  iteration numbers, use the -I switch to disable them.  This switch alson");
  2474.   printf(
  2475. "  disables output of .node and .poly files to prevent your input files fromn"
  2476. );
  2477.   printf(
  2478. "  being overwritten.  (If the input is a .poly file that contains its ownn");
  2479.   printf(
  2480. "  points, a .node file is written.  This can be quite convenient forn");
  2481.   printf("  computing CDTs or quality meshes.)nn");
  2482.   printf("Examples of How to Use Triangle:nn");
  2483.   printf(
  2484. "  `triangle dots' reads vertices from dots.node, and writes their Delaunayn"
  2485. );
  2486.   printf(
  2487. "  triangulation to dots.1.node and dots.1.ele.  (dots.1.node is identicaln");
  2488.   printf(
  2489. "  to dots.node.)  `triangle -I dots' writes the triangulation to dots.elen");
  2490.   printf(
  2491. "  instead.  (No additional .node file is needed, so none is written.)n");
  2492.   printf("n");
  2493.   printf(
  2494. "  `triangle -pe object.1' reads a PSLG from object.1.poly (and possiblyn");
  2495.   printf(
  2496. "  object.1.node, if the vertices are omitted from object.1.poly) and writesn"
  2497. );
  2498.   printf(
  2499. "  its constrained Delaunay triangulation to object.2.node and object.2.ele.n"
  2500. );
  2501.   printf(
  2502. "  The segments are copied to object.2.poly, and all edges are written ton");
  2503.   printf("  object.2.edge.nn");
  2504.   printf(
  2505. "  `triangle -pq31.5a.1 object' reads a PSLG from object.poly (and possiblyn"
  2506. );
  2507.   printf(
  2508. "  object.node), generates a mesh whose angles are all between 31.5 and 117n"
  2509. );
  2510.   printf(
  2511. "  degrees and whose triangles all have areas of 0.1 or less, and writes then"
  2512. );
  2513.   printf(
  2514. "  mesh to object.1.node and object.1.ele.  Each segment may be broken upn");
  2515.   printf("  into multiple subsegments; these are written to object.1.poly.n");
  2516.   printf("n");
  2517.   printf(
  2518. "  Here is a sample file `box.poly' describing a square with a square hole:n"
  2519. );
  2520.   printf("n");
  2521.   printf(
  2522. "    # A box with eight vertices in 2D, no attributes, one boundary marker.n"
  2523. );
  2524.   printf("    8 2 0 1n");
  2525.   printf("     # Outer box has these vertices:n");
  2526.   printf("     1   0 0   0n");
  2527.   printf("     2   0 3   0n");
  2528.   printf("     3   3 0   0n");
  2529.   printf("     4   3 3   33     # A special marker for this vertex.n");
  2530.   printf("     # Inner square has these vertices:n");
  2531.   printf("     5   1 1   0n");
  2532.   printf("     6   1 2   0n");
  2533.   printf("     7   2 1   0n");
  2534.   printf("     8   2 2   0n");
  2535.   printf("    # Five segments with boundary markers.n");
  2536.   printf("    5 1n");
  2537.   printf("     1   1 2   5      # Left side of outer box.n");
  2538.   printf("     # Square hole has these segments:n");
  2539.   printf("     2   5 7   0n");
  2540.   printf("     3   7 8   0n");
  2541.   printf("     4   8 6   10n");
  2542.   printf("     5   6 5   0n");
  2543.   printf("    # One hole in the middle of the inner square.n");
  2544.   printf("    1n");
  2545.   printf("     1   1.5 1.5n");
  2546.   printf("n");
  2547.   printf(
  2548. "  Note that some segments are missing from the outer square, so you mustn");
  2549.   printf(
  2550. "  use the `-c' switch.  After `triangle -pqc box.poly', here is the outputn"
  2551. );
  2552.   printf(
  2553. "  file `box.1.node', with twelve vertices.  The last four vertices weren");
  2554.   printf(
  2555. "  added to meet the angle constraint.  Vertices 1, 2, and 9 have markersn");
  2556.   printf(
  2557. "  from segment 1.  Vertices 6 and 8 have markers from segment 4.  All then");
  2558.   printf(
  2559. "  other vertices but 4 have been marked to indicate that they lie on an");
  2560.   printf("  boundary.nn");
  2561.   printf("    12  2  0  1n");
  2562.   printf("       1    0   0      5n");
  2563.   printf("       2    0   3      5n");
  2564.   printf("       3    3   0      1n");
  2565.   printf("       4    3   3     33n");
  2566.   printf("       5    1   1      1n");
  2567.   printf("       6    1   2     10n");
  2568.   printf("       7    2   1      1n");
  2569.   printf("       8    2   2     10n");
  2570.   printf("       9    0   1.5    5n");
  2571.   printf("      10    1.5   0    1n");
  2572.   printf("      11    3   1.5    1n");
  2573.   printf("      12    1.5   3    1n");
  2574.   printf("    # Generated by triangle -pqc box.polyn");
  2575.   printf("n");
  2576.   printf("  Here is the output file `box.1.ele', with twelve triangles.n");
  2577.   printf("n");
  2578.   printf("    12  3  0n");
  2579.   printf("       1     5   6   9n");
  2580.   printf("       2    10   3   7n");
  2581.   printf("       3     6   8  12n");
  2582.   printf("       4     9   1   5n");
  2583.   printf("       5     6   2   9n");
  2584.   printf("       6     7   3  11n");
  2585.   printf("       7    11   4   8n");
  2586.   printf("       8     7   5  10n");
  2587.   printf("       9    12   2   6n");
  2588.   printf("      10     8   7  11n");
  2589.   printf("      11     5   1  10n");
  2590.   printf("      12     8   4  12n");
  2591.   printf("    # Generated by triangle -pqc box.polynn");
  2592.   printf(
  2593. "  Here is the output file `box.1.poly'.  Note that segments have been addedn"
  2594. );
  2595.   printf(
  2596. "  to represent the convex hull, and some segments have been subdivided byn");
  2597.   printf(
  2598. "  newly added vertices.  Note also that <# of vertices> is set to zero ton");
  2599.   printf("  indicate that the vertices should be read from the .node file.n");
  2600.   printf("n");
  2601.   printf("    0  2  0  1n");
  2602.   printf("    12  1n");
  2603.   printf("       1     1   9     5n");
  2604.   printf("       2     5   7     1n");
  2605.   printf("       3     8   7     1n");
  2606.   printf("       4     6   8    10n");
  2607.   printf("       5     5   6     1n");
  2608.   printf("       6     3  10     1n");
  2609.   printf("       7     4  11     1n");
  2610.   printf("       8     2  12     1n");
  2611.   printf("       9     9   2     5n");
  2612.   printf("      10    10   1     1n");
  2613.   printf("      11    11   3     1n");
  2614.   printf("      12    12   4     1n");
  2615.   printf("    1n");
  2616.   printf("       1   1.5 1.5n");
  2617.   printf("    # Generated by triangle -pqc box.polyn");
  2618.   printf("n");
  2619.   printf("Refinement and Area Constraints:n");
  2620.   printf("n");
  2621.   printf(
  2622. "  The -r switch causes a mesh (.node and .ele files) to be read andn");
  2623.   printf(
  2624. "  refined.  If the -p switch is also used, a .poly file is read and used ton"
  2625. );
  2626.   printf(
  2627. "  specify edges that are constrained and cannot be eliminated (althoughn");
  2628.   printf(
  2629. "  they can be subdivided into smaller edges) by the refinement process.n");
  2630.   printf("n");
  2631.   printf(
  2632. "  When you refine a mesh, you generally want to impose tighter constraints.n"
  2633. );
  2634.   printf(
  2635. "  One way to accomplish this is to use -q with a larger angle, or -an");
  2636.   printf(
  2637. "  followed by a smaller area than you used to generate the mesh you aren");
  2638.   printf(
  2639. "  refining.  Another way to do this is to create an .area file, whichn");
  2640.   printf(
  2641. "  specifies a maximum area for each triangle, and use the -a switchn");
  2642.   printf(
  2643. "  (without a number following).  Each triangle's area constraint is appliedn"
  2644. );
  2645.   printf(
  2646. "  to that triangle.  Area constraints tend to diffuse as the mesh isn");
  2647.   printf(
  2648. "  refined, so if there are large variations in area constraint betweenn");
  2649.   printf(
  2650. "  adjacent triangles, you may not get the results you want.  In that case,n"
  2651. );
  2652.   printf(
  2653. "  consider instead using the -u switch and writing a C procedure thatn");
  2654.   printf("  determines which triangles are too large.nn");
  2655.   printf(
  2656. "  If you are refining a mesh composed of linear (three-node) elements, then"
  2657. );
  2658.   printf(
  2659. "  output mesh contains all the nodes present in the input mesh, in the samen"
  2660. );
  2661.   printf(
  2662. "  order, with new nodes added at the end of the .node file.  However, then");
  2663.   printf(
  2664. "  refinement is not hierarchical: there is no guarantee that each outputn");
  2665.   printf(
  2666. "  element is contained in a single input element.  Often, an output elementn"
  2667. );
  2668.   printf(
  2669. "  can overlap two or three input elements, and some input edges are notn");
  2670.   printf(
  2671. "  present in the output mesh.  Hence, a sequence of refined meshes forms an"
  2672. );
  2673.   printf(
  2674. "  hierarchy of nodes, but not a hierarchy of elements.  If you refine an");
  2675.   printf(
  2676. "  mesh of higher-order elements, the hierarchical property applies only ton"
  2677. );
  2678.   printf(
  2679. "  the nodes at the corners of an element; the midpoint nodes on each edgen");
  2680.   printf("  are discarded before the mesh is refined.nn");
  2681.   printf(
  2682. "  Maximum area constraints in .poly files operate differently from those inn"
  2683. );
  2684.   printf(
  2685. "  .area files.  A maximum area in a .poly file applies to the wholen");
  2686.   printf(
  2687. "  (segment-bounded) region in which a point falls, whereas a maximum arean");
  2688.   printf(
  2689. "  in an .area file applies to only one triangle.  Area constraints in .polyn"
  2690. );
  2691.   printf(
  2692. "  files are used only when a mesh is first generated, whereas arean");
  2693.   printf(
  2694. "  constraints in .area files are used only to refine an existing mesh, andn"
  2695. );
  2696.   printf(
  2697. "  are typically based on a posteriori error estimates resulting from an");
  2698.   printf("  finite element simulation on that mesh.nn");
  2699.   printf(
  2700. "  `triangle -rq25 object.1' reads object.1.node and object.1.ele, thenn");
  2701.   printf(
  2702. "  refines the triangulation to enforce a 25 degree minimum angle, and thenn"
  2703. );
  2704.   printf(
  2705. "  writes the refined triangulation to object.2.node and object.2.ele.n");
  2706.   printf("n");
  2707.   printf(
  2708. "  `triangle -rpaa6.2 z.3' reads z.3.node, z.3.ele, z.3.poly, and z.3.area.n"
  2709. );
  2710.   printf(
  2711. "  After reconstructing the mesh and its subsegments, Triangle refines then");
  2712.   printf(
  2713. "  mesh so that no triangle has area greater than 6.2, and furthermore then");
  2714.   printf(
  2715. "  triangles satisfy the maximum area constraints in z.3.area.  No anglen");
  2716.   printf(
  2717. "  bound is imposed at all.  The output is written to z.4.node, z.4.ele, andn"
  2718. );
  2719.   printf("  z.4.poly.nn");
  2720.   printf(
  2721. "  The sequence `triangle -qa1 x', `triangle -rqa.3 x.1', `triangle -rqa.1n");
  2722.   printf(
  2723. "  x.2' creates a sequence of successively finer meshes x.1, x.2, and x.3,n");
  2724.   printf("  suitable for multigrid.nn");
  2725.   printf("Convex Hulls and Mesh Boundaries:nn");
  2726.   printf(
  2727. "  If the input is a vertex set (not a PSLG), Triangle produces its convexn");
  2728.   printf(
  2729. "  hull as a by-product in the output .poly file if you use the -c switch.n");
  2730.   printf(
  2731. "  There are faster algorithms for finding a two-dimensional convex hulln");
  2732.   printf("  than triangulation, of course, but this one comes for free.nn");
  2733.   printf(
  2734. "  If the input is an unconstrained mesh (you are using the -r switch butn");
  2735.   printf(
  2736. "  not the -p switch), Triangle produces a list of its boundary edgesn");
  2737.   printf(
  2738. "  (including hole boundaries) as a by-product when you use the -c switch.n");
  2739.   printf(
  2740. "  If you also use the -p switch, the output .poly file contains all then");
  2741.   printf("  segments from the input .poly file as well.nn");
  2742.   printf("Voronoi Diagrams:nn");
  2743.   printf(
  2744. "  The -v switch produces a Voronoi diagram, in files suffixed .v.node andn");
  2745.   printf(
  2746. "  .v.edge.  For example, `triangle -v points' reads points.node, producesn");
  2747.   printf(
  2748. "  its Delaunay triangulation in points.1.node and points.1.ele, andn");
  2749.   printf(
  2750. "  produces its Voronoi diagram in points.1.v.node and points.1.v.edge.  Then"
  2751. );
  2752.   printf(
  2753. "  .v.node file contains a list of all Voronoi vertices, and the .v.edgen");
  2754.   printf(
  2755. "  file contains a list of all Voronoi edges, some of which may be infiniten"
  2756. );
  2757.   printf(
  2758. "  rays.  (The choice of filenames makes it easy to run the set of Voronoin");
  2759.   printf("  vertices through Triangle, if so desired.)nn");
  2760.   printf(
  2761. "  This implementation does not use exact arithmetic to compute the Voronoin"
  2762. );
  2763.   printf(
  2764. "  vertices, and does not check whether neighboring vertices are identical.n"
  2765. );
  2766.   printf(
  2767. "  Be forewarned that if the Delaunay triangulation is degenerate orn");
  2768.   printf(
  2769. "  near-degenerate, the Voronoi diagram may have duplicate vertices orn");
  2770.   printf("  crossing edges.nn");
  2771.   printf(
  2772. "  The result is a valid Voronoi diagram only if Triangle's output is a truen"
  2773. );
  2774.   printf(
  2775. "  Delaunay triangulation.  The Voronoi output is usually meaningless (andn");
  2776.   printf(
  2777. "  may contain crossing edges and other pathology) if the output is a CDT orn"
  2778. );
  2779.   printf(
  2780. "  CCDT, or if it has holes or concavities.  If the triangulated domain isn");
  2781.   printf(
  2782. "  convex and has no holes, you can use -D switch to force Triangle ton");
  2783.   printf(
  2784. "  construct a conforming Delaunay triangulation instead of a CCDT, so then");
  2785.   printf("  Voronoi diagram will be valid.nn");
  2786.   printf("Mesh Topology:nn");
  2787.   printf(
  2788. "  You may wish to know which triangles are adjacent to a certain Delaunayn");
  2789.   printf(
  2790. "  edge in an .edge file, which Voronoi cells are adjacent to a certainn");
  2791.   printf(
  2792. "  Voronoi edge in a .v.edge file, or which Voronoi cells are adjacent ton");
  2793.   printf(
  2794. "  each other.  All of this information can be found by cross-referencingn");
  2795.   printf(
  2796. "  output files with the recollection that the Delaunay triangulation andn");
  2797.   printf("  the Voronoi diagram are planar duals.nn");
  2798.   printf(
  2799. "  Specifically, edge i of an .edge file is the dual of Voronoi edge i ofn");
  2800.   printf(
  2801. "  the corresponding .v.edge file, and is rotated 90 degrees counterclock-n");
  2802.   printf(
  2803. "  wise from the Voronoi edge.  Triangle j of an .ele file is the dual ofn");
  2804.   printf(
  2805. "  vertex j of the corresponding .v.node file.  Voronoi cell k is the dualn");
  2806.   printf("  of vertex k of the corresponding .node file.nn");
  2807.   printf(
  2808. "  Hence, to find the triangles adjacent to a Delaunay edge, look at then");
  2809.   printf(
  2810. "  vertices of the corresponding Voronoi edge.  If the endpoints of an");
  2811.   printf(
  2812. "  Voronoi edge are Voronoi vertices 2 and 6 respectively, then triangles 2n"
  2813. );
  2814.   printf(
  2815. "  and 6 adjoin the left and right sides of the corresponding Delaunay edge,n"
  2816. );
  2817.   printf(
  2818. "  respectively.  To find the Voronoi cells adjacent to a Voronoi edge, lookn"
  2819. );
  2820.   printf(
  2821. "  at the endpoints of the corresponding Delaunay edge.  If the endpoints ofn"
  2822. );
  2823.   printf(
  2824. "  a Delaunay edge are input vertices 7 and 12, then Voronoi cells 7 and 12n"
  2825. );
  2826.   printf(
  2827. "  adjoin the right and left sides of the corresponding Voronoi edge,n");
  2828.   printf(
  2829. "  respectively.  To find which Voronoi cells are adjacent to each other,n");
  2830.   printf("  just read the list of Delaunay edges.nn");
  2831.   printf(
  2832. "  Triangle does not write a list of the edges adjoining each Voronoi cell,n"
  2833. );
  2834.   printf(
  2835. "  but you can reconstructed it straightforwardly.  For instance, to findn");
  2836.   printf(
  2837. "  all the edges of Voronoi cell 1, search the output .edge file for everyn");
  2838.   printf(
  2839. "  edge that has input vertex 1 as an endpoint.  The corresponding dualn");
  2840.   printf(
  2841. "  edges in the output .v.edge file form the boundary of Voronoi cell 1.n");
  2842.   printf("n");
  2843.   printf(
  2844. "  For each Voronoi vertex, the .neigh file gives a list of the threen");
  2845.   printf(
  2846. "  Voronoi vertices attached to it.  You might find this more convenientn");
  2847.   printf("  than the .v.edge file.nn");
  2848.   printf("Quadratic Elements:nn");
  2849.   printf(
  2850. "  Triangle generates meshes with subparametric quadratic elements if then");
  2851.   printf(
  2852. "  -o2 switch is specified.  Quadratic elements have six nodes per element,n"
  2853. );
  2854.   printf(
  2855. "  rather than three.  `Subparametric' means that the edges of the trianglesn"
  2856. );
  2857.   printf(
  2858. "  are always straight, so that subparametric quadratic elements aren");
  2859.   printf(
  2860. "  geometrically identical to linear elements, even though they can be usedn"
  2861. );
  2862.   printf(
  2863. "  with quadratic interpolating functions.  The three extra nodes of ann");
  2864.   printf(
  2865. "  element fall at the midpoints of the three edges, with the fourth, fifth,n"
  2866. );
  2867.   printf(
  2868. "  and sixth nodes appearing opposite the first, second, and third cornersn");
  2869.   printf("  respectively.nn");
  2870.   printf("Domains with Small Angles:nn");
  2871.   printf(
  2872. "  If two input segments adjoin each other at a small angle, clearly the -qn"
  2873. );
  2874.   printf(
  2875. "  switch cannot remove the small angle.  Moreover, Triangle may have non");
  2876.   printf(
  2877. "  choice but to generate additional triangles whose smallest angles aren");
  2878.   printf(
  2879. "  smaller than the specified bound.  However, these triangles only appearn");
  2880.   printf(
  2881. "  between input segments separated by small angles.  Moreover, if youn");
  2882.   printf(
  2883. "  request a minimum angle of theta degrees, Triangle will generally producen"
  2884. );
  2885.   printf(
  2886. "  no angle larger than 180 - 2 theta, even if it is forced to compromise onn"
  2887. );
  2888.   printf("  the minimum angle.nn");
  2889.   printf("Statistics:nn");
  2890.   printf(
  2891. "  After generating a mesh, Triangle prints a count of entities in then");
  2892.   printf(
  2893. "  output mesh, including the number of vertices, triangles, edges, exteriorn"
  2894. );
  2895.   printf(
  2896. "  boundary edges (i.e. subsegments on the boundary of the triangulation,n");
  2897.   printf(
  2898. "  including hole boundaries), interior boundary edges (i.e. subsegments ofn"
  2899. );
  2900.   printf(
  2901. "  input segments not on the boundary), and total subsegments.  If you'ven");
  2902.   printf(
  2903. "  forgotten the statistics for an existing mesh, run Triangle on that meshn"
  2904. );
  2905.   printf(
  2906. "  with the -rNEP switches to read the mesh and print the statistics withoutn"
  2907. );
  2908.   printf(
  2909. "  writing any files.  Use -rpNEP if you've got a .poly file for the mesh.n");
  2910.   printf("n");
  2911.   printf(
  2912. "  The -V switch produces extended statistics, including a rough estimaten");
  2913.   printf(
  2914. "  of memory use, the number of calls to geometric predicates, andn");
  2915.   printf(
  2916. "  histograms of the angles and the aspect ratios of the triangles in then");
  2917.   printf("  mesh.nn");
  2918.   printf("Exact Arithmetic:nn");
  2919.   printf(
  2920. "  Triangle uses adaptive exact arithmetic to perform what computationaln");
  2921.   printf(
  2922. "  geometers call the `orientation' and `incircle' tests.  If the floating-n"
  2923. );
  2924.   printf(
  2925. "  point arithmetic of your machine conforms to the IEEE 754 standard (asn");
  2926.   printf(
  2927. "  most workstations do), and does not use extended precision internaln");
  2928.   printf(
  2929. "  floating-point registers, then your output is guaranteed to be ann");
  2930.   printf(
  2931. "  absolutely true Delaunay or constrained Delaunay triangulation, roundoffn"
  2932. );
  2933.   printf(
  2934. "  error notwithstanding.  The word `adaptive' implies that these arithmeticn"
  2935. );
  2936.   printf(
  2937. "  routines compute the result only to the precision necessary to guaranteen"
  2938. );
  2939.   printf(
  2940. "  correctness, so they are usually nearly as fast as their approximaten");
  2941.   printf("  counterparts.nn");
  2942.   printf(
  2943. "  May CPUs, including Intel x86 processors, have extended precisionn");
  2944.   printf(
  2945. "  floating-point registers.  These must be reconfigured so their precisionn"
  2946. );
  2947.   printf(
  2948. "  is reduced to memory precision.  Triangle does this if it is compiledn");
  2949.   printf("  correctly.  See the makefile for details.nn");
  2950.   printf(
  2951. "  The exact tests can be disabled with the -X switch.  On most inputs, thisn"
  2952. );
  2953.   printf(
  2954. "  switch reduces the computation time by about eight percent--it's notn");
  2955.   printf(
  2956. "  worth the risk.  There are rare difficult inputs (having many collinearn");
  2957.   printf(
  2958. "  and cocircular vertices), however, for which the difference in speedn");
  2959.   printf(
  2960. "  could be a factor of two.  Be forewarned that these are precisely then");
  2961.   printf(
  2962. "  inputs most likely to cause errors if you use the -X switch.  Hence, then"
  2963. );
  2964.   printf("  -X switch is not recommended.nn");
  2965.   printf(
  2966. "  Unfortunately, the exact tests don't solve every numerical problem.n");
  2967.   printf(
  2968. "  Exact arithmetic is not used to compute the positions of new vertices,n");
  2969.   printf(
  2970. "  because the bit complexity of vertex coordinates would grow withoutn");
  2971.   printf(
  2972. "  bound.  Hence, segment intersections aren't computed exactly; in veryn");
  2973.   printf(
  2974. "  unusual cases, roundoff error in computing an intersection point mightn");
  2975.   printf(
  2976. "  actually lead to an inverted triangle and an invalid triangulation.n");
  2977.   printf(
  2978. "  (This is one reason to specify your own intersection points in your .polyn"
  2979. );
  2980.   printf(
  2981. "  files.)  Similarly, exact arithmetic is not used to compute the verticesn"
  2982. );
  2983.   printf("  of the Voronoi diagram.nn");
  2984.   printf(
  2985. "  Another pair of problems not solved by the exact arithmetic routines isn");
  2986.   printf(
  2987. "  underflow and overflow.  If Triangle is compiled for double precisionn");
  2988.   printf(
  2989. "  arithmetic, I believe that Triangle's geometric predicates work correctlyn"
  2990. );
  2991.   printf(
  2992. "  if the exponent of every input coordinate falls in the range [-148, 201].n"
  2993. );
  2994.   printf(
  2995. "  Underflow can silently prevent the orientation and incircle tests fromn");
  2996.   printf(
  2997. "  being performed exactly, while overflow typically causes a floatingn");
  2998.   printf("  exception.nn");
  2999.   printf("Calling Triangle from Another Program:nn");
  3000.   printf("  Read the file triangle.h for details.nn");
  3001.   printf("Troubleshooting:nn");
  3002.   printf("  Please read this section before mailing me bugs.nn");
  3003.   printf("  `My output mesh has no triangles!'nn");
  3004.   printf(
  3005. "    If you're using a PSLG, you've probably failed to specify a proper setn"
  3006. );
  3007.   printf(
  3008. "    of bounding segments, or forgotten to use the -c switch.  Or you mayn");
  3009.   printf(
  3010. "    have placed a hole badly, thereby eating all your triangles.  To testn");
  3011.   printf("    these possibilities, try again with the -c and -O switches.n");
  3012.   printf(
  3013. "    Alternatively, all your input vertices may be collinear, in which casen"
  3014. );
  3015.   printf("    you can hardly expect to triangulate them.nn");
  3016.   printf("  `Triangle doesn't terminate, or just crashes.'nn");
  3017.   printf(
  3018. "    Bad things can happen when triangles get so small that the distancen");
  3019.   printf(
  3020. "    between their vertices isn't much larger than the precision of yourn");
  3021.   printf(
  3022. "    machine's arithmetic.  If you've compiled Triangle for single-precisionn"
  3023. );
  3024.   printf(
  3025. "    arithmetic, you might do better by recompiling it for double-precision.n"
  3026. );
  3027.   printf(
  3028. "    Then again, you might just have to settle for more lenient constraintsn"
  3029. );
  3030.   printf(
  3031. "    on the minimum angle and the maximum area than you had planned.n");
  3032.   printf("n");
  3033.   printf(
  3034. "    You can minimize precision problems by ensuring that the origin liesn");
  3035.   printf(
  3036. "    inside your vertex set, or even inside the densest part of yourn");
  3037.   printf(
  3038. "    mesh.  If you're triangulating an object whose x-coordinates all falln");
  3039.   printf(
  3040. "    between 6247133 and 6247134, you're not leaving much floating-pointn");
  3041.   printf("    precision for Triangle to work with.nn");
  3042.   printf(
  3043. "    Precision problems can occur covertly if the input PSLG contains twon");
  3044.   printf(
  3045. "    segments that meet (or intersect) at an extremely small angle, or ifn");
  3046.   printf(
  3047. "    such an angle is introduced by the -c switch.  If you don't realizen");
  3048.   printf(
  3049. "    that a tiny angle is being formed, you might never discover whyn");
  3050.   printf(
  3051. "    Triangle is crashing.  To check for this possibility, use the -S switchn"
  3052. );
  3053.   printf(
  3054. "    (with an appropriate limit on the number of Steiner points, found byn");
  3055.   printf(
  3056. "    trial-and-error) to stop Triangle early, and view the output .poly filen"
  3057. );
  3058.   printf(
  3059. "    with Show Me (described below).  Look carefully for regions where densen"
  3060. );
  3061.   printf(
  3062. "    clusters of vertices are forming and for small angles between segments.n"
  3063. );
  3064.   printf(
  3065. "    Zoom in closely, as such segments might look like a single segment fromn"
  3066. );
  3067.   printf("    a distance.nn");
  3068.   printf(
  3069. "    If some of the input values are too large, Triangle may suffer an");
  3070.   printf(
  3071. "    floating exception due to overflow when attempting to perform ann");
  3072.   printf(
  3073. "    orientation or incircle test.  (Read the section on exact arithmeticn");
  3074.   printf(
  3075. "    above.)  Again, I recommend compiling Triangle for double (rathern");
  3076.   printf("    than single) precision arithmetic.nn");
  3077.   printf(
  3078. "    Unexpected problems can arise if you use quality meshing (-q, -a, orn");
  3079.   printf(
  3080. "    -u) with an input that is not segment-bounded--that is, if your inputn");
  3081.   printf(
  3082. "    is a vertex set, or you're using the -c switch.  If the convex hull ofn"
  3083. );
  3084.   printf(
  3085. "    your input vertices has collinear vertices on its boundary, an inputn");
  3086.   printf(
  3087. "    vertex that you think lies on the convex hull might actually lie justn");
  3088.   printf(
  3089. "    inside the convex hull.  If so, the vertex and the nearby convex hulln");
  3090.   printf(
  3091. "    edge form an extremely thin triangle.  When Triangle tries to refinen");
  3092.   printf(
  3093. "    the mesh to enforce angle and area constraints, Triangle might generaten"
  3094. );
  3095.   printf(
  3096. "    extremely tiny triangles, or it might fail because of insufficientn");
  3097.   printf("    floating-point precision.nn");
  3098.   printf(
  3099. "  `The numbering of the output vertices doesn't match the input vertices.'n"
  3100. );
  3101.   printf("n");
  3102.   printf(
  3103. "    You may have had duplicate input vertices, or you may have eaten somen");
  3104.   printf(
  3105. "    of your input vertices with a hole, or by placing them outside the arean"
  3106. );
  3107.   printf(
  3108. "    enclosed by segments.  In any case, you can solve the problem by notn");
  3109.   printf("    using the -j switch.nn");
  3110.   printf(
  3111. "  `Triangle executes without incident, but when I look at the resultingn");
  3112.   printf(
  3113. "  mesh, it has overlapping triangles or other geometric inconsistencies.'n");
  3114.   printf("n");
  3115.   printf(
  3116. "    If you select the -X switch, Triangle occasionally makes mistakes duen");
  3117.   printf(
  3118. "    to floating-point roundoff error.  Although these errors are rare,n");
  3119.   printf(
  3120. "    don't use the -X switch.  If you still have problems, please report then"
  3121. );
  3122.   printf("    bug.nn");
  3123.   printf(
  3124. "  `Triangle executes without incident, but when I look at the resultingn");
  3125.   printf("  Voronoi diagram, it has overlapping edges or other geometricn");
  3126.   printf("  inconsistencies.'n");
  3127.   printf("n");
  3128.   printf(
  3129. "    If your input is a PSLG (-p), you can only expect a meaningful Voronoin"
  3130. );
  3131.   printf(
  3132. "    diagram if the domain you are triangulating is convex and free ofn");
  3133.   printf(
  3134. "    holes, and you use the -D switch to construct a conforming Delaunayn");
  3135.   printf("    triangulation (instead of a CDT or CCDT).nn");
  3136.   printf(
  3137. "  Strange things can happen if you've taken liberties with your PSLG.  Don");
  3138.   printf(
  3139. "  you have a vertex lying in the middle of a segment?  Triangle sometimesn");
  3140.   printf(
  3141. "  copes poorly with that sort of thing.  Do you want to lay out a collinearn"
  3142. );
  3143.   printf(
  3144. "  row of evenly spaced, segment-connected vertices?  Have you simplyn");
  3145.   printf(
  3146. "  defined one long segment connecting the leftmost vertex to the rightmostn"
  3147. );
  3148.   printf(
  3149. "  vertex, and a bunch of vertices lying along it?  This method occasionallyn"
  3150. );
  3151.   printf(
  3152. "  works, especially with horizontal and vertical lines, but often itn");
  3153.   printf(
  3154. "  doesn't, and you'll have to connect each adjacent pair of vertices with an"
  3155. );
  3156.   printf("  separate segment.  If you don't like it, tough.nn");
  3157.   printf(
  3158. "  Furthermore, if you have segments that intersect other than at theirn");
  3159.   printf(
  3160. "  endpoints, try not to let the intersections fall extremely close to PSLGn"
  3161. );
  3162.   printf("  vertices or each other.nn");
  3163.   printf(
  3164. "  If you have problems refining a triangulation not produced by Triangle:n");
  3165.   printf(
  3166. "  Are you sure the triangulation is geometrically valid?  Is it formattedn");
  3167.   printf(
  3168. "  correctly for Triangle?  Are the triangles all listed so the first threen"
  3169. );
  3170.   printf(
  3171. "  vertices are their corners in counterclockwise order?  Are all of then");
  3172.   printf(
  3173. "  triangles constrained Delaunay?  Triangle's Delaunay refinement algorithmn"
  3174. );
  3175.   printf("  assumes that it starts with a CDT.nn");
  3176.   printf("Show Me:nn");
  3177.   printf(
  3178. "  Triangle comes with a separate program named `Show Me', whose primaryn");
  3179.   printf(
  3180. "  purpose is to draw meshes on your screen or in PostScript.  Its secondaryn"
  3181. );
  3182.   printf(
  3183. "  purpose is to check the validity of your input files, and do so moren");
  3184.   printf(
  3185. "  thoroughly than Triangle does.  Unlike Triangle, Show Me requires thatn");
  3186.   printf(
  3187. "  you have the X Windows system.  Sorry, Microsoft Windows users.n");
  3188.   printf("n");
  3189.   printf("Triangle on the Web:n");
  3190.   printf("n");
  3191.   printf("  To see an illustrated version of these instructions, check outn");
  3192.   printf("n");
  3193.   printf("    http://www.cs.cmu.edu/~quake/triangle.htmln");
  3194.   printf("n");
  3195.   printf("A Brief Plea:n");
  3196.   printf("n");
  3197.   printf(
  3198. "  If you use Triangle, and especially if you use it to accomplish realn");
  3199.   printf(
  3200. "  work, I would like very much to hear from you.  A short letter or emailn");
  3201.   printf(
  3202. "  (to jrs@cs.berkeley.edu) describing how you use Triangle will mean a lotn"