triangle.cpp

一个2D电磁场FEM计算的VC++源程序

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;

long incirclecount;                   /* Number of incircle tests performed. */
long counterclockcount;       /* Number of counterclockwise 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. */

/* 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.               */
/*   regionattrib: -A switch.  convex: -c 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.   */
/*     steinerleft: number of Steiner points not yet used.                   */
/*   incremental: -i switch.  sweepline: -F switch.                          */
/*   dwyer: inverse of -l switch.                                            */
/*   splitseg: -s switch.                                                    */
/*   docheck: -C switch.                                                     */
/*   quiet: -Q switch.  verbose: count of how often -V switch is selected.   */
/*   useshelles: -p, -r, -q, or -c switch; determines whether shell edges    */
/*     are used at all.                                                      */
/*                                                                           */
/* Read the instructions to find out the meaning of these switches.          */

int poly, refine, quality, vararea, fixedarea, regionattrib, convex;
int firstnumber;
int edgesout, voronoi, neighbors, geomview;
int nobound, nopolywritten, nonodewritten, noelewritten, noiterationnum;
int noholes, noexact;
int incremental, sweepline, dwyer;
int splitseg;
int docheck;
int quiet, verbose;
int useshelles;
int order;
int nobisect;
int steiner, steinerleft;
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 */

/* Triangular bounding box points.                                           */

point infpoint1, infpoint2, infpoint3;

/* 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 shell edge.  Referenced by any triangle or     */
/*   shell edge that isn't really connected to a shell edge at that          */
/*   location.                                                               */

shelle *dummysh;
shelle *dummyshbase;         /* Keep base address so we can free() it later. */

/* Pointer to a recently visited triangle.  Improves point location if       */
/*   proximate points are inserted sequentially.                             */

struct triedge recenttri;

/*****************************************************************************/
/*                                                                           */
/*  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.*/
/*                                                                           */
/*  Shell edges are manipulated similarly.  A pointer to a shell edge        */
/*  carries both an address and an orientation in the range 0 to 1.          */
/*                                                                           */
/*  The other primitives take an oriented triangle or oriented shell edge,   */
/*  and return an oriented triangle or oriented shell edge or point; or they */
/*  change the connections in the data structure.                            */
/*                                                                           */
/*****************************************************************************/

/********* 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, triedge)                                                  \
  (triedge).orient = (int) ((unsigned long) (ptr) & (unsigned long) 3l);      \
  (triedge).tri = (triangle *)                                                \
                  ((unsigned long) (ptr) ^ (unsigned long) (triedge).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 (triedge).tri are zero.*/

#define encode(triedge)                                                       \
  (triangle) ((unsigned long) (triedge).tri | (unsigned long) (triedge).orient)

/* The following edge manipulation primitives are all described by Guibas    */
/*   and Stolfi.  However, they use an edge-based data structure, whereas I  */
/*   am using a triangle-based data structure.                               */

/* sym() finds the abutting triangle, on the same edge.  Note that the       */
/*   edge direction is necessarily reversed, because triangle/edge handles   */
/*   are always directed counterclockwise around the triangle.               */

#define sym(triedge1, triedge2)                                               \
  ptr = (triedge1).tri[(triedge1).orient];                                    \
  decode(ptr, triedge2);

#define symself(triedge)                                                      \
  ptr = (triedge).tri[(triedge).orient];                                      \
  decode(ptr, triedge);

/* lnext() finds the next edge (counterclockwise) of a triangle.             */

#define lnext(triedge1, triedge2)                                             \
  (triedge2).tri = (triedge1).tri;                                            \
  (triedge2).orient = plus1mod3[(triedge1).orient]

#define lnextself(triedge)                                                    \
  (triedge).orient = plus1mod3[(triedge).orient]

/* lprev() finds the previous edge (clockwise) of a triangle.                */

#define lprev(triedge1, triedge2)                                             \
  (triedge2).tri = (triedge1).tri;                                            \
  (triedge2).orient = minus1mod3[(triedge1).orient]

#define lprevself(triedge)                                                    \
  (triedge).orient = minus1mod3[(triedge).orient]

/* onext() spins counterclockwise around a point; that is, it finds the next */
/*   edge with the same origin in the counterclockwise direction.  This edge */
/*   will be part of a different triangle.                                   */

#define onext(triedge1, triedge2)                                             \
  lprev(triedge1, triedge2);                                                  \
  symself(triedge2);

#define onextself(triedge)                                                    \
  lprevself(triedge);                                                         \
  symself(triedge);

/* oprev() spins clockwise around a point; that is, it finds the next edge   */
/*   with the same origin in the clockwise direction.  This edge will be     */
/*   part of a different triangle.                                           */

#define oprev(triedge1, triedge2)                                             \
  sym(triedge1, triedge2);                                                    \
  lnextself(triedge2);

#define oprevself(triedge)                                                    \
  symself(triedge);                                                           \
  lnextself(triedge);

/* dnext() spins counterclockwise around a point; that is, it finds the next */
/*   edge with the same destination in the counterclockwise direction.  This */
/*   edge will be part of a different triangle.                              */

#define dnext(triedge1, triedge2)                                             \
  sym(triedge1, triedge2);                                                    \
  lprevself(triedge2);

#define dnextself(triedge)                                                    \
  symself(triedge);                                                           \
  lprevself(triedge);

/* dprev() spins clockwise around a point; that is, it finds the next edge   */
/*   with the same destination in the clockwise direction.  This edge will   */