/*****************************************************************************/ /* */ /* 888888888 ,o, / 888 */ /* 888 88o88o " o8888o 88o8888o o88888o 888 o88888o */ /* 888 888 888 88b 888 888 888 888 888 d888 88b */ /* 888 888 888 o88^o888 888 888 "88888" 888 8888oo888 */ /* 888 888 888 C888 888 888 888 / 888 q888 */ /* 888 888 888 "88o^888 888 888 Cb 888 "88oooo" */ /* "8oo8D */ /* */ /* A Two-Dimensional Quality Mesh Generator and Delaunay Triangulator. */ /* (triangle.c) */ /* */ /* Version 1.4 */ /* November 1, 2002 */ /* */ /* Copyright 1993, 1995, 1997, 1998, 2002 */ /* Jonathan Richard Shewchuk */ /* 2360 Woolsey #H */ /* Berkeley, California 94705-1927 */ /* jrs@cs.berkeley.edu */ /* */ /* This program may be freely redistributed under the condition that the */ /* copyright notices (including this entire header and the copyright */ /* notice printed when the `-h' switch is selected) are not removed, and */ /* no compensation is received. Private, research, and institutional */ /* use is free. You may distribute modified versions of this code UNDER */ /* THE CONDITION THAT THIS CODE AND ANY MODIFICATIONS MADE TO IT IN THE */ /* SAME FILE REMAIN UNDER COPYRIGHT OF THE ORIGINAL AUTHOR, BOTH SOURCE */ /* AND OBJECT CODE ARE MADE FREELY AVAILABLE WITHOUT CHARGE, AND CLEAR */ /* NOTICE IS GIVEN OF THE MODIFICATIONS. Distribution of this code as */ /* part of a commercial system is permissible ONLY BY DIRECT ARRANGEMENT */ /* WITH THE AUTHOR. (If you are not directly supplying this code to a */ /* customer, and you are instead telling them how they can obtain it for */ /* free, then you are not required to make any arrangement with me.) */ /* */ /* Hypertext instructions for Triangle are available on the Web at */ /* */ /* http://www.cs.cmu.edu/~quake/triangle.html */ /* */ /* Some of the references listed below are marked with an asterisk. [*] */ /* These references are available for downloading from the Web page */ /* */ /* http://www.cs.cmu.edu/~quake/triangle.research.html */ /* */ /* Three papers discussing aspects of Triangle are available. A short */ /* overview appears in "Triangle: Engineering a 2D Quality Mesh */ /* Generator and Delaunay Triangulator," in Applied Computational */ /* Geometry: Towards Geometric Engineering, Ming C. Lin and Dinesh */ /* Manocha, editors, Lecture Notes in Computer Science volume 1148, */ /* pages 203-222, Springer-Verlag, Berlin, May 1996 (from the First ACM */ /* Workshop on Applied Computational Geometry). [*] */ /* */ /* The algorithms are discussed in the greatest detail in "Delaunay */ /* Refinement Algorithms for Triangular Mesh Generation," Computational */ /* Geometry: Theory and Applications 22(1-3):21-74, May 2002. [*] */ /* */ /* More detail about the data structures may be found in my dissertation: */ /* "Delaunay Refinement Mesh Generation," Ph.D. thesis, Technical Report */ /* CMU-CS-97-137, School of Computer Science, Carnegie Mellon University, */ /* Pittsburgh, Pennsylvania, 18 May 1997. [*] */ /* */ /* Triangle was created as part of the Archimedes project in the School of */ /* Computer Science at Carnegie Mellon University. Archimedes is a */ /* system for compiling parallel finite element solvers. For further */ /* information, see Hesheng Bao, Jacobo Bielak, Omar Ghattas, Loukas F. */ /* Kallivokas, David R. O'Hallaron, Jonathan R. Shewchuk, and Jifeng Xu, */ /* "Large-scale Simulation of Elastic Wave Propagation in Heterogeneous */ /* Media on Parallel Computers," Computer Methods in Applied Mechanics */ /* and Engineering 152(1-2):85-102, 22 January 1998. */ /* */ /* Triangle's Delaunay refinement algorithm for quality mesh generation is */ /* a hybrid of one due to Jim Ruppert, "A Delaunay Refinement Algorithm */ /* for Quality 2-Dimensional Mesh Generation," Journal of Algorithms */ /* 18(3):548-585, May 1995 [*], and one due to L. Paul Chew, "Guaranteed- */ /* Quality Mesh Generation for Curved Surfaces," Proceedings of the Ninth */ /* Annual Symposium on Computational Geometry (San Diego, California), */ /* pages 274-280, Association for Computing Machinery, May 1993. */ /* */ /* The Delaunay refinement algorithm has been modified so that it */ /* consistently meshes domains with small input angles, as described in */ /* my lengthy journal article listed above, or in abbreviated form in */ /* Jonathan Richard Shewchuk, "Mesh Generation for Domains with Small */ /* Angles," Proceedings of the Sixteenth Annual Symposium on */ /* Computational Geometry (Hong Kong), pages 1-10, Association for */ /* Computing Machinery, June 2000. [*] */ /* */ /* My implementation of the divide-and-conquer and incremental Delaunay */ /* triangulation algorithms follows closely the presentation of Guibas */ /* and Stolfi, even though I use a triangle-based data structure instead */ /* of their quad-edge data structure. (In fact, I originally implemented */ /* Triangle using the quad-edge data structure, but the switch to a */ /* triangle-based data structure sped Triangle by a factor of two.) The */ /* mesh manipulation primitives and the two aforementioned Delaunay */ /* triangulation algorithms are described by Leonidas J. Guibas and Jorge */ /* Stolfi, "Primitives for the Manipulation of General Subdivisions and */ /* the Computation of Voronoi Diagrams," ACM Transactions on Graphics */ /* 4(2):74-123, April 1985. */ /* */ /* Their O(n log n) divide-and-conquer algorithm is adapted from Der-Tsai */ /* Lee and Bruce J. Schachter, "Two Algorithms for Constructing the */ /* Delaunay Triangulation," International Journal of Computer and */ /* Information Science 9(3):219-242, 1980. Triangle's improvement of the */ /* divide-and-conquer algorithm by alternating between vertical and */ /* horizontal cuts was introduced by Rex A. Dwyer, "A Faster Divide-and- */ /* Conquer Algorithm for Constructing Delaunay Triangulations," */ /* Algorithmica 2(2):137-151, 1987. */ /* */ /* The incremental insertion algorithm was first proposed by C. L. Lawson, */ /* "Software for C1 Surface Interpolation," in Mathematical Software III, */ /* John R. Rice, editor, Academic Press, New York, pp. 161-194, 1977. */ /* For point location, I use the algorithm of Ernst P. Mucke, Isaac */ /* Saias, and Binhai Zhu, "Fast Randomized Point Location Without */ /* Preprocessing in Two- and Three-Dimensional Delaunay Triangulations," */ /* Proceedings of the Twelfth Annual Symposium on Computational Geometry, */ /* ACM, May 1996. [*] If I were to randomize the order of vertex */ /* insertion (I currently don't bother), their result combined with the */ /* result of Kenneth L. Clarkson and Peter W. Shor, "Applications of */ /* Random Sampling in Computational Geometry II," Discrete & */ /* Computational Geometry 4(1):387-421, 1989, would yield an expected */ /* O(n^{4/3}) bound on running time. */ /* */ /* The O(n log n) sweepline Delaunay triangulation algorithm is taken from */ /* Steven Fortune, "A Sweepline Algorithm for Voronoi Diagrams", */ /* Algorithmica 2(2):153-174, 1987. A random sample of edges on the */ /* boundary of the triangulation are maintained in a splay tree for the */ /* purpose of point location. Splay trees are described by Daniel */ /* Dominic Sleator and Robert Endre Tarjan, "Self-Adjusting Binary Search */ /* Trees," Journal of the ACM 32(3):652-686, July 1985. */ /* */ /* The algorithms for exact computation of the signs of determinants are */ /* described in Jonathan Richard Shewchuk, "Adaptive Precision Floating- */ /* Point Arithmetic and Fast Robust Geometric Predicates," Discrete & */ /* Computational Geometry 18(3):305-363, October 1997. (Also available */ /* as Technical Report CMU-CS-96-140, School of Computer Science, */ /* Carnegie Mellon University, Pittsburgh, Pennsylvania, May 1996.) [*] */ /* An abbreviated version appears as Jonathan Richard Shewchuk, "Robust */ /* Adaptive Floating-Point Geometric Predicates," Proceedings of the */ /* Twelfth Annual Symposium on Computational Geometry, ACM, May 1996. [*] */ /* Many of the ideas for my exact arithmetic routines originate with */ /* Douglas M. Priest, "Algorithms for Arbitrary Precision Floating Point */ /* Arithmetic," Tenth Symposium on Computer Arithmetic, pp. 132-143, IEEE */ /* Computer Society Press, 1991. [*] Many of the ideas for the correct */ /* evaluation of the signs of determinants are taken from Steven Fortune */ /* and Christopher J. Van Wyk, "Efficient Exact Arithmetic for Computa- */ /* tional Geometry," Proceedings of the Ninth Annual Symposium on */ /* Computational Geometry, ACM, pp. 163-172, May 1993, and from Steven */ /* Fortune, "Numerical Stability of Algorithms for 2D Delaunay Triangu- */ /* lations," International Journal of Computational Geometry & Applica- */ /* tions 5(1-2):193-213, March-June 1995. */ /* */ /* For definitions of and results involving Delaunay triangulations, */ /* constrained and conforming versions thereof, and other aspects of */ /* triangular mesh generation, see the excellent survey by Marshall Bern */ /* and David Eppstein, "Mesh Generation and Optimal Triangulation," in */ /* Computing and Euclidean Geometry, Ding-Zhu Du and Frank Hwang, */ /* editors, World Scientific, Singapore, pp. 23-90, 1992. [*] */ /* */ /* The time for incrementally adding PSLG (planar straight line graph) */ /* segments to create a constrained Delaunay triangulation is probably */ /* O(t^2) per segment in the worst case and O(t) per segment in the */ /* common case, where t is the number of triangles that intersect the */ /* segment before it is inserted. This doesn't count point location, */ /* which can be much more expensive. I could improve this to O(d log d) */ /* time, but d is usually quite small, so it's not worth the bother. */ /* (This note does not apply to conforming Delaunay triangulations, for */ /* which a different method is used to insert segments.) */ /* */ /* The time for adding segments to a conforming Delaunay triangulation is */ /* not clear, but does not depend upon t alone. In some cases, very */ /* small features (like a vertex lying next to a segment) can cause a */ /* single segment to be split an arbitrary number of times. Of course, */ /* floating-point precision is a practical barrier to how much this can */ /* happen. */ /* */ /* The time for deleting a vertex from a Delaunay triangulation is O(d^2) */ /* in the worst case and O(d) in the common case, where d is the degree */ /* of the vertex being deleted. I could improve this to O(d log d) time, */ /* but d is usually quite small, so it's not worth the bother. */ /* */ /* Ruppert's Delaunay refinement algorithm typically generates triangles */ /* at a linear rate (constant time per triangle) after the initial */ /* triangulation is formed. There may be pathological cases where */ /* quadratic time is required, but these never arise in practice. */ /* */ /* The geometric predicates (circumcenter calculations, segment */ /* intersection formulae, etc.) appear in my "Lecture Notes on Geometric */ /* Robustness" at http://www.cs.berkeley.edu/~jrs/mesh.html . */ /* */ /* If you make any improvements to this code, please please please let me */ /* know, so that I may obtain the improvements. Even if you don't change */ /* the code, I'd still love to hear what it's being used for. */ /* */ /* Disclaimer: Neither I nor Carnegie Mellon warrant this code in any way */ /* whatsoever. This code is provided "as-is". Use at your own risk. */ /* */ /*****************************************************************************/ /* For single precision (which will save some memory and reduce paging), */ /* define the symbol SINGLE by using the -DSINGLE compiler switch or by */ /* writing "#define SINGLE" below. */ /* */ /* For double precision (which will allow you to refine meshes to a smaller */ /* edge length), leave SINGLE undefined. */ /* */ /* Double precision uses more memory, but improves the resolution of the */ /* meshes you can generate with Triangle. It also reduces the likelihood */ /* of a floating exception due to overflow. Finally, it is much faster */ /* than single precision on 64-bit architectures like the DEC Alpha. I */ /* recommend double precision unless you want to generate a mesh for which */ /* you do not have enough memory. */ /* #define SINGLE */ #ifdef SINGLE #define REAL float #else /* not SINGLE */ #define REAL double #endif /* not SINGLE */ /* If yours is not a Unix system, define the NO_TIMER compiler switch to */ /* remove the Unix-specific timing code. */ /* #define NO_TIMER */ /* To insert lots of self-checks for internal errors, define the SELF_CHECK */ /* symbol. This will slow down the program significantly. It is best to */ /* define the symbol using the -DSELF_CHECK compiler switch, but you could */ /* write "#define SELF_CHECK" below. If you are modifying this code, I */ /* recommend you turn self-checks on until your work is debugged. */ /* #define SELF_CHECK */ /* To compile Triangle as a callable object library (triangle.o), define the */ /* TRILIBRARY symbol. Read the file triangle.h for details on how to call */ /* the procedure triangulate() that results. */ /* #define TRILIBRARY */ /* It is possible to generate a smaller version of Triangle using one or */ /* both of the following symbols. Define the REDUCED symbol to eliminate */ /* all features that are primarily of research interest; specifically, the */ /* -i, -F, -s, and -C switches. Define the CDT_ONLY symbol to eliminate */ /* all meshing algorithms above and beyond constrained Delaunay */ /* triangulation; specifically, the -r, -q, -a, -S, and -s switches. */ /* These reductions are most likely to be useful when generating an object */ /* library (triangle.o) by defining the TRILIBRARY symbol. */ /* #define REDUCED */ /* #define CDT_ONLY */ /* On some machines, my exact arithmetic routines might be defeated by the */ /* use of internal extended precision floating-point registers. The best */ /* way to solve this problem is to set the floating-point registers to use */ /* single or double precision internally. On 80x86 processors, this may */ /* be accomplished by setting the CPU86 symbol for the Microsoft C */ /* compiler, or the LINUX symbol for the gcc compiler running on Linux. */ /* */ /* An inferior solution is to declare certain values as `volatile', thus */ /* forcing them to be stored to memory and rounded off. Unfortunately, */ /* this solution might slow Triangle down quite a bit. To use volatile */ /* values, write "#define INEXACT volatile" below. Normally, however, */ /* INEXACT should be defined to be nothing. ("#define INEXACT".) */ /* */ /* For more discussion, see http://www.cs.cmu.edu/~quake/robust.pc.html . */ /* For yet more discussion, see Section 5 of my paper, "Adaptive Precision */ /* Floating-Point Arithmetic and Fast Robust Geometric Predicates" (also */ /* available as Section 6.6 of my dissertation). */ /* #define CPU86 */ /* #define LINUX */ #define INEXACT /* Nothing */ /* #define INEXACT volatile */ /* Maximum number of characters in a file name (including the null). */ #define FILENAMESIZE 512 /* Maximum number of characters in a line read from a file (including the */ /* null). */ #define INPUTLINESIZE 512 /* For efficiency, a variety of data structures are allocated in bulk. The */ /* following constants determine how many of each structure is allocated */ /* at once. */ #define TRIPERBLOCK 4092 /* Number of triangles allocated at once. */ #define SUBSEGPERBLOCK 508 /* Number of subsegments allocated at once. */ #define VERTEXPERBLOCK 4092 /* Number of vertices allocated at once. */ #define VIRUSPERBLOCK 1020 /* Number of virus triangles allocated at once. */ /* Number of encroached subsegments allocated at once. */ #define BADSUBSEGPERBLOCK 252 /* Number of skinny triangles allocated at once. */ #define BADTRIPERBLOCK 4092 /* Number of flipped triangles allocated at once. */ #define FLIPSTACKERPERBLOCK 252 /* Number of splay tree nodes allocated at once. */ #define SPLAYNODEPERBLOCK 508 /* The vertex types. A DEADVERTEX has been deleted entirely. An */ /* UNDEADVERTEX is not part of the mesh, but is written to the output */ /* .node file and affects the node indexing in the other output files. */ #define INPUTVERTEX 0 #define SEGMENTVERTEX 1 #define FREEVERTEX 2 #define DEADVERTEX -32768 #define UNDEADVERTEX -32767 /* The next line is used to outsmart some very stupid compilers. If your */ /* compiler is smarter, feel free to replace the "int" with "void". */ /* Not that it matters. */ #define VOID int /* Two constants for algorithms based on random sampling. Both constants */ /* have been chosen empirically to optimize their respective algorithms. */ /* Used for the point location scheme of Mucke, Saias, and Zhu, to decide */ /* how large a random sample of triangles to inspect. */ #define SAMPLEFACTOR 11 /* Used in Fortune's sweepline Delaunay algorithm to determine what fraction */ /* of boundary edges should be maintained in the splay tree for point */ /* location on the front. */ #define SAMPLERATE 10 /* A number that speaks for itself, every kissable digit. */ #define PI 3.141592653589793238462643383279502884197169399375105820974944592308 /* Another fave. */ #define SQUAREROOTTWO 1.4142135623730950488016887242096980785696718753769480732 /* And here's one for those of you who are intimidated by math. */ #define ONETHIRD 0.333333333333333333333333333333333333333333333333333333333333 #include #include #include #include #ifndef NO_TIMER #include #endif /* not NO_TIMER */ #ifdef CPU86 #include #endif /* CPU86 */ #ifdef LINUX #include #endif /* LINUX */ #ifdef TRILIBRARY #include "triangle.h" #endif /* TRILIBRARY */ /* A few forward declarations. */ #ifndef TRILIBRARY char *readline(); char *findfield(); #endif /* not TRILIBRARY */ /* Labels that signify whether a record consists primarily of pointers or of */ /* floating-point words. Used to make decisions about data alignment. */ enum wordtype {POINTER, FLOATINGPOINT}; /* Labels that signify the result of point location. The result of a */ /* search indicates that the point falls in the interior of a triangle, on */ /* an edge, on a vertex, or outside the mesh. */ enum locateresult {INTRIANGLE, ONEDGE, ONVERTEX, OUTSIDE}; /* Labels that signify the result of vertex insertion. The result indicates */ /* that the vertex was inserted with complete success, was inserted but */ /* encroaches upon a subsegment, was not inserted because it lies on a */ /* segment, or was not inserted because another vertex occupies the same */ /* location. */ enum insertvertexresult {SUCCESSFULVERTEX, ENCROACHINGVERTEX, VIOLATINGVERTEX, DUPLICATEVERTEX}; /* Labels that signify the result of direction finding. The result */ /* indicates that a segment connecting the two query points falls within */ /* the direction triangle, along the left edge of the direction triangle, */ /* or along the right edge of the direction triangle. */ enum finddirectionresult {WITHIN, LEFTCOLLINEAR, RIGHTCOLLINEAR}; /*****************************************************************************/ /* */ /* The basic mesh data structures */ /* */ /* There are three: vertices, triangles, and subsegments (abbreviated */ /* `subseg'). These three data structures, linked by pointers, comprise */ /* the mesh. A vertex simply represents a mesh vertex and its properties. */ /* A triangle is a triangle. A subsegment is a special data structure used */ /* to represent an impenetrable edge of the mesh (perhaps on the outer */ /* boundary, on the boundary of a hole, or part of an internal boundary */ /* separating two triangulated regions). Subsegments represent boundaries, */ /* defined by the user, that triangles may not lie across. */ /* */ /* A triangle consists of a list of three vertices, a list of three */ /* adjoining triangles, a list of three adjoining subsegments (when */ /* segments exist), an arbitrary number of optional user-defined */ /* floating-point attributes, and an optional area constraint. The latter */ /* is an upper bound on the permissible area of each triangle in a region, */ /* used for mesh refinement. */ /* */ /* For a triangle on a boundary of the mesh, some or all of the neighboring */ /* triangles may not be present. For a triangle in the interior of the */ /* mesh, often no neighboring subsegments are present. Such absent */ /* triangles and subsegments are never represented by NULL pointers; they */ /* are represented by two special records: `dummytri', the triangle that */ /* fills "outer space", and `dummysub', the omnipresent subsegment. */ /* `dummytri' and `dummysub' are used for several reasons; for instance, */ /* they can be dereferenced and their contents examined without violating */ /* protected memory. */ /* */ /* However, it is important to understand that a triangle includes other */ /* information as well. The pointers to adjoining vertices, triangles, and */ /* subsegments are ordered in a way that indicates their geometric relation */ /* to each other. Furthermore, each of these pointers contains orientation */ /* information. Each pointer to an adjoining triangle indicates which face */ /* of that triangle is contacted. Similarly, each pointer to an adjoining */ /* subsegment indicates which side of that subsegment is contacted, and how */ /* the subsegment is oriented relative to the triangle. */ /* */ /* The data structure representing a subsegment may be thought to be */ /* abutting the edge of one or two triangle data structures: either */ /* sandwiched between two triangles, or resting against one triangle on an */ /* exterior boundary or hole boundary. */ /* */ /* A subsegment consists of a list of two vertices, a list of two */ /* adjoining subsegments, and a list of two adjoining triangles. One of */ /* the two adjoining triangles may not be present (though there should */ /* always be one), and neighboring subsegments might not be present. */ /* Subsegments also store a user-defined integer "boundary marker". */ /* Typically, this integer is used to indicate what boundary conditions are */ /* to be applied at that location in a finite element simulation. */ /* */ /* Like triangles, subsegments maintain information about the relative */ /* orientation of neighboring objects. */ /* */ /* Vertices are relatively simple. A vertex is a list of floating-point */ /* numbers, starting with the x, and y coordinates, followed by an */ /* arbitrary number of optional user-defined floating-point attributes, */ /* followed by an integer boundary marker. During the segment insertion */ /* phase, there is also a pointer from each vertex to a triangle that may */ /* contain it. Each pointer is not always correct, but when one is, it */ /* speeds up segment insertion. These pointers are assigned values once */ /* at the beginning of the segment insertion phase, and are not used or */ /* updated except during this phase. Edge flipping during segment */ /* insertion will render some of them incorrect. Hence, don't rely upon */ /* them for anything. */ /* */ /* Other than the exception mentioned above, vertices have no information */ /* about what triangles, subfacets, or subsegments they are linked to. */ /* */ /*****************************************************************************/ /*****************************************************************************/ /* */ /* Handles */ /* */ /* The oriented triangle (`otri') and oriented subsegment (`osub') data */ /* structures defined below do not themselves store any part of the mesh. */ /* The mesh itself is made of `triangle's, `subseg's, and `vertex's. */ /* */ /* Oriented triangles and oriented subsegments will usually be referred to */ /* as "handles." A handle is essentially a pointer into the mesh; it */ /* allows you to "hold" one particular part of the mesh. Handles are used */ /* to specify the regions in which one is traversing and modifying the mesh.*/ /* A single `triangle' may be held by many handles, or none at all. (The */ /* latter case is not a memory leak, because the triangle is still */ /* connected to other triangles in the mesh.) */ /* */ /* An `otri' is a handle that holds a triangle. It holds a specific edge */ /* of the triangle. An `osub' is a handle that holds a subsegment. It */ /* holds either the left or right side of the subsegment. */ /* */ /* Navigation about the mesh is accomplished through a set of mesh */ /* manipulation primitives, further below. Many of these primitives take */ /* a handle and produce a new handle that holds the mesh near the first */ /* handle. Other primitives take two handles and glue the corresponding */ /* parts of the mesh together. The orientation of the handles is */ /* important. For instance, when two triangles are glued together by the */ /* bond() primitive, they are glued at the edges on which the handles lie. */ /* */ /* Because vertices have no information about which triangles they are */ /* attached to, I commonly represent a vertex by use of a handle whose */ /* origin is the vertex. A single handle can simultaneously represent a */ /* triangle, an edge, and a vertex. */ /* */ /*****************************************************************************/ /* The triangle data structure. Each triangle contains three pointers to */ /* adjoining triangles, plus three pointers to vertices, plus three */ /* pointers to subsegments (declared below; these pointers are usually */ /* `dummysub'). It may or may not also contain user-defined attributes */ /* and/or a floating-point "area constraint." It may also contain extra */ /* pointers for nodes, when the user asks for high-order elements. */ /* Because the size and structure of a `triangle' is not decided until */ /* runtime, I haven't simply declared the type `triangle' as a struct. */ typedef REAL **triangle; /* Really: typedef triangle *triangle */ /* An oriented triangle: includes a pointer to a triangle and orientation. */ /* The orientation denotes an edge of the triangle. Hence, there are */ /* three possible orientations. By convention, each edge always points */ /* counterclockwise about the corresponding triangle. */ struct otri { triangle *tri; int orient; /* Ranges from 0 to 2. */ }; /* The subsegment data structure. Each subsegment contains two pointers to */ /* adjoining subsegments, plus two pointers to vertices, plus two pointers */ /* to adjoining triangles, plus one boundary marker. */ typedef REAL **subseg; /* Really: typedef subseg *subseg */ /* An oriented subsegment: includes a pointer to a subsegment and an */ /* orientation. The orientation denotes a side of the edge. Hence, there */ /* are two possible orientations. By convention, the edge is always */ /* directed so that the "side" denoted is the right side of the edge. */ struct osub { subseg *ss; int ssorient; /* Ranges from 0 to 1. */ }; /* The vertex data structure. Each vertex is actually an array of REALs. */ /* The number of REALs is unknown until runtime. An integer boundary */ /* marker, and sometimes a pointer to a triangle, is appended after the */ /* REALs. */ typedef REAL *vertex; /* A queue used to store encroached subsegments. Each subsegment's vertices */ /* are stored so that we can check whether a subsegment is still the same. */ struct badsubseg { subseg encsubseg; /* An encroached subsegment. */ vertex subsegorg, subsegdest; /* Its two vertices. */ }; /* A queue used to store bad triangles. The key is the square of the cosine */ /* of the smallest angle of the triangle. Each triangle's vertices are */ /* stored so that one can check whether a triangle is still the same. */ struct badtriang { triangle poortri; /* A skinny or too-large triangle. */ REAL key; /* cos^2 of smallest (apical) angle. */ vertex triangorg, triangdest, triangapex; /* Its three vertices. */ struct badtriang *nexttriang; /* Pointer to next bad triangle. */ }; /* A stack of triangles flipped during the most recent vertex insertion. */ /* The stack is used to undo the vertex insertion if the vertex encroaches */ /* upon a subsegment. */ struct flipstacker { triangle flippedtri; /* A recently flipped triangle. */ struct flipstacker *prevflip; /* Previous flip in the stack. */ }; /* A node in a heap used to store events for the sweepline Delaunay */ /* algorithm. Nodes do not point directly to their parents or children in */ /* the heap. Instead, each node knows its position in the heap, and can */ /* look up its parent and children in a separate array. The `eventptr' */ /* points either to a `vertex' or to a triangle (in encoded format, so */ /* that an orientation is included). In the latter case, the origin of */ /* the oriented triangle is the apex of a "circle event" of the sweepline */ /* algorithm. To distinguish site events from circle events, all circle */ /* events are given an invalid (smaller than `xmin') x-coordinate `xkey'. */ struct event { REAL xkey, ykey; /* Coordinates of the event. */ VOID *eventptr; /* Can be a vertex or the location of a circle event. */ int heapposition; /* Marks this event's position in the heap. */ }; /* A node in the splay tree. Each node holds an oriented ghost triangle */ /* that represents a boundary edge of the growing triangulation. When a */ /* circle event covers two boundary edges with a triangle, so that they */ /* are no longer boundary edges, those edges are not immediately deleted */ /* from the tree; rather, they are lazily deleted when they are next */ /* encountered. (Since only a random sample of boundary edges are kept */ /* in the tree, lazy deletion is faster.) `keydest' is used to verify */ /* that a triangle is still the same as when it entered the splay tree; if */ /* it has been rotated (due to a circle event), it no longer represents a */ /* boundary edge and should be deleted. */ struct splaynode { struct otri keyedge; /* Lprev of an edge on the front. */ vertex keydest; /* Used to verify that splay node is still live. */ struct splaynode *lchild, *rchild; /* Children in splay tree. */ }; /* A type used to allocate memory. firstblock is the first block of items. */ /* nowblock is the block from which items are currently being allocated. */ /* nextitem points to the next slab of free memory for an item. */ /* deaditemstack is the head of a linked list (stack) of deallocated items */ /* that can be recycled. unallocateditems is the number of items that */ /* remain to be allocated from nowblock. */ /* */ /* Traversal is the process of walking through the entire list of items, and */ /* is separate from allocation. Note that a traversal will visit items on */ /* the "deaditemstack" stack as well as live items. pathblock points to */ /* the block currently being traversed. pathitem points to the next item */ /* to be traversed. pathitemsleft is the number of items that remain to */ /* be traversed in pathblock. */ /* */ /* itemwordtype is set to POINTER or FLOATINGPOINT, and is used to suggest */ /* what sort of word the record is primarily made up of. alignbytes */ /* determines how new records should be aligned in memory. itembytes and */ /* itemwords are the length of a record in bytes (after rounding up) and */ /* words. itemsperblock is the number of items allocated at once in a */ /* single block. items is the number of currently allocated items. */ /* maxitems is the maximum number of items that have been allocated at */ /* once; it is the current number of items plus the number of records kept */ /* on deaditemstack. */ struct memorypool { VOID **firstblock, **nowblock; VOID *nextitem; VOID *deaditemstack; VOID **pathblock; VOID *pathitem; enum wordtype itemwordtype; int alignbytes; int itembytes, itemwords; int itemsperblock; long items, maxitems; int unallocateditems; int pathitemsleft; }; /* Global constants. */ REAL splitter; /* Used to split REAL factors for exact multiplication. */ REAL epsilon; /* Floating-point machine epsilon. */ REAL resulterrbound; REAL ccwerrboundA, ccwerrboundB, ccwerrboundC; REAL iccerrboundA, iccerrboundB, iccerrboundC; REAL o3derrboundA, o3derrboundB, o3derrboundC; /* Random number seed is not constant, but I've made it global anyway. */ unsigned long randomseed; /* Current random number seed. */ /* Mesh data structure. Triangle operates on only one mesh, but the mesh */ /* structure is used (instead of global variables) to allow reentrancy. */ struct mesh { /* Variables used to allocate memory for triangles, subsegments, vertices, */ /* viri (triangles being eaten), encroached segments, bad (skinny or too */ /* large) triangles, and splay tree nodes. */ struct memorypool triangles; struct memorypool subsegs; struct memorypool vertices; struct memorypool viri; struct memorypool badsubsegs; struct memorypool badtriangles; struct memorypool flipstackers; struct memorypool splaynodes; /* Variables that maintain the bad triangle queues. The queues are */ /* ordered from 63 (highest priority) to 0 (lowest priority). */ struct badtriang *queuefront[64]; struct badtriang *queuetail[64]; int nextnonemptyq[64]; int firstnonemptyq; /* Variable that maintains the stack of recently flipped triangles. */ struct flipstacker *lastflip; /* Other variables. */ REAL xmin, xmax, ymin, ymax; /* x and y bounds. */ REAL xminextreme; /* Nonexistent x value used as a flag in sweepline. */ int invertices; /* Number of input vertices. */ int inelements; /* Number of input triangles. */ int insegments; /* Number of input segments. */ int holes; /* Number of input holes. */ int regions; /* Number of input regions. */ int undeads; /* Number of input vertices that don't appear in the mesh. */ long edges; /* Number of output edges. */ int mesh_dim; /* Dimension (ought to be 2). */ int nextras; /* Number of attributes per vertex. */ int eextras; /* Number of attributes per triangle. */ long hullsize; /* Number of edges in convex hull. */ int steinerleft; /* Number of Steiner points not yet used. */ int vertexmarkindex; /* Index to find boundary marker of a vertex. */ int vertex2triindex; /* Index to find a triangle adjacent to a vertex. */ int highorderindex; /* Index to find extra nodes for high-order elements. */ int elemattribindex; /* Index to find attributes of a triangle. */ int areaboundindex; /* Index to find area bound of a triangle. */ int checksegments; /* Are there segments in the triangulation yet? */ int checkquality; /* Has quality triangulation begun yet? */ int readnodefile; /* Has a .node file been read? */ long samples; /* Number of random samples for point location. */ long incirclecount; /* Number of incircle tests performed. */ long counterclockcount; /* Number of counterclockwise tests performed. */ long orient3dcount; /* Number of 3D orientation tests performed. */ long hyperbolacount; /* Number of right-of-hyperbola tests performed. */ long circumcentercount; /* Number of circumcenter calculations performed. */ long circletopcount; /* Number of circle top calculations performed. */ /* Triangular bounding box vertices. */ vertex infvertex1, infvertex2, infvertex3; /* Pointer to the `triangle' that occupies all of "outer space." */ triangle *dummytri; triangle *dummytribase; /* Keep base address so we can free() it later. */ /* Pointer to the omnipresent subsegment. Referenced by any triangle or */ /* subsegment that isn't really connected to a subsegment at that */ /* location. */ subseg *dummysub; subseg *dummysubbase; /* Keep base address so we can free() it later. */ /* Pointer to a recently visited triangle. Improves point location if */ /* proximate vertices are inserted sequentially. */ struct otri recenttri; }; /* End of `struct mesh'. */ /* Data structure for command line switches and file names. This structure */ /* is used (instead of global variables) to allow reentrancy. */ struct behavior { /* Switches for the triangulator. */ /* poly: -p switch. refine: -r switch. */ /* quality: -q switch. */ /* minangle: minimum angle bound, specified after -q switch. */ /* goodangle: cosine squared of minangle. */ /* vararea: -a switch without number. */ /* fixedarea: -a switch with number. */ /* maxarea: maximum area bound, specified after -a switch. */ /* usertest: -u switch. */ /* regionattrib: -A switch. convex: -c switch. */ /* weighted: 1 for -w switch, 2 for -W switch. jettison: -j switch */ /* firstnumber: inverse of -z switch. All items are numbered starting */ /* from `firstnumber'. */ /* edgesout: -e switch. voronoi: -v switch. */ /* neighbors: -n switch. geomview: -g switch. */ /* nobound: -B switch. nopolywritten: -P switch. */ /* nonodewritten: -N switch. noelewritten: -E switch. */ /* noiterationnum: -I switch. noholes: -O switch. */ /* noexact: -X switch. */ /* order: element order, specified after -o switch. */ /* nobisect: count of how often -Y switch is selected. */ /* steiner: maximum number of Steiner points, specified after -S switch. */ /* incremental: -i switch. sweepline: -F switch. */ /* dwyer: inverse of -l switch. */ /* splitseg: -s switch. */ /* nolenses: -L switch. docheck: -C switch. */ /* quiet: -Q switch. verbose: count of how often -V switch is selected. */ /* usesegments: -p, -r, -q, or -c switch; determines whether segments are */ /* used at all. */ /* */ /* Read the instructions to find out the meaning of these switches. */ int poly, refine, quality, vararea, fixedarea, usertest; int regionattrib, convex, weighted, jettison; int firstnumber; int edgesout, voronoi, neighbors, geomview; int nobound, nopolywritten, nonodewritten, noelewritten, noiterationnum; int noholes, noexact, nolenses; int incremental, sweepline, dwyer; int splitseg; int docheck; int quiet, verbose; int usesegments; int order; int nobisect; int steiner; REAL minangle, goodangle; REAL maxarea; /* Variables for file names. */ #ifndef TRILIBRARY char innodefilename[FILENAMESIZE]; char inelefilename[FILENAMESIZE]; char inpolyfilename[FILENAMESIZE]; char areafilename[FILENAMESIZE]; char outnodefilename[FILENAMESIZE]; char outelefilename[FILENAMESIZE]; char outpolyfilename[FILENAMESIZE]; char edgefilename[FILENAMESIZE]; char vnodefilename[FILENAMESIZE]; char vedgefilename[FILENAMESIZE]; char neighborfilename[FILENAMESIZE]; char offfilename[FILENAMESIZE]; #endif /* not TRILIBRARY */ }; /* End of `struct behavior'. */ /*****************************************************************************/ /* */ /* Mesh manipulation primitives. Each triangle contains three pointers to */ /* other triangles, with orientations. Each pointer points not to the */ /* first byte of a triangle, but to one of the first three bytes of a */ /* triangle. It is necessary to extract both the triangle itself and the */ /* orientation. To save memory, I keep both pieces of information in one */ /* pointer. To make this possible, I assume that all triangles are aligned */ /* to four-byte boundaries. The decode() routine below decodes a pointer, */ /* extracting an orientation (in the range 0 to 2) and a pointer to the */ /* beginning of a triangle. The encode() routine compresses a pointer to a */ /* triangle and an orientation into a single pointer. My assumptions that */ /* triangles are four-byte-aligned and that the `unsigned long' type is */ /* long enough to hold a pointer are two of the few kludges in this program.*/ /* */ /* Subsegments are manipulated similarly. A pointer to a subsegment */ /* carries both an address and an orientation in the range 0 to 1. */ /* */ /* The other primitives take an oriented triangle or oriented subsegment, */ /* and return an oriented triangle or oriented subsegment or vertex; or */ /* they change the connections in the data structure. */ /* */ /* Below, triangles and subsegments are denoted by their vertices. The */ /* triangle abc has origin (org) a, destination (dest) b, and apex (apex) */ /* c. These vertices occur in counterclockwise order about the triangle. */ /* The handle abc may simultaneously denote vertex a, edge ab, and triangle */ /* abc. */ /* */ /* Similarly, the subsegment ab has origin (sorg) a and destination (sdest) */ /* b. If ab is thought to be directed upward (with b directly above a), */ /* then the handle ab is thought to grasp the right side of ab, and may */ /* simultaneously denote vertex a and edge ab. */ /* */ /* An asterisk (*) denotes a vertex whose identity is unknown. */ /* */ /* Given this notation, a partial list of mesh manipulation primitives */ /* follows. */ /* */ /* */ /* For triangles: */ /* */ /* sym: Find the abutting triangle; same edge. */ /* sym(abc) -> ba* */ /* */ /* lnext: Find the next edge (counterclockwise) of a triangle. */ /* lnext(abc) -> bca */ /* */ /* lprev: Find the previous edge (clockwise) of a triangle. */ /* lprev(abc) -> cab */ /* */ /* onext: Find the next edge counterclockwise with the same origin. */ /* onext(abc) -> ac* */ /* */ /* oprev: Find the next edge clockwise with the same origin. */ /* oprev(abc) -> a*b */ /* */ /* dnext: Find the next edge counterclockwise with the same destination. */ /* dnext(abc) -> *ba */ /* */ /* dprev: Find the next edge clockwise with the same destination. */ /* dprev(abc) -> cb* */ /* */ /* rnext: Find the next edge (counterclockwise) of the adjacent triangle. */ /* rnext(abc) -> *a* */ /* */ /* rprev: Find the previous edge (clockwise) of the adjacent triangle. */ /* rprev(abc) -> b** */ /* */ /* org: Origin dest: Destination apex: Apex */ /* org(abc) -> a dest(abc) -> b apex(abc) -> c */ /* */ /* bond: Bond two triangles together at the resepective handles. */ /* bond(abc, bad) */ /* */ /* */ /* For subsegments: */ /* */ /* ssym: Reverse the orientation of a subsegment. */ /* ssym(ab) -> ba */ /* */ /* spivot: Find adjoining subsegment with the same origin. */ /* spivot(ab) -> a* */ /* */ /* snext: Find next subsegment in sequence. */ /* snext(ab) -> b* */ /* */ /* sorg: Origin sdest: Destination */ /* sorg(ab) -> a sdest(ab) -> b */ /* */ /* sbond: Bond two subsegments together at the respective origins. */ /* sbond(ab, ac) */ /* */ /* */ /* For interacting tetrahedra and subfacets: */ /* */ /* tspivot: Find a subsegment abutting a triangle. */ /* tspivot(abc) -> ba */ /* */ /* stpivot: Find a triangle abutting a subsegment. */ /* stpivot(ab) -> ba* */ /* */ /* tsbond: Bond a triangle to a subsegment. */ /* tsbond(abc, ba) */ /* */ /*****************************************************************************/ /********* Mesh manipulation primitives begin here *********/ /** **/ /** **/ /* Fast lookup arrays to speed some of the mesh manipulation primitives. */ int plus1mod3[3] = {1, 2, 0}; int minus1mod3[3] = {2, 0, 1}; /********* Primitives for triangles *********/ /* */ /* */ /* decode() converts a pointer to an oriented triangle. The orientation is */ /* extracted from the two least significant bits of the pointer. */ #define decode(ptr, otri) \ (otri).orient = (int) ((unsigned long) (ptr) & (unsigned long) 3l); \ (otri).tri = (triangle *) \ ((unsigned long) (ptr) ^ (unsigned long) (otri).orient) /* encode() compresses an oriented triangle into a single pointer. It */ /* relies on the assumption that all triangles are aligned to four-byte */ /* boundaries, so the two least significant bits of (otri).tri are zero. */ #define encode(otri) \ (triangle) ((unsigned long) (otri).tri | (unsigned long) (otri).orient) /* The following handle manipulation primitives are all described by Guibas */ /* and Stolfi. However, Guibas and Stolfi use an edge-based data */ /* structure, whereas I use a triangle-based data structure. */ /* sym() finds the abutting triangle, on the same edge. Note that the edge */ /* direction is necessarily reversed, because the handle specified by an */ /* oriented triangle is directed counterclockwise around the triangle. */ #define sym(otri1, otri2) \ ptr = (otri1).tri[(otri1).orient]; \ decode(ptr, otri2); #define symself(otri) \ ptr = (otri).tri[(otri).orient]; \ decode(ptr, otri); /* lnext() finds the next edge (counterclockwise) of a triangle. */ #define lnext(otri1, otri2) \ (otri2).tri = (otri1).tri; \ (otri2).orient = plus1mod3[(otri1).orient] #define lnextself(otri) \ (otri).orient = plus1mod3[(otri).orient] /* lprev() finds the previous edge (clockwise) of a triangle. */ #define lprev(otri1, otri2) \ (otri2).tri = (otri1).tri; \ (otri2).orient = minus1mod3[(otri1).orient] #define lprevself(otri) \ (otri).orient = minus1mod3[(otri).orient] /* onext() spins counterclockwise around a vertex; that is, it finds the */ /* next edge with the same origin in the counterclockwise direction. This */ /* edge is part of a different triangle. */ #define onext(otri1, otri2) \ lprev(otri1, otri2); \ symself(otri2); #define onextself(otri) \ lprevself(otri); \ symself(otri); /* oprev() spins clockwise around a vertex; that is, it finds the next edge */ /* with the same origin in the clockwise direction. This edge is part of */ /* a different triangle. */ #define oprev(otri1, otri2) \ sym(otri1, otri2); \ lnextself(otri2); #define oprevself(otri) \ symself(otri); \ lnextself(otri); /* dnext() spins counterclockwise around a vertex; that is, it finds the */ /* next edge with the same destination in the counterclockwise direction. */ /* This edge is part of a different triangle. */ #define dnext(otri1, otri2) \ sym(otri1, otri2); \ lprevself(otri2); #define dnextself(otri) \ symself(otri); \ lprevself(otri); /* dprev() spins clockwise around a vertex; that is, it finds the next edge */ /* with the same destination in the clockwise direction. This edge is */ /* part of a different triangle. */ #define dprev(otri1, otri2) \ lnext(otri1, otri2); \ symself(otri2); #define dprevself(otri) \ lnextself(otri); \ symself(otri); /* rnext() moves one edge counterclockwise about the adjacent triangle. */ /* (It's best understood by reading Guibas and Stolfi. It involves */ /* changing triangles twice.) */ #define rnext(otri1, otri2) \ sym(otri1, otri2); \ lnextself(otri2); \ symself(otri2); #define rnextself(otri) \ symself(otri); \ lnextself(otri); \ symself(otri); /* rprev() moves one edge clockwise about the adjacent triangle. */ /* (It's best understood by reading Guibas and Stolfi. It involves */ /* changing triangles twice.) */ #define rprev(otri1, otri2) \ sym(otri1, otri2); \ lprevself(otri2); \ symself(otri2); #define rprevself(otri) \ symself(otri); \ lprevself(otri); \ symself(otri); /* These primitives determine or set the origin, destination, or apex of a */ /* triangle. */ #define org(otri, vertexptr) \ vertexptr = (vertex) (otri).tri[plus1mod3[(otri).orient] + 3] #define dest(otri, vertexptr) \ vertexptr = (vertex) (otri).tri[minus1mod3[(otri).orient] + 3] #define apex(otri, vertexptr) \ vertexptr = (vertex) (otri).tri[(otri).orient + 3] #define setorg(otri, vertexptr) \ (otri).tri[plus1mod3[(otri).orient] + 3] = (triangle) vertexptr #define setdest(otri, vertexptr) \ (otri).tri[minus1mod3[(otri).orient] + 3] = (triangle) vertexptr #define setapex(otri, vertexptr) \ (otri).tri[(otri).orient + 3] = (triangle) vertexptr /* Bond two triangles together. */ #define bond(otri1, otri2) \ (otri1).tri[(otri1).orient] = encode(otri2); \ (otri2).tri[(otri2).orient] = encode(otri1) /* Dissolve a bond (from one side). Note that the other triangle will still */ /* think it's connected to this triangle. Usually, however, the other */ /* triangle is being deleted entirely, or bonded to another triangle, so */ /* it doesn't matter. */ #define dissolve(otri) \ (otri).tri[(otri).orient] = (triangle) m->dummytri /* Copy an oriented triangle. */ #define otricopy(otri1, otri2) \ (otri2).tri = (otri1).tri; \ (otri2).orient = (otri1).orient /* Test for equality of oriented triangles. */ #define otriequal(otri1, otri2) \ (((otri1).tri == (otri2).tri) && \ ((otri1).orient == (otri2).orient)) /* Primitives to infect or cure a triangle with the virus. These rely on */ /* the assumption that all subsegments are aligned to four-byte boundaries.*/ #define infect(otri) \ (otri).tri[6] = (triangle) \ ((unsigned long) (otri).tri[6] | (unsigned long) 2l) #define uninfect(otri) \ (otri).tri[6] = (triangle) \ ((unsigned long) (otri).tri[6] & ~ (unsigned long) 2l) /* Test a triangle for viral infection. */ #define infected(otri) \ (((unsigned long) (otri).tri[6] & (unsigned long) 2l) != 0l) /* Check or set a triangle's attributes. */ #define elemattribute(otri, attnum) \ ((REAL *) (otri).tri)[m->elemattribindex + (attnum)] #define setelemattribute(otri, attnum, value) \ ((REAL *) (otri).tri)[m->elemattribindex + (attnum)] = value /* Check or set a triangle's maximum area bound. */ #define areabound(otri) ((REAL *) (otri).tri)[m->areaboundindex] #define setareabound(otri, value) \ ((REAL *) (otri).tri)[m->areaboundindex] = value /* Check or set a triangle's deallocation. Its second pointer is set to */ /* NULL to indicate that it is not allocated. (Its first pointer is used */ /* for the stack of dead items.) Its fourth pointer (its first vertex) */ /* is set to NULL in case a `badtriang' structure points to it. */ #define deadtri(tria) ((tria)[1] == (triangle) NULL) #define killtri(tria) \ (tria)[1] = (triangle) NULL; \ (tria)[3] = (triangle) NULL /********* Primitives for subsegments *********/ /* */ /* */ /* sdecode() converts a pointer to an oriented subsegment. The orientation */ /* is extracted from the least significant bit of the pointer. The two */ /* least significant bits (one for orientation, one for viral infection) */ /* are masked out to produce the real pointer. */ #define sdecode(sptr, osub) \ (osub).ssorient = (int) ((unsigned long) (sptr) & (unsigned long) 1l); \ (osub).ss = (subseg *) \ ((unsigned long) (sptr) & ~ (unsigned long) 3l) /* sencode() compresses an oriented subsegment into a single pointer. It */ /* relies on the assumption that all subsegments are aligned to two-byte */ /* boundaries, so the least significant bit of (osub).ss is zero. */ #define sencode(osub) \ (subseg) ((unsigned long) (osub).ss | (unsigned long) (osub).ssorient) /* ssym() toggles the orientation of a subsegment. */ #define ssym(osub1, osub2) \ (osub2).ss = (osub1).ss; \ (osub2).ssorient = 1 - (osub1).ssorient #define ssymself(osub) \ (osub).ssorient = 1 - (osub).ssorient /* spivot() finds the other subsegment (from the same segment) that shares */ /* the same origin. */ #define spivot(osub1, osub2) \ sptr = (osub1).ss[(osub1).ssorient]; \ sdecode(sptr, osub2) #define spivotself(osub) \ sptr = (osub).ss[(osub).ssorient]; \ sdecode(sptr, osub) /* snext() finds the next subsegment (from the same segment) in sequence; */ /* one whose origin is the input subsegment's destination. */ #define snext(osub1, osub2) \ sptr = (osub1).ss[1 - (osub1).ssorient]; \ sdecode(sptr, osub2) #define snextself(osub) \ sptr = (osub).ss[1 - (osub).ssorient]; \ sdecode(sptr, osub) /* These primitives determine or set the origin or destination of a */ /* subsegment. */ #define sorg(osub, vertexptr) \ vertexptr = (vertex) (osub).ss[2 + (osub).ssorient] #define sdest(osub, vertexptr) \ vertexptr = (vertex) (osub).ss[3 - (osub).ssorient] #define setsorg(osub, vertexptr) \ (osub).ss[2 + (osub).ssorient] = (subseg) vertexptr #define setsdest(osub, vertexptr) \ (osub).ss[3 - (osub).ssorient] = (subseg) vertexptr /* These primitives read or set a boundary marker. Boundary markers are */ /* used to hold user-defined tags for setting boundary conditions in */ /* finite element solvers. */ #define mark(osub) (* (int *) ((osub).ss + 6)) #define setmark(osub, value) \ * (int *) ((osub).ss + 6) = value /* Bond two subsegments together. */ #define sbond(osub1, osub2) \ (osub1).ss[(osub1).ssorient] = sencode(osub2); \ (osub2).ss[(osub2).ssorient] = sencode(osub1) /* Dissolve a subsegment bond (from one side). Note that the other */ /* subsegment will still think it's connected to this subsegment. */ #define sdissolve(osub) \ (osub).ss[(osub).ssorient] = (subseg) m->dummysub /* Copy a subsegment. */ #define subsegcopy(osub1, osub2) \ (osub2).ss = (osub1).ss; \ (osub2).ssorient = (osub1).ssorient /* Test for equality of subsegments. */ #define subsegequal(osub1, osub2) \ (((osub1).ss == (osub2).ss) && \ ((osub1).ssorient == (osub2).ssorient)) /* Check or set a subsegment's deallocation. Its second pointer is set to */ /* NULL to indicate that it is not allocated. (Its first pointer is used */ /* for the stack of dead items.) Its third pointer (its first vertex) */ /* is set to NULL in case a `badsubseg' structure points to it. */ #define deadsubseg(sub) ((sub)[1] == (subseg) NULL) #define killsubseg(sub) \ (sub)[1] = (subseg) NULL; \ (sub)[2] = (subseg) NULL /********* Primitives for interacting triangles and subsegments *********/ /* */ /* */ /* tspivot() finds a subsegment abutting a triangle. */ #define tspivot(otri, osub) \ sptr = (subseg) (otri).tri[6 + (otri).orient]; \ sdecode(sptr, osub) /* stpivot() finds a triangle abutting a subsegment. It requires that the */ /* variable `ptr' of type `triangle' be defined. */ #define stpivot(osub, otri) \ ptr = (triangle) (osub).ss[4 + (osub).ssorient]; \ decode(ptr, otri) /* Bond a triangle to a subsegment. */ #define tsbond(otri, osub) \ (otri).tri[6 + (otri).orient] = (triangle) sencode(osub); \ (osub).ss[4 + (osub).ssorient] = (subseg) encode(otri) /* Dissolve a bond (from the triangle side). */ #define tsdissolve(otri) \ (otri).tri[6 + (otri).orient] = (triangle) m->dummysub /* Dissolve a bond (from the subsegment side). */ #define stdissolve(osub) \ (osub).ss[4 + (osub).ssorient] = (subseg) m->dummytri /********* Primitives for vertices *********/ /* */ /* */ #define vertexmark(vx) ((int *) (vx))[m->vertexmarkindex] #define setvertexmark(vx, value) \ ((int *) (vx))[m->vertexmarkindex] = value #define vertextype(vx) ((int *) (vx))[m->vertexmarkindex + 1] #define setvertextype(vx, value) \ ((int *) (vx))[m->vertexmarkindex + 1] = value #define vertex2tri(vx) ((triangle *) (vx))[m->vertex2triindex] #define setvertex2tri(vx, value) \ ((triangle *) (vx))[m->vertex2triindex] = value /** **/ /** **/ /********* Mesh manipulation primitives end here *********/ /********* User-defined triangle evaluation routine begins here *********/ /** **/ /** **/ /*****************************************************************************/ /* */ /* triunsuitable() Determine if a triangle is unsuitable, and thus must */ /* be further refined. */ /* */ /* You may write your own procedure that decides whether or not a selected */ /* triangle is too big (and needs to be refined). There are two ways to do */ /* this. */ /* */ /* (1) Modify the procedure `triunsuitable' below, then recompile */ /* Triangle. */ /* */ /* (2) Define the symbol EXTERNAL_TEST (either by adding the definition */ /* to this file, or by using the appropriate compiler switch). This way, */ /* you can compile triangle.c separately from your test. Write your own */ /* `triunsuitable' procedure in a separate C file (using the same prototype */ /* as below). Compile it and link the object code with triangle.o. */ /* */ /* This procedure returns 1 if the triangle is too large and should be */ /* refined; 0 otherwise. */ /* */ /*****************************************************************************/ #ifdef EXTERNAL_TEST #ifdef ANSI_DECLARATORS extern int triunsuitable(vertex triorg, vertex tridest, vertex triapex, REAL area); #else /* not ANSI_DECLARATORS */ extern int triunsuitable(); #endif /* not ANSI_DECLARATORS */ #else /* not EXTERNAL_TEST */ #ifdef ANSI_DECLARATORS int triunsuitable(vertex triorg, vertex tridest, vertex triapex, REAL area) #else /* not ANSI_DECLARATORS */ int triunsuitable(triorg, tridest, triapex, area) vertex triorg; /* The triangle's origin vertex. */ vertex tridest; /* The triangle's destination vertex. */ vertex triapex; /* The triangle's apex vertex. */ REAL area; /* The area of the triangle. */ #endif /* not ANSI_DECLARATORS */ { REAL dxoa, dxda, dxod; REAL dyoa, dyda, dyod; REAL oalen, dalen, odlen; REAL maxlen; dxoa = triorg[0] - triapex[0]; dyoa = triorg[1] - triapex[1]; dxda = tridest[0] - triapex[0]; dyda = tridest[1] - triapex[1]; dxod = triorg[0] - tridest[0]; dyod = triorg[1] - tridest[1]; /* Find the squares of the lengths of the triangle's three edges. */ oalen = dxoa * dxoa + dyoa * dyoa; dalen = dxda * dxda + dyda * dyda; odlen = dxod * dxod + dyod * dyod; /* Find the square of the length of the longest edge. */ maxlen = (dalen > oalen) ? dalen : oalen; maxlen = (odlen > maxlen) ? odlen : maxlen; if (maxlen > 0.05 * (triorg[0] * triorg[0] + triorg[1] * triorg[1]) + 0.02) { return 1; } else { return 0; } } #endif /* not EXTERNAL_TEST */ /** **/ /** **/ /********* User-defined triangle evaluation routine ends here *********/ /********* Memory allocation wrappers begin here *********/ /** **/ /** **/ #ifdef ANSI_DECLARATORS VOID *trimalloc(int size) #else /* not ANSI_DECLARATORS */ VOID *trimalloc(size) int size; #endif /* not ANSI_DECLARATORS */ { VOID *memptr; memptr = malloc(size); if (memptr == (VOID *) NULL) { printf("Error: Out of memory.\n"); exit(1); } return(memptr); } #ifdef ANSI_DECLARATORS void trifree(VOID *memptr) #else /* not ANSI_DECLARATORS */ void trifree(memptr) VOID *memptr; #endif /* not ANSI_DECLARATORS */ { free(memptr); } /** **/ /** **/ /********* Memory allocation wrappers end here *********/ /********* User interaction routines begin here *********/ /** **/ /** **/ /*****************************************************************************/ /* */ /* syntax() Print list of command line switches. */ /* */ /*****************************************************************************/ #ifndef TRILIBRARY void syntax() { #ifdef CDT_ONLY #ifdef REDUCED printf("triangle [-pAcjevngBPNEIOXzo_lQVh] input_file\n"); #else /* not REDUCED */ printf("triangle [-pAcjevngBPNEIOXzo_iFlCQVh] input_file\n"); #endif /* not REDUCED */ #else /* not CDT_ONLY */ #ifdef REDUCED printf("triangle [-prq__a__uAcjevngBPNEIOXzo_YS__LlQVh] input_file\n"); #else /* not REDUCED */ printf("triangle [-prq__a__uAcjevngBPNEIOXzo_YS__LiFlsCQVh] input_file\n"); #endif /* not REDUCED */ #endif /* not CDT_ONLY */ printf(" -p Triangulates a Planar Straight Line Graph (.poly file).\n"); #ifndef CDT_ONLY printf(" -r Refines a previously generated mesh.\n"); printf( " -q Quality mesh generation. A minimum angle may be specified.\n"); printf(" -a Applies a maximum triangle area constraint.\n"); printf(" -u Applies a user-defined triangle constraint.\n"); #endif /* not CDT_ONLY */ printf( " -A Applies attributes to identify triangles in certain regions.\n"); printf(" -c Encloses the convex hull with segments.\n"); printf(" -w Weighted Delaunay triangulation.\n"); printf(" -W Regular triangulation (lower hull of a height field).\n"); printf(" -j Jettison unused vertices from output .node file.\n"); printf(" -e Generates an edge list.\n"); printf(" -v Generates a Voronoi diagram.\n"); printf(" -n Generates a list of triangle neighbors.\n"); printf(" -g Generates an .off file for Geomview.\n"); printf(" -B Suppresses output of boundary information.\n"); printf(" -P Suppresses output of .poly file.\n"); printf(" -N Suppresses output of .node file.\n"); printf(" -E Suppresses output of .ele file.\n"); printf(" -I Suppresses mesh iteration numbers.\n"); printf(" -O Ignores holes in .poly file.\n"); printf(" -X Suppresses use of exact arithmetic.\n"); printf(" -z Numbers all items starting from zero (rather than one).\n"); printf(" -o2 Generates second-order subparametric elements.\n"); #ifndef CDT_ONLY printf(" -Y Suppresses boundary segment splitting.\n"); printf(" -S Specifies maximum number of added Steiner points.\n"); printf(" -L Uses equatorial circles, not equatorial lenses.\n"); #endif /* not CDT_ONLY */ #ifndef REDUCED printf(" -i Uses incremental method, rather than divide-and-conquer.\n"); printf(" -F Uses Fortune's sweepline algorithm, rather than d-and-c.\n"); #endif /* not REDUCED */ printf(" -l Uses vertical cuts only, rather than alternating cuts.\n"); #ifndef REDUCED #ifndef CDT_ONLY printf( " -s Force segments into mesh by splitting (instead of using CDT).\n"); printf(" -L Uses Ruppert's diametral spheres, not diametral lenses.\n"); #endif /* not CDT_ONLY */ printf(" -C Check consistency of final mesh.\n"); #endif /* not REDUCED */ printf(" -Q Quiet: No terminal output except errors.\n"); printf(" -V Verbose: Detailed information on what I'm doing.\n"); printf(" -h Help: Detailed instructions for Triangle.\n"); exit(0); } #endif /* not TRILIBRARY */ /*****************************************************************************/ /* */ /* info() Print out complete instructions. */ /* */ /*****************************************************************************/ #ifndef TRILIBRARY void info() { printf("Triangle\n"); printf( "A Two-Dimensional Quality Mesh Generator and Delaunay Triangulator.\n"); printf("Version 1.4\n\n"); printf("Copyright 1993, 1995, 1997, 1998, 2002 Jonathan Richard Shewchuk\n"); printf("2360 Woolsey #H / Berkeley, California 94705-1927\n"); printf("Bugs/comments to jrs@cs.berkeley.edu\n"); printf( "Created as part of the Archimedes project (tools for parallel FEM).\n"); printf( "Supported in part by NSF Grant CMS-9318163 and an NSERC 1967 Scholarship.\n"); printf("There is no warranty whatsoever. Use at your own risk.\n"); #ifdef SINGLE printf("This executable is compiled for single precision arithmetic.\n\n\n"); #else /* not SINGLE */ printf("This executable is compiled for double precision arithmetic.\n\n\n"); #endif /* not SINGLE */ printf( "Triangle generates exact Delaunay triangulations, constrained Delaunay\n"); printf( "triangulations, Voronoi diagrams, and quality conforming Delaunay\n"); printf( "triangulations. The latter can be generated with no small angles, and are\n" ); printf( "thus suitable for finite element analysis. If no command line switches are\n" ); printf( "specified, your .node input file is read, and the Delaunay triangulation is\n" ); printf("returned in .node and .ele output files. The command syntax is:\n"); printf("\n"); printf("triangle [-prq__a__uAcjevngBPNEIOXzo_YS__LiFlsCQVh] input_file\n"); printf("\n"); printf( "Underscores indicate that numbers may optionally follow certain switches.\n"); printf( "Do not leave any space between a switch and its numeric parameter.\n"); printf( "input_file must be a file with extension .node, or extension .poly if the\n"); printf( "-p switch is used. If -r is used, you must supply .node and .ele files,\n"); printf( "and possibly a .poly file and an .area file as well. The formats of these\n" ); printf("files are described below.\n\n"); printf("Command Line Switches:\n\n"); printf( " -p Reads a Planar Straight Line Graph (.poly file), which can specify\n" ); printf(" vertices, segments, holes, regional attributes, and area\n"); printf( " constraints. Generates a constrained Delaunay triangulation (CDT)\n" ); printf( " fitting the input; or, if -s, -q, -a, or -u is used, a conforming\n"); printf( " constrained Delaunay triangulation (CCDT). If -p is not used,\n"); printf(" Triangle reads a .node file by default.\n"); printf( " -r Refines a previously generated mesh. The mesh is read from a .node\n" ); printf( " file and an .ele file. If -p is also used, a .poly file is read\n"); printf( " and used to constrain segments in the mesh. If -a is also used\n"); printf( " (with no number following), an .area file is read and used to\n"); printf( " impose area constraints on the mesh. Further details on refinement\n" ); printf(" are given below.\n"); printf( " -q Quality mesh generation by my variant of Jim Ruppert's Delaunay\n"); printf( " refinement algorithm. Adds vertices to the mesh to ensure that no\n" ); printf( " angles smaller than 20 degrees occur. An alternative minimum angle\n" ); printf( " may be specified after the `q'. If the minimum angle is 20.7\n"); printf( " degrees or smaller, the triangulation algorithm is mathematically\n"); printf( " guaranteed to terminate (assuming infinite precision arithmetic--\n"); printf( " Triangle may fail to terminate if you run out of precision). In\n"); printf( " practice, the algorithm often succeeds for minimum angles up to\n"); printf( " 33.8 degrees. For some meshes, however, it may be necessary to\n"); printf( " reduce the minimum angle to avoid problems associated with\n"); printf( " insufficient floating-point precision. The specified angle may\n"); printf(" include a decimal point.\n"); printf( " -a Imposes a maximum triangle area. If a number follows the `a', no\n"); printf( " triangle is generated whose area is larger than that number. If no\n" ); printf( " number is specified, an .area file (if -r is used) or .poly file\n"); printf( " (if -r is not used) specifies a set of maximum area constraints.\n"); printf( " An .area file contains a separate area constraint for each\n"); printf( " triangle, and is useful for refining a finite element mesh based on\n" ); printf( " a posteriori error estimates. A .poly file can optionally contain\n" ); printf( " an area constraint for each segment-bounded region, thereby\n"); printf( " controlling triangle densities in a first triangulation of a PSLG.\n" ); printf( " You can impose both a fixed area constraint and a varying area\n"); printf( " constraint by invoking the -a switch twice, once with and once\n"); printf( " without a number following. Each area specified may include a\n"); printf(" decimal point.\n"); printf( " -u Imposes a user-defined constraint on triangle size. There are two\n" ); printf( " ways to use this feature. One is to edit the triunsuitable()\n"); printf( " procedure in triangle.c to encode any constraint you like, then\n"); printf( " recompile Triangle. The other is to compile triangle.c with the\n"); printf( " EXTERNAL_TEST symbol set (compiler switch -DEXTERNAL_TEST), then\n"); printf( " link Triangle against a separate object file that implements\n"); printf( " triunsuitable(). In either case, the -u switch causes the user-\n"); printf(" defined test to be applied to every triangle.\n"); printf( " -A Assigns an additional attribute to each triangle that identifies\n"); printf( " what segment-bounded region each triangle belongs to. Attributes\n"); printf( " are assigned to regions by the .poly file. If a region is not\n"); printf( " explicitly marked by the .poly file, triangles in that region are\n"); printf( " assigned an attribute of zero. The -A switch has an effect only\n"); printf(" when the -p switch is used and the -r switch is not.\n"); printf( " -c Creates segments on the convex hull of the triangulation. If you\n"); printf( " are triangulating a vertex set, this switch causes a .poly file to\n" ); printf( " be written, containing all edges in the convex hull. If you are\n"); printf( " triangulating a PSLG, this switch specifies that the whole convex\n"); printf( " hull of the PSLG should be triangulated, regardless of what\n"); printf( " segments the PSLG has. If you do not use this switch when\n"); printf( " triangulating a PSLG, it is assumed that you have identified the\n"); printf( " region to be triangulated by surrounding it with segments of the\n"); printf( " input PSLG. Beware: if you are not careful, this switch can cause\n" ); printf( " the introduction of an extremely thin angle between a PSLG segment\n" ); printf( " and a convex hull segment, which can cause overrefinement (and\n"); printf( " possibly failure if Triangle runs out of precision). If you are\n"); printf( " refining a mesh, the -c switch works differently; it generates the\n" ); printf( " set of boundary edges of the mesh (useful if no .poly file was\n"); printf(" read).\n"); printf( " -j Jettisons vertices that are not part of the final triangulation\n"); printf( " from the output .node file. By default, Triangle copies all\n"); printf( " vertices in the input .node file to the output .node file, in the\n"); printf( " same order, so their indices do not change. The -j switch prevents\n" ); printf( " duplicated input vertices from appearing in the output .node file;\n" ); printf( " hence, if two input vertices have exactly the same coordinates,\n"); printf( " only the first appears in the output. If any vertices are\n"); printf( " jettisoned, the vertex numbering in the output .node file differs\n"); printf(" from that of the input .node file.\n"); printf( " -e Outputs (to an .edge file) a list of edges of the triangulation.\n"); printf( " -v Outputs the Voronoi diagram associated with the triangulation.\n"); printf( " Does not attempt to detect degeneracies, so some Voronoi vertices\n"); printf( " may be duplicated. See the discussion of Voronoi diagrams below.\n"); printf( " -n Outputs (to a .neigh file) a list of triangles neighboring each\n"); printf(" triangle.\n"); printf( " -g Outputs the mesh to an Object File Format (.off) file, suitable for\n" ); printf(" viewing with the Geometry Center's Geomview package.\n"); printf( " -B No boundary markers in the output .node, .poly, and .edge output\n"); printf( " files. See the detailed discussion of boundary markers below.\n"); printf( " -P No output .poly file. Saves disk space, but you lose the ability\n"); printf( " to maintain constraining segments on later refinements of the mesh.\n" ); printf(" -N No output .node file.\n"); printf(" -E No output .ele file.\n"); printf( " -I No iteration numbers. Suppresses the output of .node and .poly\n"); printf( " files, so your input files won't be overwritten. (If your input is\n" ); printf( " a .poly file only, a .node file is written.) Cannot be used with\n"); printf( " the -r switch, because that would overwrite your input .ele file.\n"); printf( " Shouldn't be used with the -q, -a, -u, or -s switch if you are\n"); printf( " using a .node file for input, because no .node file is written, so\n" ); printf(" there is no record of any added Steiner points.\n"); printf(" -O No holes. Ignores the holes in the .poly file.\n"); printf( " -X No exact arithmetic. Normally, Triangle uses exact floating-point\n" ); printf( " arithmetic for certain tests if it thinks the inexact tests are not\n" ); printf( " accurate enough. Exact arithmetic ensures the robustness of the\n"); printf( " triangulation algorithms, despite floating-point roundoff error.\n"); printf( " Disabling exact arithmetic with the -X switch causes a small\n"); printf( " improvement in speed and creates the possibility (albeit small)\n"); printf( " that Triangle will fail to produce a valid mesh. Not recommended.\n" ); printf( " -z Numbers all items starting from zero (rather than one). Note that\n" ); printf( " this switch is normally overrided by the value used to number the\n"); printf( " first vertex of the input .node or .poly file. However, this\n"); printf( " switch is useful when calling Triangle from another program.\n"); printf( " -o2 Generates second-order subparametric elements with six nodes each.\n" ); printf( " -Y No new vertices on the boundary. This switch is useful when the\n"); printf( " mesh boundary must be preserved so that it conforms to some\n"); printf( " adjacent mesh. Be forewarned that you will probably sacrifice some\n" ); printf( " of the quality of the mesh; Triangle will try, but the resulting\n"); printf( " mesh may contain triangles of poor aspect ratio. Works well if all\n" ); printf( " the boundary vertices are closely spaced. Specify this switch\n"); printf( " twice (`-YY') to prevent all segment splitting, including internal\n" ); printf(" boundaries.\n"); printf( " -S Specifies the maximum number of Steiner points (vertices that are\n"); printf( " not in the input, but are added to meet the constraints on minimum\n" ); printf( " angle and maximum area). The default is to allow an unlimited\n"); printf( " number. If you specify this switch with no number after it,\n"); printf( " the limit is set to zero. Triangle always adds vertices at segment\n" ); printf( " intersections, even if it needs to use more vertices than the limit\n" ); printf( " you set. When Triangle inserts segments by splitting (-s), it\n"); printf( " always adds enough vertices to ensure that all the segments of the\n" ); printf(" PLSG are recovered, ignoring the limit if necessary.\n"); printf( " -L Do not use diametral lenses to determine whether subsegments are\n"); printf( " encroached; use diametral circles instead (as in Ruppert's\n"); printf( " algorithm). Use this switch if you want all triangles in the mesh\n" ); printf( " to be Delaunay, and not just constrained Delaunay; or if you want\n"); printf( " to ensure that all Voronoi vertices lie within the triangulation.\n"); printf( " (Applications such as some finite volume methods may have this\n"); printf( " requirement.) This switch may increase the number of vertices in\n"); printf(" the mesh to meet these constraints.\n"); printf( " -i Uses an incremental rather than divide-and-conquer algorithm to\n"); printf( " form a Delaunay triangulation. Try it if the divide-and-conquer\n"); printf(" algorithm fails.\n"); printf( " -F Uses Steven Fortune's sweepline algorithm to form a Delaunay\n"); printf( " triangulation. Warning: does not use exact arithmetic for all\n"); printf(" calculations. An exact result is not guaranteed.\n"); printf( " -l Uses only vertical cuts in the divide-and-conquer algorithm. By\n"); printf( " default, Triangle uses alternating vertical and horizontal cuts,\n"); printf( " which usually improve the speed except with vertex sets that are\n"); printf( " small or short and wide. This switch is primarily of theoretical\n"); printf(" interest.\n"); printf( " -s Specifies that segments should be forced into the triangulation by\n" ); printf( " recursively splitting them at their midpoints, rather than by\n"); printf( " generating a constrained Delaunay triangulation. Segment splitting\n" ); printf( " is true to Ruppert's original algorithm, but can create needlessly\n" ); printf( " small triangles. This switch is primarily of theoretical interest.\n" ); printf( " -C Check the consistency of the final mesh. Uses exact arithmetic for\n" ); printf( " checking, even if the -X switch is used. Useful if you suspect\n"); printf(" Triangle is buggy.\n"); printf( " -Q Quiet: Suppresses all explanation of what Triangle is doing,\n"); printf(" unless an error occurs.\n"); printf( " -V Verbose: Gives detailed information about what Triangle is doing.\n" ); printf( " Add more `V's for increasing amount of detail. `-V' gives\n"); printf( " information on algorithmic progress and more detailed statistics.\n"); printf( " `-VV' gives vertex-by-vertex details, and prints so much that\n"); printf( " Triangle runs much more slowly. `-VVVV' gives information only\n"); printf(" a debugger could love.\n"); printf(" -h Help: Displays these instructions.\n"); printf("\n"); printf("Definitions:\n"); printf("\n"); printf( " A Delaunay triangulation of a vertex set is a triangulation whose\n"); printf( " vertices are the vertex set, wherein no vertex in the vertex set falls in\n" ); printf( " the interior of the circumcircle (circle that passes through all three\n"); printf(" vertices) of any triangle in the triangulation.\n\n"); printf( " A Voronoi diagram of a vertex set is a subdivision of the plane into\n"); printf( " polygonal regions (some of which may be infinite), where each region is\n"); printf( " the set of points in the plane that are closer to some input vertex than\n" ); printf( " to any other input vertex. (The Voronoi diagram is the geometric dual of\n" ); printf(" the Delaunay triangulation.)\n\n"); printf( " A Planar Straight Line Graph (PSLG) is a set of vertices and segments.\n"); printf( " Segments are simply edges, whose endpoints are vertices in the PSLG.\n"); printf( " Segments may intersect each other only at their endpoints. The file\n"); printf(" format for PSLGs (.poly files) is described below.\n\n"); printf( " A constrained Delaunay triangulation (CDT) of a PSLG is similar to a\n"); printf( " Delaunay triangulation, but each PSLG segment is present as a single edge\n" ); printf( " in the triangulation. (A constrained Delaunay triangulation is not truly\n" ); printf( " a Delaunay triangulation.) By definition, a CDT does not have any\n"); printf(" vertices other than those specified in the input PSLG.\n\n"); printf( " A conforming Delaunay triangulation of a PSLG is a true Delaunay\n"); printf( " triangulation in which each PSLG segment is represented by a linear\n"); printf( " contiguous sequence of edges in the triangulation. Each input segment\n"); printf( " may have been subdivided into shorter subsegments by the insertion of\n"); printf( " additional vertices. These inserted vertices are necessary to maintain\n"); printf( " the Delaunay property while ensuring that every segment is represented.\n"); printf("\n"); printf("File Formats:\n"); printf("\n"); printf( " All files may contain comments prefixed by the character '#'. Vertices,\n" ); printf( " triangles, edges, holes, and maximum area constraints must be numbered\n"); printf( " consecutively, starting from either 1 or 0. Whichever you choose, all\n"); printf( " input files must be consistent; if the vertices are numbered from 1, so\n"); printf( " must be all other objects. Triangle automatically detects your choice\n"); printf( " while reading the .node (or .poly) file. (When calling Triangle from\n"); printf( " another program, use the -z switch if you wish to number objects from\n"); printf(" zero.) Examples of these file formats are given below.\n\n"); printf(" .node files:\n"); printf( " First line: <# of vertices> <# of attributes>\n" ); printf( " <# of boundary markers (0 or 1)>\n" ); printf( " Remaining lines: [attributes] [boundary marker]\n"); printf("\n"); printf( " The attributes, which are typically floating-point values of physical\n"); printf( " quantities (such as mass or conductivity) associated with the nodes of\n" ); printf( " a finite element mesh, are copied unchanged to the output mesh. If -q,\n" ); printf( " -a, -u, or -s is selected, each new Steiner point added to the mesh\n"); printf(" has attributes assigned to it by linear interpolation.\n\n"); printf( " If the fourth entry of the first line is `1', the last column of the\n"); printf( " remainder of the file is assumed to contain boundary markers. Boundary\n" ); printf( " markers are used to identify boundary vertices and vertices resting on\n" ); printf( " PSLG segments; a complete description appears in a section below. The\n" ); printf( " .node file produced by Triangle contains boundary markers in the last\n"); printf(" column unless they are suppressed by the -B switch.\n\n"); printf(" .ele files:\n"); printf( " First line: <# of triangles> <# of attributes>\n"); printf( " Remaining lines: ... [attributes]\n"); printf("\n"); printf( " Nodes are indices into the corresponding .node file. The first three\n"); printf( " nodes are the corner vertices, and are listed in counterclockwise order\n" ); printf( " around each triangle. (The remaining nodes, if any, depend on the type\n" ); printf(" of finite element used.)\n\n"); printf( " The attributes are just like those of .node files. Because there is no\n" ); printf( " simple mapping from input to output triangles, an attempt is made to\n"); printf( " interpolate attributes, which may result in a good deal of diffusion of\n" ); printf( " attributes among nearby triangles as the triangulation is refined.\n"); printf( " Attributes do not diffuse across segments, so attributes used to\n"); printf(" identify segment-bounded regions remain intact.\n\n"); printf( " In .ele files produced by Triangle, each triangular element has three\n"); printf( " nodes (vertices) unless the -o2 switch is used, in which case\n"); printf( " subparametric quadratic elements with six nodes each are generated.\n"); printf( " The first three nodes are the corners in counterclockwise order, and\n"); printf( " the fourth, fifth, and sixth nodes lie on the midpoints of the edges\n"); printf( " opposite the first, second, and third vertices, respectively.\n"); printf("\n"); printf(" .poly files:\n"); printf( " First line: <# of vertices> <# of attributes>\n" ); printf( " <# of boundary markers (0 or 1)>\n" ); printf( " Following lines: [attributes] [boundary marker]\n"); printf(" One line: <# of segments> <# of boundary markers (0 or 1)>\n"); printf( " Following lines: [boundary marker]\n"); printf(" One line: <# of holes>\n"); printf(" Following lines: \n"); printf( " Optional line: <# of regional attributes and/or area constraints>\n"); printf( " Optional following lines: \n"); printf("\n"); printf( " A .poly file represents a PSLG, as well as some additional information.\n" ); printf( " The first section lists all the vertices, and is identical to the\n"); printf( " format of .node files. <# of vertices> may be set to zero to indicate\n" ); printf( " that the vertices are listed in a separate .node file; .poly files\n"); printf( " produced by Triangle always have this format. A vertex set represented\n" ); printf( " this way has the advantage that it may easily be triangulated with or\n"); printf( " without segments (depending on whether the .poly or .node file is\n"); printf(" read).\n\n"); printf( " The second section lists the segments. Segments are edges whose\n"); printf( " presence in the triangulation is enforced (although each segment may be\n" ); printf( " subdivided into smaller edges). Each segment is specified by listing\n"); printf( " the indices of its two endpoints. This means that you must include its\n" ); printf( " endpoints in the vertex list. Each segment, like each point, may have\n" ); printf(" a boundary marker.\n\n"); printf( " If -q, -a, -u, and -s are not selected, Triangle produces a constrained\n" ); printf( " Delaunay triangulation (CDT), in which each segment appears as a single\n" ); printf( " edge in the triangulation. If -q, -a, -u, or -s is selected, Triangle\n" ); printf( " produces a conforming constrained Delaunay triangulation (CCDT), in\n"); printf( " which segments may be subdivided into smaller edges. If -L is selected\n" ); printf( " as well, Triangle produces a conforming Delaunay triangulation, so\n"); printf( " every triangle is Delaunay, and not just constrained Delaunay.\n"); printf("\n"); printf( " The third section lists holes (and concavities, if -c is selected) in\n"); printf( " the triangulation. Holes are specified by identifying a point inside\n"); printf( " each hole. After the triangulation is formed, Triangle creates holes\n"); printf( " by eating triangles, spreading out from each hole point until its\n"); printf( " progress is blocked by PSLG segments; you must be careful to enclose\n"); printf( " each hole in segments, or your whole triangulation might be eaten away.\n" ); printf( " If the two triangles abutting a segment are eaten, the segment itself\n"); printf( " is also eaten. Do not place a hole directly on a segment; if you do,\n"); printf(" Triangle chooses one side of the segment arbitrarily.\n\n"); printf( " The optional fourth section lists regional attributes (to be assigned\n"); printf( " to all triangles in a region) and regional constraints on the maximum\n"); printf( " triangle area. Triangle reads this section only if the -A switch is\n"); printf( " used or the -a switch is used without a number following it, and the -r\n" ); printf( " switch is not used. Regional attributes and area constraints are\n"); printf( " propagated in the same manner as holes; you specify a point for each\n"); printf( " attribute and/or constraint, and the attribute and/or constraint\n"); printf( " affects the whole region (bounded by segments) containing the point.\n"); printf( " If two values are written on a line after the x and y coordinate, the\n"); printf( " first such value is assumed to be a regional attribute (but is only\n"); printf( " applied if the -A switch is selected), and the second value is assumed\n" ); printf( " to be a regional area constraint (but is only applied if the -a switch\n" ); printf( " is selected). You may specify just one value after the coordinates,\n"); printf( " which can serve as both an attribute and an area constraint, depending\n" ); printf( " on the choice of switches. If you are using the -A and -a switches\n"); printf( " simultaneously and wish to assign an attribute to some region without\n"); printf(" imposing an area constraint, use a negative maximum area.\n\n"); printf( " When a triangulation is created from a .poly file, you must either\n"); printf( " enclose the entire region to be triangulated in PSLG segments, or\n"); printf( " use the -c switch, which encloses the convex hull of the input vertex\n"); printf( " set. If you do not use the -c switch, Triangle eats all triangles that\n" ); printf( " are not enclosed by segments; if you are not careful, your whole\n"); printf( " triangulation may be eaten away. If you do use the -c switch, you can\n" ); printf( " still produce concavities by the appropriate placement of holes just\n"); printf(" within the convex hull.\n\n"); printf( " An ideal PSLG has no intersecting segments, nor any vertices that lie\n"); printf( " upon segments (except, of course, the endpoints of each segment.) You\n" ); printf( " aren't required to make your .poly files ideal, but you should be aware\n" ); printf( " of what can go wrong. Segment intersections are relatively safe--\n"); printf( " Triangle calculates the intersection points for you and adds them to\n"); printf( " the triangulation--as long as your machine's floating-point precision\n"); printf( " doesn't become a problem. You are tempting the fates if you have three\n" ); printf( " segments that cross at the same location, and expect Triangle to figure\n" ); printf( " out where the intersection point is. Thanks to floating-point roundoff\n" ); printf( " error, Triangle will probably decide that the three segments intersect\n" ); printf( " at three different points, and you will find a minuscule triangle in\n"); printf( " your output--unless Triangle tries to refine the tiny triangle, uses\n"); printf( " up the last bit of machine precision, and fails to terminate at all.\n"); printf( " You're better off putting the intersection point in the input files,\n"); printf( " and manually breaking up each segment into two. Similarly, if you\n"); printf( " place a vertex at the middle of a segment, and hope that Triangle will\n" ); printf( " break up the segment at that vertex, you might get lucky. On the other\n" ); printf( " hand, Triangle might decide that the vertex doesn't lie precisely on\n"); printf( " the segment, and you'll have a needle-sharp triangle in your output--or\n" ); printf(" a lot of tiny triangles if you're generating a quality mesh.\n"); printf("\n"); printf( " When Triangle reads a .poly file, it also writes a .poly file, which\n"); printf( " includes all edges that are parts of input segments. If the -c switch\n" ); printf( " is used, the output .poly file also includes all of the edges on the\n"); printf( " convex hull. Hence, the output .poly file is useful for finding edges\n" ); printf( " associated with input segments and for setting boundary conditions in\n"); printf( " finite element simulations. Moreover, you will need it if you plan to\n" ); printf( " refine the output mesh, and don't want segments to be missing in later\n" ); printf(" triangulations.\n\n"); printf(" .area files:\n"); printf(" First line: <# of triangles>\n"); printf(" Following lines: \n\n"); printf( " An .area file associates with each triangle a maximum area that is used\n" ); printf( " for mesh refinement. As with other file formats, every triangle must\n"); printf( " be represented, and they must be numbered consecutively. A triangle\n"); printf( " may be left unconstrained by assigning it a negative maximum area.\n"); printf("\n"); printf(" .edge files:\n"); printf(" First line: <# of edges> <# of boundary markers (0 or 1)>\n"); printf( " Following lines: [boundary marker]\n"); printf("\n"); printf( " Endpoints are indices into the corresponding .node file. Triangle can\n" ); printf( " produce .edge files (use the -e switch), but cannot read them. The\n"); printf( " optional column of boundary markers is suppressed by the -B switch.\n"); printf("\n"); printf( " In Voronoi diagrams, one also finds a special kind of edge that is an\n"); printf( " infinite ray with only one endpoint. For these edges, a different\n"); printf(" format is used:\n\n"); printf(" -1 \n\n"); printf( " The `direction' is a floating-point vector that indicates the direction\n" ); printf(" of the infinite ray.\n\n"); printf(" .neigh files:\n"); printf( " First line: <# of triangles> <# of neighbors per triangle (always 3)>\n" ); printf( " Following lines: \n"); printf("\n"); printf( " Neighbors are indices into the corresponding .ele file. An index of -1\n" ); printf( " indicates no neighbor (because the triangle is on an exterior\n"); printf( " boundary). The first neighbor of triangle i is opposite the first\n"); printf(" corner of triangle i, and so on.\n\n"); printf( " Triangle can produce .neigh files (use the -n switch), but cannot read\n" ); printf(" them.\n\n"); printf("Boundary Markers:\n\n"); printf( " Boundary markers are tags used mainly to identify which output vertices\n"); printf( " and edges are associated with which PSLG segment, and to identify which\n"); printf( " vertices and edges occur on a boundary of the triangulation. A common\n"); printf( " use is to determine where boundary conditions should be applied to a\n"); printf( " finite element mesh. You can prevent boundary markers from being written\n" ); printf(" into files produced by Triangle by using the -B switch.\n\n"); printf( " The boundary marker associated with each segment in an output .poly file\n" ); printf(" and each edge in an output .edge file is chosen as follows:\n"); printf( " - If an output edge is part or all of a PSLG segment with a nonzero\n"); printf( " boundary marker, then the edge is assigned the same marker.\n"); printf( " - Otherwise, if the edge occurs on a boundary of the triangulation\n"); printf( " (including boundaries of holes), then the edge is assigned the marker\n" ); printf(" one (1).\n"); printf(" - Otherwise, the edge is assigned the marker zero (0).\n"); printf( " The boundary marker associated with each vertex in an output .node file\n"); printf(" is chosen as follows:\n"); printf( " - If a vertex is assigned a nonzero boundary marker in the input file,\n" ); printf( " then it is assigned the same marker in the output .node file.\n"); printf( " - Otherwise, if the vertex lies on a PSLG segment (including the\n"); printf( " segment's endpoints) with a nonzero boundary marker, then the vertex\n" ); printf( " is assigned the same marker. If the vertex lies on several such\n"); printf(" segments, one of the markers is chosen arbitrarily.\n"); printf( " - Otherwise, if the vertex occurs on a boundary of the triangulation,\n"); printf(" then the vertex is assigned the marker one (1).\n"); printf(" - Otherwise, the vertex is assigned the marker zero (0).\n"); printf("\n"); printf( " If you want Triangle to determine for you which vertices and edges are on\n" ); printf( " the boundary, assign them the boundary marker zero (or use no markers at\n" ); printf( " all) in your input files. In the output files, all boundary vertices,\n"); printf(" edges, and segments are assigned the value one.\n\n"); printf("Triangulation Iteration Numbers:\n\n"); printf( " Because Triangle can read and refine its own triangulations, input\n"); printf( " and output files have iteration numbers. For instance, Triangle might\n"); printf( " read the files mesh.3.node, mesh.3.ele, and mesh.3.poly, refine the\n"); printf( " triangulation, and output the files mesh.4.node, mesh.4.ele, and\n"); printf(" mesh.4.poly. Files with no iteration number are treated as if\n"); printf( " their iteration number is zero; hence, Triangle might read the file\n"); printf( " points.node, triangulate it, and produce the files points.1.node and\n"); printf(" points.1.ele.\n\n"); printf( " Iteration numbers allow you to create a sequence of successively finer\n"); printf( " meshes suitable for multigrid methods. They also allow you to produce a\n" ); printf( " sequence of meshes using error estimate-driven mesh refinement.\n"); printf("\n"); printf( " If you're not using refinement or quality meshing, and you don't like\n"); printf( " iteration numbers, use the -I switch to disable them. This switch also\n"); printf( " disables output of .node and .poly files to prevent your input files from\n" ); printf( " being overwritten. (If the input is a .poly file that contains its own\n"); printf(" points, a .node file is written.)\n\n"); printf("Examples of How to Use Triangle:\n\n"); printf( " `triangle dots' reads vertices from dots.node, and writes their Delaunay\n" ); printf( " triangulation to dots.1.node and dots.1.ele. (dots.1.node is identical\n"); printf( " to dots.node.) `triangle -I dots' writes the triangulation to dots.ele\n"); printf( " instead. (No additional .node file is needed, so none is written.)\n"); printf("\n"); printf( " `triangle -pe object.1' reads a PSLG from object.1.poly (and possibly\n"); printf( " object.1.node, if the vertices are omitted from object.1.poly) and writes\n" ); printf( " its constrained Delaunay triangulation to object.2.node and object.2.ele.\n" ); printf( " The segments are copied to object.2.poly, and all edges are written to\n"); printf(" object.2.edge.\n\n"); printf( " `triangle -pq31.5a.1 object' reads a PSLG from object.poly (and possibly\n" ); printf( " object.node), generates a mesh whose angles are all 31.5 degrees or\n"); printf( " greater and whose triangles all have areas of 0.1 or less, and writes the\n" ); printf( " mesh to object.1.node and object.1.ele. Each segment may be broken up\n"); printf(" into multiple subsegments; these are written to object.1.poly.\n"); printf("\n"); printf( " Here is a sample file `box.poly' describing a square with a square hole:\n" ); printf("\n"); printf( " # A box with eight vertices in 2D, no attributes, one boundary marker.\n" ); printf(" 8 2 0 1\n"); printf(" # Outer box has these vertices:\n"); printf(" 1 0 0 0\n"); printf(" 2 0 3 0\n"); printf(" 3 3 0 0\n"); printf(" 4 3 3 33 # A special marker for this vertex.\n"); printf(" # Inner square has these vertices:\n"); printf(" 5 1 1 0\n"); printf(" 6 1 2 0\n"); printf(" 7 2 1 0\n"); printf(" 8 2 2 0\n"); printf(" # Five segments with boundary markers.\n"); printf(" 5 1\n"); printf(" 1 1 2 5 # Left side of outer box.\n"); printf(" # Square hole has these segments:\n"); printf(" 2 5 7 0\n"); printf(" 3 7 8 0\n"); printf(" 4 8 6 10\n"); printf(" 5 6 5 0\n"); printf(" # One hole in the middle of the inner square.\n"); printf(" 1\n"); printf(" 1 1.5 1.5\n"); printf("\n"); printf( " Note that some segments are missing from the outer square, so one must\n"); printf( " use the `-c' switch. After `triangle -pqc box.poly', here is the output\n" ); printf( " file `box.1.node', with twelve vertices. The last four vertices were\n"); printf( " added to meet the angle constraint. Vertices 1, 2, and 9 have markers\n"); printf( " from segment 1. Vertices 6 and 8 have markers from segment 4. All the\n"); printf( " other vertices but 4 have been marked to indicate that they lie on a\n"); printf(" boundary.\n\n"); printf(" 12 2 0 1\n"); printf(" 1 0 0 5\n"); printf(" 2 0 3 5\n"); printf(" 3 3 0 1\n"); printf(" 4 3 3 33\n"); printf(" 5 1 1 1\n"); printf(" 6 1 2 10\n"); printf(" 7 2 1 1\n"); printf(" 8 2 2 10\n"); printf(" 9 0 1.5 5\n"); printf(" 10 1.5 0 1\n"); printf(" 11 3 1.5 1\n"); printf(" 12 1.5 3 1\n"); printf(" # Generated by triangle -pqc box.poly\n"); printf("\n"); printf(" Here is the output file `box.1.ele', with twelve triangles.\n"); printf("\n"); printf(" 12 3 0\n"); printf(" 1 5 6 9\n"); printf(" 2 10 3 7\n"); printf(" 3 6 8 12\n"); printf(" 4 9 1 5\n"); printf(" 5 6 2 9\n"); printf(" 6 7 3 11\n"); printf(" 7 11 4 8\n"); printf(" 8 7 5 10\n"); printf(" 9 12 2 6\n"); printf(" 10 8 7 11\n"); printf(" 11 5 1 10\n"); printf(" 12 8 4 12\n"); printf(" # Generated by triangle -pqc box.poly\n\n"); printf( " Here is the output file `box.1.poly'. Note that segments have been added\n" ); printf( " to represent the convex hull, and some segments have been split by newly\n" ); printf( " added vertices. Note also that <# of vertices> is set to zero to\n"); printf(" indicate that the vertices should be read from the .node file.\n"); printf("\n"); printf(" 0 2 0 1\n"); printf(" 12 1\n"); printf(" 1 1 9 5\n"); printf(" 2 5 7 1\n"); printf(" 3 8 7 1\n"); printf(" 4 6 8 10\n"); printf(" 5 5 6 1\n"); printf(" 6 3 10 1\n"); printf(" 7 4 11 1\n"); printf(" 8 2 12 1\n"); printf(" 9 9 2 5\n"); printf(" 10 10 1 1\n"); printf(" 11 11 3 1\n"); printf(" 12 12 4 1\n"); printf(" 1\n"); printf(" 1 1.5 1.5\n"); printf(" # Generated by triangle -pqc box.poly\n"); printf("\n"); printf("Refinement and Area Constraints:\n"); printf("\n"); printf( " The -r switch causes a mesh (.node and .ele files) to be read and\n"); printf( " refined. If the -p switch is also used, a .poly file is read and used to\n" ); printf( " specify edges that are constrained and cannot be eliminated (although\n"); printf( " they can be divided into smaller edges) by the refinement process.\n"); printf("\n"); printf( " When you refine a mesh, you generally want to impose tighter quality\n"); printf( " constraints. One way to accomplish this is to use -q with a larger\n"); printf( " angle, or -a followed by a smaller area than you used to generate the\n"); printf( " mesh you are refining. Another way to do this is to create an .area\n"); printf( " file, which specifies a maximum area for each triangle, and use the -a\n"); printf( " switch (without a number following). Each triangle's area constraint is\n" ); printf( " applied to that triangle. Area constraints tend to diffuse as the mesh\n"); printf( " is refined, so if there are large variations in area constraint between\n"); printf(" adjacent triangles, you may not get the results you want.\n\n"); printf( " If you are refining a mesh composed of linear (three-node) elements, the\n" ); printf( " output mesh contains all the nodes present in the input mesh, in the same\n" ); printf( " order, with new nodes added at the end of the .node file. However, the\n"); printf( " refinement is not hierarchical: there is no guarantee that each output\n"); printf( " element is contained in a single input element. Often, output elements\n"); printf( " overlap two input elements, and some input edges are not present in the\n"); printf( " output mesh. Hence, a sequence of refined meshes forms a hierarchy of\n"); printf( " nodes, but not a hierarchy of elements. If you refine a mesh of higher-\n" ); printf( " order elements, the hierarchical property applies only to the nodes at\n"); printf( " the corners of an element; other nodes may not be present in the refined\n" ); printf(" mesh.\n\n"); printf( " Maximum area constraints in .poly files operate differently from those in\n" ); printf( " .area files. A maximum area in a .poly file applies to the whole\n"); printf( " (segment-bounded) region in which a point falls, whereas a maximum area\n"); printf( " in an .area file applies to only one triangle. Area constraints in .poly\n" ); printf( " files are used only when a mesh is first generated, whereas area\n"); printf( " constraints in .area files are used only to refine an existing mesh, and\n" ); printf( " are typically based on a posteriori error estimates resulting from a\n"); printf(" finite element simulation on that mesh.\n\n"); printf( " `triangle -rq25 object.1' reads object.1.node and object.1.ele, then\n"); printf( " refines the triangulation to enforce a 25 degree minimum angle, and then\n" ); printf( " writes the refined triangulation to object.2.node and object.2.ele.\n"); printf("\n"); printf( " `triangle -rpaa6.2 z.3' reads z.3.node, z.3.ele, z.3.poly, and z.3.area.\n" ); printf( " After reconstructing the mesh and its subsegments, Triangle refines the\n"); printf( " mesh so that no triangle has area greater than 6.2, and furthermore the\n"); printf( " triangles satisfy the maximum area constraints in z.3.area. No angle\n"); printf( " bound is imposed at all. The output is written to z.4.node, z.4.ele, and\n" ); printf(" z.4.poly.\n\n"); printf( " The sequence `triangle -qa1 x', `triangle -rqa.3 x.1', `triangle -rqa.1\n"); printf( " x.2' creates a sequence of successively finer