triangle.cpp

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

/*  at the beginning of the segment insertion phase, and are not used or     */
/*  updated at any other time.  Edge swapping during segment insertion will  */
/*  render some of them incorrect.  Hence, don't rely upon them for          */
/*  anything.  For the most part, points do not have any information about   */
/*  what triangles or shell edges they are linked to.                        */
/*                                                                           */
/*****************************************************************************/

/*****************************************************************************/
/*                                                                           */
/*  Handles                                                                  */
/*                                                                           */
/*  The oriented triangle (`triedge') and oriented shell edge (`edge') data  */
/*  structures defined below do not themselves store any part of the mesh.   */
/*  The mesh itself is made of `triangle's, `shelle's, and `point's.         */
/*                                                                           */
/*  Oriented triangles and oriented shell edges 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.)                               */
/*                                                                           */
/*  A `triedge' is a handle that holds a triangle.  It holds a specific side */
/*  of the triangle.  An `edge' is a handle that holds a shell edge.  It     */
/*  holds either the left or right side of the edge.                         */
/*                                                                           */
/*  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 exact position of the handles is        */
/*  important.  For instance, when two triangles are glued together by the   */
/*  bond() primitive, they are glued by the sides on which the handles lie.  */
/*                                                                           */
/*  Because points have no information about which triangles they are        */
/*  attached to, I commonly represent a point by use of a handle whose       */
/*  origin is the point.  A single handle can simultaneously represent a     */
/*  triangle, an edge, and a point.                                          */
/*                                                                           */
/*****************************************************************************/

/* The triangle data structure.  Each triangle contains three pointers to    */
/*   adjoining triangles, plus three pointers to vertex points, plus three   */
/*   pointers to shell edges (defined below; these pointers are usually      */
/*   `dummysh').  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 defined the type `triangle' to be 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 is always        */
/*   directed to point counterclockwise about the corresponding triangle.    */

struct triedge {
  triangle *tri;
  int orient;                                         /* Ranges from 0 to 2. */
};

/* The shell data structure.  Each shell edge contains two pointers to       */
/*   adjoining shell edges, plus two pointers to vertex points, plus two     */
/*   pointers to adjoining triangles, plus one shell marker.                 */

typedef REAL **shelle;                  /* Really:  typedef shelle *shelle   */

/* An oriented shell edge:  includes a pointer to a shell edge 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 edge {
  shelle *sh;
  int shorient;                                       /* Ranges from 0 to 1. */
};

/* The point data structure.  Each point 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 *point;

/* A queue used to store encroached segments.  Each segment's vertices are   */
/*   stored so that one can check whether a segment is still the same.       */

struct badsegment {
  struct edge encsegment;                          /* An encroached segment. */
  point segorg, segdest;                                /* The two vertices. */
  struct badsegment *nextsegment;     /* Pointer to next encroached segment. */
};

/* 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 badface {
  struct triedge badfacetri;                              /* A bad triangle. */
  REAL key;                             /* cos^2 of smallest (apical) angle. */
  point faceorg, facedest, faceapex;                  /* The three vertices. */
  struct badface *nextface;                 /* Pointer to next bad triangle. */
};

/* 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 `point' 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 point 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 triedge keyedge;                  /* Lprev of an edge on the front. */
  point 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;
};

/* Variables used to allocate memory for triangles, shell edges, points,     */
/*   viri (triangles being eaten), bad (encroached) segments, bad (skinny    */
/*   or too large) triangles, and splay tree nodes.                          */

struct memorypool triangles;
struct memorypool shelles;
struct memorypool points;
struct memorypool viri;
struct memorypool badsegments;
struct memorypool badtriangles;
struct memorypool splaynodes;

/* Variables that maintain the bad triangle queues.  The tails are pointers  */
/*   to the pointers that have to be filled in to enqueue an item.           */

struct badface *queuefront[64];
struct badface **queuetail[64];

REAL xmin, xmax, ymin, ymax;                              /* x and y bounds. */
REAL xminextreme;        /* Nonexistent x value used as a flag in sweepline. */
int inpoints;                                     /* Number of input points. */
int inelements;                                /* Number of input triangles. */
int insegments;                                 /* Number of input segments. */
int holes;                                         /* Number of input holes. */
int regions;                                     /* Number of input regions. */
long edges;                                       /* Number of output edges. */
int mesh_dim;                                  /* Dimension (ought to be 2). */
int nextras;                              /* Number of attributes per point. */
int eextras;                           /* Number of attributes per triangle. */
long hullsize;                            /* Number of edges of convex hull. */
int triwords;                                   /* Total words per triangle. */
int shwords;                                  /* Total words per shell edge. */
int pointmarkindex;             /* Index to find boundary marker of a point. */
int point2triindex;         /* Index to find a triangle adjacent to a point. */
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 readnodefile;                             /* Has a .node file been read? */
long samples;                /* Number of random samples for point location. */
unsigned long randomseed;                     /* Current random number seed. */