vtm/vtm-android/jni/triangle/triangle_private.h
Hannes Janetzek 83cd73156a split up
2013-10-09 01:56:08 +02:00

1068 lines
57 KiB
C

/* 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, -u, -D, -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 volatile" below. Normally, however, */
/* should be defined to be nothing. ("#define ".) */
/* */
/* 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 /* Nothing */
/* #define volatile */
/* Maximum number of characters in a file name (including the null). */
#define FILENAMESIZE 2048
/* Maximum number of characters in a line read from a file (including the */
/* null). */
#define INPUTLINESIZE 1024
/* 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
#define VOID void
/* 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. */
/* 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 <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#ifndef NO_TIMER
#include <sys/time.h>
#endif /* not NO_TIMER */
#ifdef CPU86
#include <float.h>
#endif /* CPU86 */
#ifdef LINUX
#include <fpu_control.h>
#endif /* LINUX */
#include "triangle.h"
#include <android/log.h>
#define printf(...) __android_log_print(ANDROID_LOG_DEBUG, "Triangle", __VA_ARGS__)
/* A few forward declarations. */
/* 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 four vertices--the vertices of the */
/* subsegment, and the vertices of the segment it is a part of--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 four pointers to vertices, plus two */
/* pointers to adjoining triangles, plus one boundary marker, plus one */
/* segment number. */
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. */
/* */
/* alignbytes determines how new records should be aligned in memory. */
/* itembytes is the length of a record in bytes (after rounding up). */
/* itemsperblock is the number of items allocated at once in a single */
/* block. itemsfirstblock is the number of items in the first block, */
/* which can vary from the others. 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;
int alignbytes;
int itembytes;
int itemsperblock;
int itemsfirstblock;
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 4095 (highest priority) to 0 (lowest priority). */
struct badtriang *queuefront[4096];
struct badtriang *queuetail[4096];
int nextnonemptyq[4096];
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'. */
/*****************************************************************************/
/* */
/* 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. */
extern int plus1mod3[]; // = { 1, 2, 0 };
extern int minus1mod3[]; // = { 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 or the segment that includes it. */
#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
#define segorg(osub, vertexptr) \
vertexptr = (vertex) (osub).ss[4 + (osub).ssorient]
#define segdest(osub, vertexptr) \
vertexptr = (vertex) (osub).ss[5 - (osub).ssorient]
#define setsegorg(osub, vertexptr) \
(osub).ss[4 + (osub).ssorient] = (subseg) vertexptr
#define setsegdest(osub, vertexptr) \
(osub).ss[5 - (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 + 8))
#define setmark(osub, value) \
* (int *) ((osub).ss + 8) = 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[6 + (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[6 + (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[6 + (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 *********/
void printtriangle(struct mesh *m, struct behavior *b, struct otri *t);
void printsubseg(struct mesh *m, struct behavior *b, struct osub *s);
void quality_statistics(struct mesh *m, struct behavior *b);
void statistics(struct mesh *m, struct behavior *b);
// provided by triangle.c
void traversalinit(struct memorypool *pool);
REAL counterclockwise(struct mesh *m, struct behavior *b, vertex pa, vertex pb, vertex pc);
triangle *triangletraverse(struct mesh *m);