tetgenBR.h

Go to the documentation of this file.

// This file is part of Tetgen/BR, the Boundary Recovery code of TetGen
//
// Copyright 2015-2016 Weierstrass Institute for Applied Analysis and
// Stochastics
//
// This file is relicensed under the terms of LICENSE.txt for use in Gmsh thanks
// to a Software License Agreement between Weierstrass Institute for Applied
// Analysis and Stochastics and GMESH SPRL.
 
//                                                                           //
// TetGen                                                                    //
//                                                                           //
// A Quality Tetrahedral Mesh Generator and A 3D Delaunay Triangulator       //
//                                                                           //
// Version 1.5                                                               //
// May 31, 2014                                                              //
//                                                                           //
// Copyright (C) 2002--2016                                                  //
//                                                                           //
// TetGen is freely available through the website: http://www.tetgen.org.    //
//   It may be copied, modified, and redistributed for non-commercial use.   //
//   Please consult the file LICENSE for the detailed copyright notices.     //
//                                                                           //
 
// The following code of this file was automatically generated.  Do not edit!
// TetGenBR -- The Boundary Recovery code of TetGen.
 
#ifndef tetgenBRH
#define tetgenBRH
 
//                                                                           //
// tetgenbehavior                                                            //
//                                                                           //
// A structure for maintaining the switches and parameters used by TetGen's  //
// mesh data structure and algorithms.                                       //
//                                                                           //
// All switches and parameters are initialized with default values. They can //
// be set by the command line arguments (a list of strings) of TetGen.       //
//                                                                           //
// NOTE: Some of the switches are incompatible. While some may depend on     //
// other switches.  The routine parse_commandline() sets the switches from   //
// the command line (a list of strings) and checks the consistency of the    //
// applied switches.                                                         //
//                                                                           //
 
class tetgenbehavior {
public:
  // Switches of TetGen.
  int plc; // '-p', 0.
  int psc; // '-s', 0.
  int refine; // '-r', 0.
  int quality; // '-q', 0.
  int nobisect; // '-Y', 0.
  int coarsen; // '-R', 0.
  int weighted; // '-w', 0.
  int brio_hilbert; // '-b', 1.
  int incrflip; // '-l', 0.
  int flipinsert; // '-L', 0.
  int metric; // '-m', 0.
  int varvolume; // '-a', 0.
  int fixedvolume; // '-a', 0.
  int regionattrib; // '-A', 0.
  int cdtrefine; // '-D', 0.
  int insertaddpoints; // '-i', 0.
  int diagnose; // '-d', 0.
  int convex; // '-c', 0.
  int nomergefacet; // '-M', 0.
  int nomergevertex; // '-M', 0.
  int noexact; // '-X', 0.
  int nostaticfilter; // '-X', 0.
  int zeroindex; // '-z', 0.
  int facesout; // '-f', 0.
  int edgesout; // '-e', 0.
  int neighout; // '-n', 0.
  int voroout; // '-v', 0.
  int meditview; // '-g', 0.
  int vtkview; // '-k', 0.
  int nobound; // '-B', 0.
  int nonodewritten; // '-N', 0.
  int noelewritten; // '-E', 0.
  int nofacewritten; // '-F', 0.
  int noiterationnum; // '-I', 0.
  int nojettison; // '-J', 0.
  int docheck; // '-C', 0.
  int quiet; // '-Q', 0.
  int verbose; // '-V', 0.
 
  // Parameters of TetGen.
  int vertexperblock; // '-x', 4092.
  int tetrahedraperblock; // '-x', 8188.
  int shellfaceperblock; // '-x', 2044.
  int nobisect_nomerge; // '-Y', 1.
  int supsteiner_level; // '-Y/', 2.
  int addsteiner_algo; // '-Y//', 1.
  int coarsen_param; // '-R', 0.
  int weighted_param; // '-w', 0.
  int fliplinklevel; // -1.
  int flipstarsize; // -1.
  int fliplinklevelinc; //  1.
  int reflevel; // '-D', 3.
  int optlevel; // '-O', 2.
  int optscheme; // '-O', 7.
  int delmaxfliplevel; // 1.
  int order; // '-o', 1.
  int reversetetori; // '-o/', 0.
  int steinerleft; // '-S', 0.
  int no_sort; // 0.
  int hilbert_order; // '-b///', 52.
  int hilbert_limit; // '-b//'  8.
  int brio_threshold; // '-b' 64.
  REAL brio_ratio; // '-b/' 0.125.
  REAL facet_separate_ang_tol; // '-p', 179.9.
  REAL facet_overlap_ang_tol; // '-p/',  0.1.
  REAL facet_small_ang_tol; // '-p//', 15.0.
  REAL maxvolume; // '-a', -1.0.
  REAL minratio; // '-q', 0.0.
  REAL mindihedral; // '-q', 5.0.
  REAL optmaxdihedral; // 165.0.
  REAL optminsmtdihed; // 179.0.
  REAL optminslidihed; // 179.0.
  REAL epsilon; // '-T', 1.0e-8.
  REAL coarsen_percent; // -R1/#, 1.0.
 
  // Strings of command line arguments and input/output file names.
  char commandline[1024];
  char infilename[1024];
  char outfilename[1024];
  char addinfilename[1024];
  char bgmeshfilename[1024];
 
  // The input object of TetGen. They are recognized by either the input
  //   file extensions or by the specified options.
  // Currently the following objects are supported:
  //   - NODES, a list of nodes (.node);
  //   - POLY, a piecewise linear complex (.poly or .smesh);
  //   - OFF, a polyhedron (.off, Geomview's file format);
  //   - PLY, a polyhedron (.ply, file format from gatech, only ASCII);
  //   - STL, a surface mesh (.stl, stereolithography format);
  //   - MEDIT, a surface mesh (.mesh, Medit's file format);
  //   - MESH, a tetrahedral mesh (.ele).
  // If no extension is available, the imposed command line switch
  //   (-p or -r) implies the object.
  enum objecttype { NODES, POLY, OFF, PLY, STL, MEDIT, VTK, MESH } object;
 
  void syntax();
  void usage();
 
  // Command line parse routine.
  bool parse_commandline(int argc, char **argv);
  bool parse_commandline(char *switches)
  {
    return parse_commandline(0, &switches);
  }
 
  // Initialize all variables.
  tetgenbehavior()
  {
    plc = 0;
    psc = 0;
    refine = 0;
    quality = 0;
    nobisect = 0;
    coarsen = 0;
    metric = 0;
    weighted = 0;
    brio_hilbert = 1;
    incrflip = 0;
    flipinsert = 0;
    varvolume = 0;
    fixedvolume = 0;
    noexact = 0;
    nostaticfilter = 0;
    insertaddpoints = 0;
    regionattrib = 0;
    cdtrefine = 0;
    diagnose = 0;
    convex = 0;
    zeroindex = 0;
    facesout = 0;
    edgesout = 0;
    neighout = 0;
    voroout = 0;
    meditview = 0;
    vtkview = 0;
    nobound = 0;
    nonodewritten = 0;
    noelewritten = 0;
    nofacewritten = 0;
    noiterationnum = 0;
    nomergefacet = 0;
    nomergevertex = 0;
    nojettison = 0;
    docheck = 0;
    quiet = 0;
    verbose = 0;
 
    vertexperblock = 4092;
    tetrahedraperblock = 8188;
    shellfaceperblock = 4092;
    nobisect_nomerge = 1;
    supsteiner_level = 2;
    addsteiner_algo = 1;
    coarsen_param = 0;
    weighted_param = 0;
    fliplinklevel = -1;
    flipstarsize = -1;
    fliplinklevelinc = 1;
    reflevel = 3;
    optscheme = 7;
    optlevel = 2;
    delmaxfliplevel = 1;
    order = 1;
    reversetetori = 0;
    steinerleft = -1;
    no_sort = 0;
    hilbert_order = 52; //-1;
    hilbert_limit = 8;
    brio_threshold = 64;
    brio_ratio = 0.125;
    facet_separate_ang_tol = 179.9;
    facet_overlap_ang_tol = 0.1;
    facet_small_ang_tol = 15.0;
    maxvolume = -1.0;
    minratio = 2.0;
    mindihedral = 0.0;
    optmaxdihedral = 165.00;
    optminsmtdihed = 179.00;
    optminslidihed = 179.00;
    epsilon = 1.0e-8;
    coarsen_percent = 1.0;
    object = NODES;
 
    commandline[0] = '\0';
    infilename[0] = '\0';
    outfilename[0] = '\0';
    addinfilename[0] = '\0';
    bgmeshfilename[0] = '\0';
  }
 
}; // class tetgenbehavior
 
//                                                                           //
// tetgenmesh                                                                //
//                                                                           //
// A structure for creating and updating tetrahedral meshes.                 //
//                                                                           //
 
class tetgenmesh {
public:
  //                                                                           //
  // Mesh data structure //
  //                                                                           //
  // A tetrahedral mesh T of a 3D piecewise linear complex (PLC) X is a 3D //
  // simplicial complex whose underlying space is equal to the space of X.  T //
  // contains a 2D subcomplex S which is a triangular mesh of the boundary of //
  // X. S contains a 1D subcomplex L which is a linear mesh of the boundary of
  // // S. Faces and edges in S and L are respectively called subfaces and
  // segme- // nts to distinguish them from others in T.
  // //
  //                                                                           //
  // TetGen stores the tetrahedra and vertices of T. The basic structure of a //
  // tetrahedron contains pointers to its vertices and adjacent tetrahedra. A //
  // vertex stores its x-, y-, and z-coordinates, and a pointer to a tetrahed-
  // // ron containing it. Both tetrahedra and vertices may contain user data.
  // //
  //                                                                           //
  // Each face of T belongs to either two tetrahedra or one tetrahedron. In //
  // the latter case, the face is an exterior boundary face of T.  TetGen adds
  // // fictitious tetrahedra (one-to-one) at such faces, and connects them to
  // an // "infinite vertex" (which has no geometric coordinates).  One can
  // imagine  // such a vertex lies in 4D space and is visible by all exterior
  // boundary    // faces.  The extended set of tetrahedra (including the
  // infinite vertex) is // a tetrahedralization of a 3-pseudomanifold without
  // boundary.  It has the  // property that every face is shared by exactly two
  // tetrahedra.             //
  //                                                                           //
  // The current version of TetGen stores explicitly the subfaces and segments
  // // (which are in surface mesh S and the linear mesh L), respectively.
  // Extra // pointers are allocated in tetrahedra and subfaces to point each
  // others.   //
  //                                                                           //
 
  // The tetrahedron data structure.  It includes the following fields:
  //   - a list of four adjoining tetrahedra;
  //   - a list of four vertices;
  //   - a pointer to a list of four subfaces (optional, for -p switch);
  //   - a pointer to a list of six segments  (optional, for -p switch);
  //   - a list of user-defined floating-point attributes (optional);
  //   - a volume constraint (optional, for -a switch);
  //   - an integer of element marker (and flags);
  // The structure of a tetrahedron is an array of pointers.  Its actual size
  //   (the length of the array) is determined at runtime.
 
  typedef REAL **tetrahedron;
 
  // The subface data structure.  It includes the following fields:
  //   - a list of three adjoining subfaces;
  //   - a list of three vertices;
  //   - a list of three adjoining segments;
  //   - two adjoining tetrahedra;
  //   - an area constraint (optional, for -q switch);
  //   - an integer for boundary marker;
  //   - an integer for type, flags, etc.
 
  typedef REAL **shellface;
 
  // The point data structure.  It includes the following fields:
  //   - x, y and z coordinates;
  //   - a list of user-defined point attributes (optional);
  //   - u, v coordinates (optional, for -s switch);
  //   - a metric tensor (optional, for -q or -m switch);
  //   - a pointer to an adjacent tetrahedron;
  //   - a pointer to a parent (or a duplicate) point;
  //   - a pointer to an adjacent subface or segment (optional, -p switch);
  //   - a pointer to a tet in background mesh (optional, for -m switch);
  //   - an integer for boundary marker (point index);
  //   - an integer for point type (and flags).
  //   - an integer for geometry tag (optional, for -s switch).
  // The structure of a point is an array of REALs.  Its acutal size is
  //   determined at the runtime.
 
  typedef REAL *point;
 
  //                                                                           //
  // Handles //
  //                                                                           //
  // Navigation and manipulation in a tetrahedralization are accomplished by //
  // operating on structures referred as ``handles". A handle is a pair (t,v),
  // // where t is a pointer to a tetrahedron, and v is a 4-bit integer, in the
  // // range from 0 to 11. v is called the ``version'' of a tetrahedron, it
  // rep- // resents a directed edge of a specific face of the tetrahedron.
  // //
  //                                                                           //
  // There are 12 even permutations of the four vertices, each of them corres-
  // // ponds to a directed edge (a version) of the tetrahedron.  The 12
  // versions // can be grouped into 4 distinct ``edge rings'' in 4 ``oriented
  // faces'' of  // this tetrahedron.  One can encode each version (a directed
  // edge) into a   // 4-bit integer such that the two upper bits encode the
  // index (from 0 to 2) // of this edge in the edge ring, and the two lower
  // bits encode the index (  // from 0 to 3) of the oriented face which
  // contains this edge.               //
  //                                                                           //
  // The four vertices of a tetrahedron are indexed from 0 to 3 (according to //
  // their storage in the data structure).  Give each face the same index as //
  // the node opposite it in the tetrahedron.  Denote the edge connecting face
  // // i to face j as i/j. We number the twelve versions as follows:
  // //
  //                                                                           //
  //           |   edge 0     edge 1     edge 2 //
  //   --------|-------------------------------- //
  //    face 0 |   0 (0/1)    4 (0/3)    8 (0/2) // face 1 |   1 (1/2)    5
  //    (1/3)    9 (1/0)                               // face 2 |   2 (2/3)
  //    6 (2/1)   10 (2/0)                               // face 3 |   3 (3/0)
  //    7 (3/1)   11 (3/2)                               //
  //                                                                           //
  // Similarly, navigation and manipulation in a (boundary) triangulation are //
  // done by using handles of triangles. Each handle is a pair (s, v), where s
  // // is a pointer to a triangle, and v is a version in the range from 0 to 5.
  // // Each version corresponds to a directed edge of this triangle.
  // //
  //                                                                           //
  // Number the three vertices of a triangle from 0 to 2 (according to their //
  // storage in the data structure). Give each edge the same index as the node
  // // opposite it in the triangle. The six versions of a triangle are:
  // //
  //                                                                           //
  //                 | edge 0   edge 1   edge 2 //
  //  ---------------|-------------------------- //
  //   ccw orieation |   0        2        4 //
  //    cw orieation |   1        3        5 //
  //                                                                           //
  // In the following, a 'triface' is a handle of tetrahedron, and a 'face' is
  // // a handle of a triangle.
  // //
  //                                                                           //
 
  class triface {
  public:
    tetrahedron *tet;
    int ver; // Range from 0 to 11.
    triface() : tet(0), ver(0) {}
    triface &operator=(const triface &t)
    {
      tet = t.tet;
      ver = t.ver;
      return *this;
    }
  };
 
  class face {
  public:
    shellface *sh;
    int shver; // Range from 0 to 5.
    face() : sh(0), shver(0) {}
    face &operator=(const face &s)
    {
      sh = s.sh;
      shver = s.shver;
      return *this;
    }
  };
 
  //                                                                           //
  // Arraypool //
  //                                                                           //
  // A dynamic linear array. (It is written by J. Shewchuk) //
  //                                                                           //
  // Each arraypool contains an array of pointers to a number of blocks.  Each
  // // block contains the same fixed number of objects.  Each index of the
  // array // addresses a particular object in the pool. The most significant
  // bits add- // ress the index of the block containing the object. The less
  // significant   // bits address this object within the block.
  // //
  //                                                                           //
  // 'objectbytes' is the size of one object in blocks; 'log2objectsperblock' //
  // is the base-2 logarithm of 'objectsperblock'; 'objects' counts the number
  // // of allocated objects; 'totalmemory' is the total memory in bytes.
  // //
  //                                                                           //
 
  class arraypool {
  public:
    int objectbytes;
    int objectsperblock;
    int log2objectsperblock;
    int objectsperblockmark;
    int toparraylen;
    char **toparray;
    long objects;
    unsigned long totalmemory;
 
    void restart();
    void poolinit(int sizeofobject, int log2objperblk);
    char *getblock(int objectindex);
    void *lookup(int objectindex);
    int newindex(void **newptr);
 
    arraypool(int sizeofobject, int log2objperblk);
    ~arraypool();
  };
 
  // fastlookup() -- A fast, unsafe operation. Return the pointer to the object
  //   with a given index.  Note: The object's block must have been allocated,
  //   i.e., by the function newindex().
 
#define fastlookup(pool, index)                                                \
  (void *)((pool)->toparray[(index) >> (pool)->log2objectsperblock] +          \
           ((index) & (pool)->objectsperblockmark) * (pool)->objectbytes)
 
  //                                                                           //
  // Memorypool //
  //                                                                           //
  // A structure for memory allocation. (It is written by J. Shewchuk) //
  //                                                                           //
  // 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.
  //   //
  //                                                                           //
 
  class memorypool {
  public:
    void **firstblock, **nowblock;
    void *nextitem;
    void *deaditemstack;
    void **pathblock;
    void *pathitem;
    int alignbytes;
    int itembytes, itemwords;
    int itemsperblock;
    long items, maxitems;
    int unallocateditems;
    int pathitemsleft;
 
    memorypool();
    memorypool(int, int, int, int);
    ~memorypool();
 
    void poolinit(int, int, int, int);
    void restart();
    void *alloc();
    void dealloc(void *);
    void traversalinit();
    void *traverse();
  };
 
  //                                                                           //
  // badface //
  //                                                                           //
  // Despite of its name, a 'badface' can be used to represent one of the //
  // following objects: //
  //   - a face of a tetrahedron which is (possibly) non-Delaunay; //
  //   - an encroached subsegment or subface; //
  //   - a bad-quality tetrahedron, i.e, has too large radius-edge ratio; //
  //   - a sliver, i.e., has good radius-edge ratio but nearly zero volume; //
  //   - a recently flipped face (saved for undoing the flip later). //
  //                                                                           //
 
  class badface {
  public:
    triface tt;
    face ss;
    REAL key, cent[6]; // circumcenter or cos(dihedral angles) at 6 edges.
    point forg, fdest, fapex, foppo, noppo;
    badface *nextitem;
    badface()
      : key(0), forg(0), fdest(0), fapex(0), foppo(0), noppo(0), nextitem(0)
    {
    }
  };
 
  //                                                                           //
  // insertvertexflags //
  //                                                                           //
  // A collection of flags that pass to the routine insertvertex(). //
  //                                                                           //
 
  class insertvertexflags {
  public:
    int iloc; // input/output.
    int bowywat, lawson;
    int splitbdflag, validflag, respectbdflag;
    int rejflag, chkencflag, cdtflag;
    int assignmeshsize;
    int sloc, sbowywat;
 
    // Used by Delaunay refinement.
    int refineflag; // 0, 1, 2, 3
    triface refinetet;
    face refinesh;
    int smlenflag; // for useinsertradius.
    REAL smlen; // for useinsertradius.
    point parentpt;
 
    insertvertexflags()
    {
      iloc = bowywat = lawson = 0;
      splitbdflag = validflag = respectbdflag = 0;
      rejflag = chkencflag = cdtflag = 0;
      assignmeshsize = 0;
      sloc = sbowywat = 0;
 
      refineflag = 0;
      refinetet.tet = NULL;
      refinesh.sh = NULL;
      smlenflag = 0;
      smlen = 0.0;
    }
  };
 
  //                                                                           //
  // flipconstraints //
  //                                                                           //
  // A structure of a collection of data (options and parameters) which pass //
  // to the edge flip function flipnm(). //
  //                                                                           //
 
  class flipconstraints {
  public:
    // Elementary flip flags.
    int enqflag; // (= flipflag)
    int chkencflag;
 
    // Control flags
    int unflip; // Undo the performed flips.
    int collectnewtets; // Collect the new tets created by flips.
    int collectencsegflag;
 
    // Optimization flags.
    int remove_ndelaunay_edge; // Remove a non-Delaunay edge.
    REAL bak_tetprism_vol; // The value to be minimized.
    REAL tetprism_vol_sum;
    int remove_large_angle; // Remove a large dihedral angle at edge.
    REAL cosdihed_in; // The input cosine of the dihedral angle (> 0).
    REAL cosdihed_out; // The improved cosine of the dihedral angle.
 
    // Boundary recovery flags.
    int checkflipeligibility;
    point seg[2]; // A constraining edge to be recovered.
    point fac[3]; // A constraining face to be recovered.
    point remvert; // A vertex to be removed.
 
    flipconstraints()
    {
      enqflag = 0;
      chkencflag = 0;
 
      unflip = 0;
      collectnewtets = 0;
      collectencsegflag = 0;
 
      remove_ndelaunay_edge = 0;
      bak_tetprism_vol = 0.0;
      tetprism_vol_sum = 0.0;
      remove_large_angle = 0;
      cosdihed_in = 0.0;
      cosdihed_out = 0.0;
 
      checkflipeligibility = 0;
      seg[0] = NULL;
      fac[0] = NULL;
      remvert = NULL;
    }
  };
 
  //                                                                           //
  // optparameters //
  //                                                                           //
  // Optimization options and parameters. //
  //                                                                           //
 
  class optparameters {
  public:
    // The one of goals of optimization.
    int max_min_volume; // Maximize the minimum volume.
    int min_max_aspectratio; // Minimize the maximum aspect ratio.
    int min_max_dihedangle; // Minimize the maximum dihedral angle.
 
    // The initial and improved value.
    REAL initval, imprval;
 
    int numofsearchdirs;
    REAL searchstep;
    int maxiter; // Maximum smoothing iterations (disabled by -1).
    int smthiter; // Performed iterations.
 
    optparameters()
    {
      max_min_volume = 0;
      min_max_aspectratio = 0;
      min_max_dihedangle = 0;
 
      initval = imprval = 0.0;
 
      numofsearchdirs = 10;
      searchstep = 0.01;
      maxiter = -1; // Unlimited smoothing iterations.
      smthiter = 0;
    }
  };
 
  //                                                                           //
  // Labels (enumeration declarations) used by TetGen. //
  //                                                                           //
 
  // Labels that signify the type of a vertex.
  enum verttype {
    UNUSEDVERTEX,
    DUPLICATEDVERTEX,
    RIDGEVERTEX,
    ACUTEVERTEX,
    FACETVERTEX,
    VOLVERTEX,
    FREESEGVERTEX,
    FREEFACETVERTEX,
    FREEVOLVERTEX,
    NREGULARVERTEX,
    DEADVERTEX
  };
 
  // Labels that signify the result of triangle-triangle intersection test.
  enum interresult {
    DISJOINT,
    INTERSECT,
    SHAREVERT,
    SHAREEDGE,
    SHAREFACE,
    TOUCHEDGE,
    TOUCHFACE,
    ACROSSVERT,
    ACROSSEDGE,
    ACROSSFACE
  };
 
  // Labels that signify the result of point location.
  enum locateresult {
    UNKNOWN,
    OUTSIDE,
    INTETRAHEDRON,
    ONFACE,
    ONEDGE,
    ONVERTEX,
    ENCVERTEX,
    ENCSEGMENT,
    ENCSUBFACE,
    NEARVERTEX,
    NONREGULAR,
    INSTAR,
    BADELEMENT
  };
 
  //                                                                           //
  // Variables of TetGen //
  //                                                                           //
 
  // Pointer to the input data (a set of nodes, a PLC, or a mesh).
  tetgenio *in, *addin;
 
  // Pointer to the switches and parameters.
  tetgenbehavior *b;
 
  // Pointer to a background mesh (contains size specification map).
  tetgenmesh *bgm;
 
  // Memorypools to store mesh elements (points, tetrahedra, subfaces, and
  //   segments) and extra pointers between tetrahedra, subfaces, and segments.
  memorypool *tetrahedrons, *subfaces, *subsegs, *points;
  memorypool *tet2subpool, *tet2segpool;
 
  // Memorypools to store bad-quality (or encroached) elements.
  memorypool *badtetrahedrons, *badsubfacs, *badsubsegs;
 
  // A memorypool to store faces to be flipped.
  memorypool *flippool;
  arraypool *unflipqueue;
  badface *flipstack;
 
  // Arrays used for point insertion (the Bowyer-Watson algorithm).
  arraypool *cavetetlist, *cavebdrylist, *caveoldtetlist;
  arraypool *cavetetshlist, *cavetetseglist, *cavetetvertlist;
  arraypool *caveencshlist, *caveencseglist;
  arraypool *caveshlist, *caveshbdlist, *cavesegshlist;
 
  // Stacks used for CDT construction and boundary recovery.
  arraypool *subsegstack, *subfacstack, *subvertstack;
 
  // Arrays of encroached segments and subfaces (for mesh refinement).
  arraypool *encseglist, *encshlist;
 
  // The map between facets to their vertices (for mesh refinement).
  int *idx2facetlist;
  point *facetverticeslist;
 
  // The map between segments to their endpoints (for mesh refinement).
  point *segmentendpointslist;
 
  // The infinite vertex.
  point dummypoint;
  // The recently visited tetrahedron, subface.
  triface recenttet;
  face recentsh;
 
  // PI is the ratio of a circle's circumference to its diameter.
  static REAL PI;
 
  // Array (size = numberoftetrahedra * 6) for storing high-order nodes of
  //   tetrahedra (only used when -o2 switch is selected).
  point *highordertable;
 
  // Various variables.
  int numpointattrib; // Number of point attributes.
  int numelemattrib; // Number of tetrahedron attributes.
  int sizeoftensor; // Number of REALs per metric tensor.
  int pointmtrindex; // Index to find the metric tensor of a point.
  int pointparamindex; // Index to find the u,v coordinates of a point.
  int point2simindex; // Index to find a simplex adjacent to a point.
  int pointmarkindex; // Index to find boundary marker of a point.
  int pointinsradiusindex; // Index to find the insertion radius of a point.
  int elemattribindex; // Index to find attributes of a tetrahedron.
  int volumeboundindex; // Index to find volume bound of a tetrahedron.
  int elemmarkerindex; // Index to find marker of a tetrahedron.
  int shmarkindex; // Index to find boundary marker of a subface.
  int areaboundindex; // Index to find area bound of a subface.
  int checksubsegflag; // Are there segments in the tetrahedralization yet?
  int checksubfaceflag; // Are there subfaces in the tetrahedralization yet?
  int checkconstraints; // Are there variant (node, seg, facet) constraints?
  int nonconvex; // Is current mesh non-convex?
  int autofliplinklevel; // The increase of link levels, default is 1.
  int useinsertradius; // Save the insertion radius for Steiner points.
  long samples; // Number of random samples for point location.
  unsigned long randomseed; // Current random number seed.
  REAL cosmaxdihed, cosmindihed; // The cosine values of max/min dihedral.
  REAL cossmtdihed; // The cosine value of a bad dihedral to be smoothed.
  REAL cosslidihed; // The cosine value of the max dihedral of a sliver.
  REAL minfaceang, minfacetdihed; // The minimum input (dihedral) angles.
  REAL tetprism_vol_sum; // The total volume of tetrahedral-prisms (in 4D).
  REAL longest; // The longest possible edge length.
  REAL minedgelength; // = longest * b->epsion.
  REAL xmax, xmin, ymax, ymin, zmax, zmin; // Bounding box of points.
 
  // Counters.
  long insegments; // Number of input segments.
  long hullsize; // Number of exterior boundary faces.
  long meshedges; // Number of mesh edges.
  long meshhulledges; // Number of boundary mesh edges.
  long steinerleft; // Number of Steiner points not yet used.
  long dupverts; // Are there duplicated vertices?
  long unuverts; // Are there unused vertices?
  long nonregularcount; // Are there non-regular vertices?
  long st_segref_count, st_facref_count, st_volref_count; // Steiner points.
  long fillregioncount, cavitycount, cavityexpcount;
  long flip14count, flip26count, flipn2ncount;
  long flip23count, flip32count, flip44count, flip41count;
  long flip31count, flip22count;
  unsigned long totalworkmemory; // Total memory used by working arrays.
 
  //                                                                           //
  // Mesh manipulation primitives //
  //                                                                           //
 
  // Fast lookup tables for mesh manipulation primitives.
  static int bondtbl[12][12], fsymtbl[12][12];
  static int esymtbl[12], enexttbl[12], eprevtbl[12];
  static int enextesymtbl[12], eprevesymtbl[12];
  static int eorgoppotbl[12], edestoppotbl[12];
  static int facepivot1[12], facepivot2[12][12];
  static int orgpivot[12], destpivot[12], apexpivot[12], oppopivot[12];
  static int tsbondtbl[12][6], stbondtbl[12][6];
  static int tspivottbl[12][6], stpivottbl[12][6];
  static int ver2edge[12], edge2ver[6], epivot[12];
  static int sorgpivot[6], sdestpivot[6], sapexpivot[6];
  static int snextpivot[6];
 
  void inittables();
 
  // Primitives for tetrahedra.
  inline tetrahedron encode(triface &t);
  inline tetrahedron encode2(tetrahedron *ptr, int ver);
  inline void decode(tetrahedron ptr, triface &t);
  inline void bond(triface &t1, triface &t2);
  inline void dissolve(triface &t);
  inline void esym(triface &t1, triface &t2);
  inline void esymself(triface &t);
  inline void enext(triface &t1, triface &t2);
  inline void enextself(triface &t);
  inline void eprev(triface &t1, triface &t2);
  inline void eprevself(triface &t);
  inline void enextesym(triface &t1, triface &t2);
  inline void enextesymself(triface &t);
  inline void eprevesym(triface &t1, triface &t2);
  inline void eprevesymself(triface &t);
  inline void eorgoppo(triface &t1, triface &t2);
  inline void eorgoppoself(triface &t);
  inline void edestoppo(triface &t1, triface &t2);
  inline void edestoppoself(triface &t);
  inline void fsym(triface &t1, triface &t2);
  inline void fsymself(triface &t);
  inline void fnext(triface &t1, triface &t2);
  inline void fnextself(triface &t);
  inline point org(triface &t);
  inline point dest(triface &t);
  inline point apex(triface &t);
  inline point oppo(triface &t);
  inline void setorg(triface &t, point p);
  inline void setdest(triface &t, point p);
  inline void setapex(triface &t, point p);
  inline void setoppo(triface &t, point p);
  inline REAL elemattribute(tetrahedron *ptr, int attnum);
  inline void setelemattribute(tetrahedron *ptr, int attnum, REAL value);
  inline REAL volumebound(tetrahedron *ptr);
  inline void setvolumebound(tetrahedron *ptr, REAL value);
  inline int elemindex(tetrahedron *ptr);
  inline void setelemindex(tetrahedron *ptr, int value);
  inline int elemmarker(tetrahedron *ptr);
  inline void setelemmarker(tetrahedron *ptr, int value);
  inline void infect(triface &t);
  inline void uninfect(triface &t);
  inline bool infected(triface &t);
  inline void marktest(triface &t);
  inline void unmarktest(triface &t);
  inline bool marktested(triface &t);
  inline void markface(triface &t);
  inline void unmarkface(triface &t);
  inline bool facemarked(triface &t);
  inline void markedge(triface &t);
  inline void unmarkedge(triface &t);
  inline bool edgemarked(triface &t);
  inline void marktest2(triface &t);
  inline void unmarktest2(triface &t);
  inline bool marktest2ed(triface &t);
  inline int elemcounter(triface &t);
  inline void setelemcounter(triface &t, int value);
  inline void increaseelemcounter(triface &t);
  inline void decreaseelemcounter(triface &t);
  inline bool ishulltet(triface &t);
  inline bool isdeadtet(triface &t);
 
  // Primitives for subfaces and subsegments.
  inline void sdecode(shellface sptr, face &s);
  inline shellface sencode(face &s);
  inline shellface sencode2(shellface *sh, int shver);
  inline void spivot(face &s1, face &s2);
  inline void spivotself(face &s);
  inline void sbond(face &s1, face &s2);
  inline void sbond1(face &s1, face &s2);
  inline void sdissolve(face &s);
  inline point sorg(face &s);
  inline point sdest(face &s);
  inline point sapex(face &s);
  inline void setsorg(face &s, point pointptr);
  inline void setsdest(face &s, point pointptr);
  inline void setsapex(face &s, point pointptr);
  inline void sesym(face &s1, face &s2);
  inline void sesymself(face &s);
  inline void senext(face &s1, face &s2);
  inline void senextself(face &s);
  inline void senext2(face &s1, face &s2);
  inline void senext2self(face &s);
  inline REAL areabound(face &s);
  inline void setareabound(face &s, REAL value);
  inline int shellmark(face &s);
  inline void setshellmark(face &s, int value);
  inline void sinfect(face &s);
  inline void suninfect(face &s);
  inline bool sinfected(face &s);
  inline void smarktest(face &s);
  inline void sunmarktest(face &s);
  inline bool smarktested(face &s);
  inline void smarktest2(face &s);
  inline void sunmarktest2(face &s);
  inline bool smarktest2ed(face &s);
  inline void smarktest3(face &s);
  inline void sunmarktest3(face &s);
  inline bool smarktest3ed(face &s);
  inline void setfacetindex(face &f, int value);
  inline int getfacetindex(face &f);
 
  // Primitives for interacting tetrahedra and subfaces.
  inline void tsbond(triface &t, face &s);
  inline void tsdissolve(triface &t);
  inline void stdissolve(face &s);
  inline void tspivot(triface &t, face &s);
  inline void stpivot(face &s, triface &t);
 
  // Primitives for interacting tetrahedra and segments.
  inline void tssbond1(triface &t, face &seg);
  inline void sstbond1(face &s, triface &t);
  inline void tssdissolve1(triface &t);
  inline void sstdissolve1(face &s);
  inline void tsspivot1(triface &t, face &s);
  inline void sstpivot1(face &s, triface &t);
 
  // Primitives for interacting subfaces and segments.
  inline void ssbond(face &s, face &edge);
  inline void ssbond1(face &s, face &edge);
  inline void ssdissolve(face &s);
  inline void sspivot(face &s, face &edge);
 
  // Primitives for points.
  inline int pointmark(point pt);
  inline void setpointmark(point pt, int value);
  inline enum verttype pointtype(point pt);
  inline void setpointtype(point pt, enum verttype value);
  inline int pointgeomtag(point pt);
  inline void setpointgeomtag(point pt, int value);
  inline REAL pointgeomuv(point pt, int i);
  inline void setpointgeomuv(point pt, int i, REAL value);
  inline void pinfect(point pt);
  inline void puninfect(point pt);
  inline bool pinfected(point pt);
  inline void pmarktest(point pt);
  inline void punmarktest(point pt);
  inline bool pmarktested(point pt);
  inline void pmarktest2(point pt);
  inline void punmarktest2(point pt);
  inline bool pmarktest2ed(point pt);
  inline void pmarktest3(point pt);
  inline void punmarktest3(point pt);
  inline bool pmarktest3ed(point pt);
  inline tetrahedron point2tet(point pt);
  inline void setpoint2tet(point pt, tetrahedron value);
  inline shellface point2sh(point pt);
  inline void setpoint2sh(point pt, shellface value);
  inline point point2ppt(point pt);
  inline void setpoint2ppt(point pt, point value);
  inline tetrahedron point2bgmtet(point pt);
  inline void setpoint2bgmtet(point pt, tetrahedron value);
  inline void setpointinsradius(point pt, REAL value);
  inline REAL getpointinsradius(point pt);
  inline bool issteinerpoint(point pt);
 
  // Advanced primitives.
  inline void point2tetorg(point pt, triface &t);
  inline void point2shorg(point pa, face &s);
  inline point farsorg(face &seg);
  inline point farsdest(face &seg);
 
  //                                                                           //
  //  Memory managment //
  //                                                                           //
 
  void tetrahedrondealloc(tetrahedron *);
  tetrahedron *tetrahedrontraverse();
  tetrahedron *alltetrahedrontraverse();
  void shellfacedealloc(memorypool *, shellface *);
  shellface *shellfacetraverse(memorypool *);
  void pointdealloc(point);
  point pointtraverse();
 
  void makeindex2pointmap(point *&);
  void makepoint2submap(memorypool *, int *&, face *&);
  void maketetrahedron(triface *);
  void makeshellface(memorypool *, face *);
  void makepoint(point *, enum verttype);
 
  void initializepools();
 
  //                                                                           //
  // Advanced geometric predicates and calculations //
  //                                                                           //
  // TetGen uses a simplified symbolic perturbation scheme from Edelsbrunner, //
  // et al [*].  Hence the point-in-sphere test never returns a zero. The idea
  // // is to perturb the weights of vertices in the fourth dimension.  TetGen
  // // uses the indices of the vertices decide the amount of perturbation. It
  // is // implemented in the routine insphere_s().
  //                                                                           //
  // The routine tri_edge_test() determines whether or not a triangle and an //
  // edge intersect in 3D. If they intersect, their intersection type is also //
  // reported. This test is a combination of n 3D orientation tests (n is bet-
  // // ween 3 and 9). It uses the robust orient3d() test to make the branch
  // dec- // isions.  The routine tri_tri_test() determines whether or not two
  // triang- // les intersect in 3D. It also uses the robust orient3d() test.
  // //
  //                                                                           //
  // There are a number of routines to calculate geometrical quantities, e.g.,
  // // circumcenters, angles, dihedral angles, face normals, face areas, etc.
  // // They are so far done by the default floating-point arithmetics which are
  // // non-robust. They should be improved in the future.
  // //
  //                                                                           //
 
  // Symbolic perturbations (robust)
  REAL insphere_s(REAL *, REAL *, REAL *, REAL *, REAL *);
  REAL orient4d_s(REAL *, REAL *, REAL *, REAL *, REAL *, REAL, REAL, REAL,
                  REAL, REAL);
 
  // Triangle-edge intersection test (robust)
  int tri_edge_2d(point, point, point, point, point, point, int, int *, int *);
  int tri_edge_tail(point, point, point, point, point, point, REAL, REAL, int,
                    int *, int *);
  int tri_edge_test(point, point, point, point, point, point, int, int *,
                    int *);
 
  // Triangle-triangle intersection test (robust)
  int tri_edge_inter_tail(point, point, point, point, point, REAL, REAL);
  int tri_tri_inter(point, point, point, point, point, point);
 
  // Linear algebra functions
  inline REAL dot(REAL *v1, REAL *v2);
  inline void cross(REAL *v1, REAL *v2, REAL *n);
  bool lu_decmp(REAL lu[4][4], int n, int *ps, REAL *d, int N);
  void lu_solve(REAL lu[4][4], int n, int *ps, REAL *b, int N);
 
  // An embedded 2-dimensional geometric predicate (non-robust)
  REAL incircle3d(point pa, point pb, point pc, point pd);
 
  // Geometric calculations (non-robust)
  REAL orient3dfast(REAL *pa, REAL *pb, REAL *pc, REAL *pd);
  inline REAL norm2(REAL x, REAL y, REAL z);
  inline REAL distance(REAL *p1, REAL *p2);
  void facenormal(point pa, point pb, point pc, REAL *n, int pivot, REAL *lav);
  REAL shortdistance(REAL *p, REAL *e1, REAL *e2);
  REAL triarea(REAL *pa, REAL *pb, REAL *pc);
  REAL interiorangle(REAL *o, REAL *p1, REAL *p2, REAL *n);
  void projpt2edge(REAL *p, REAL *e1, REAL *e2, REAL *prj);
  void projpt2face(REAL *p, REAL *f1, REAL *f2, REAL *f3, REAL *prj);
  bool tetalldihedral(point, point, point, point, REAL *, REAL *, REAL *);
  void tetallnormal(point, point, point, point, REAL N[4][3], REAL *volume);
  REAL tetaspectratio(point, point, point, point);
  bool circumsphere(REAL *, REAL *, REAL *, REAL *, REAL *cent, REAL *radius);
  bool orthosphere(REAL *, REAL *, REAL *, REAL *, REAL, REAL, REAL, REAL,
                   REAL *, REAL *);
  void planelineint(REAL *, REAL *, REAL *, REAL *, REAL *, REAL *, REAL *);
  int linelineint(REAL *, REAL *, REAL *, REAL *, REAL *, REAL *, REAL *,
                  REAL *);
  REAL tetprismvol(REAL *pa, REAL *pb, REAL *pc, REAL *pd);
  bool calculateabovepoint(arraypool *, point *, point *, point *);
  void calculateabovepoint4(point, point, point, point);
 
  // PLC error reports.
  void report_overlapping_facets(face *, face *, REAL dihedang = 0.0);
  int report_selfint_edge(point, point, face *sedge, triface *searchtet,
                          enum interresult);
  int report_selfint_face(point, point, point, face *sface, triface *iedge,
                          int intflag, int *types, int *poss);
 
  //                                                                           //
  // Local mesh transformations //
  //                                                                           //
  // A local transformation replaces a small set of tetrahedra with another //
  // set of tetrahedra which fills the same space and the same boundaries. //
  //   In 3D, the most simplest local transformations are the elementary flips
  //   //
  // performed within the convex hull of five vertices: 2-to-3, 3-to-2,
  // 1-to-4,// and 4-to-1 flips,  where the numbers indicate the number of
  // tetrahedra    // before and after each flip.  The 1-to-4 and 4-to-1 flip
  // involve inserting // or deleting a vertex, respectively.
  // //
  //   There are complex local transformations which can be decomposed as a //
  // combination of elementary flips. For example,a 4-to-4 flip which replaces
  // // two coplanar edges can be regarded by a 2-to-3 flip and a 3-to-2 flip.
  // // Note that the first 2-to-3 flip will temporarily create a degenerate
  // tet- // rahedron which is removed immediately by the followed 3-to-2 flip.
  // More  // generally, a n-to-m flip, where n > 3, m = (n - 2) * 2, which
  // removes an  // edge can be done by first performing a sequence of (n - 3)
  // 2-to-3 flips   // followed by a 3-to-2 flip.
  // //
  //                                                                           //
  // The routines flip23(), flip32(), and flip41() perform the three element- //
  // ray flips. The flip14() is available inside the routine insertpoint(). //
  //                                                                           //
  // The routines flipnm() and flipnm_post() implement a generalized edge flip
  // // algorithm which uses a combination of elementary flips.
  // //
  //                                                                           //
  // The routine insertpoint() implements a variant of Bowyer-Watson's cavity //
  // algorithm to insert a vertex. It works for arbitrary tetrahedralization, //
  // either Delaunay, or constrained Delaunay, or non-Delaunay. //
  //                                                                           //
 
  // The elementary flips.
  void flip23(triface *, int, flipconstraints *fc);
  void flip32(triface *, int, flipconstraints *fc);
  void flip41(triface *, int, flipconstraints *fc);
 
  // A generalized edge flip.
  int flipnm(triface *, int n, int level, int, flipconstraints *fc);
  int flipnm_post(triface *, int n, int nn, int, flipconstraints *fc);
 
  // Point insertion.
  int insertpoint(point, triface *, face *, face *, insertvertexflags *);
  void insertpoint_abort(face *, insertvertexflags *);
 
  //                                                                           //
  // Delaunay tetrahedralization //
  //                                                                           //
  // The routine incrementaldelaunay() implemented two incremental algorithms //
  // for constructing Delaunay tetrahedralizations (DTs):  the Bowyer-Watson //
  // (B-W) algorithm and the incremental flip algorithm of Edelsbrunner and //
  // Shah, "Incremental topological flipping works for regular triangulation,"
  // // Algorithmica, 15:233-241, 1996.
  // //
  //                                                                           //
  // The routine incrementalflip() implements the flip algorithm of [Edelsbru-
  // // nner and Shah, 1996].  It flips a queue of locally non-Delaunay faces
  // (in // an arbitrary order).  The success is guaranteed when the Delaunay
  // tetrah- // edralization is constructed incrementally by adding one vertex
  // at a time. //
  //                                                                           //
  // The routine locate() finds a tetrahedron contains a new point in current //
  // DT.  It uses a simple stochastic walk algorithm: starting from an arbitr-
  // // ary tetrahedron in DT, it finds the destination by visit one tetrahedron
  // // at a time, randomly chooses a tetrahedron if there are more than one
  // // choices. This algorithm terminates due to Edelsbrunner's acyclic
  // theorem. //
  //   Choose a good starting tetrahedron is crucial to the speed of the walk.
  //   //
  // TetGen originally uses the "jump-and-walk" algorithm of Muecke, E.P., //
  // Saias, I., and Zhu, B. "Fast Randomized Point Location Without Preproces-
  // // sing." In Proceedings of the 12th ACM Symposium on Computational
  // Geometry,// 274-283, 1996.  It first randomly samples several tetrahedra in
  // the DT    // and then choosing the closet one to start walking.
  // //
  //   The above algorithm slows download dramatically as the number of points
  //   //
  // grows -- reported in Amenta, N., Choi, S. and Rote, G., "Incremental //
  // construction con {BRIO}," In Proceedings of 19th ACM Symposium on //
  // Computational Geometry, 211-219, 2003.  On the other hand, Liu and //
  // Snoeyink showed that the point location can be made in constant time if //
  // the points are pre-sorted so that the nearby points in space have nearby //
  // indices, then adding the points in this order. They sorted the points //
  // along the 3D Hilbert curve. //
  //                                                                           //
  // The routine hilbert_sort3() sorts a set of 3D points along the 3D Hilbert
  // // curve. It recursively splits a point set according to the Hilbert
  // indices // mapped to the subboxes of the bounding box of the point set.
  // //
  //   The Hilbert indices is calculated by Butz's algorithm in 1971.  A nice //
  // exposition of this algorithm can be found in the paper of Hamilton, C., //
  // "Compact Hilbert Indices", Technical Report CS-2006-07, Computer Science,
  // // Dalhousie University, 2006 (the Section 2). My implementation also
  // refer- // enced Steven Witham's implementation of "Hilbert walk"
  // (hopefully, it is  // still available at:
  // http://www.tiac.net/~sw/2008/10/Hilbert/).            //
  //                                                                           //
  // TetGen sorts the points using the method in the paper of
  // Boissonnat,J.-D.,// Devillers, O. and Hornus, S. "Incremental Construction
  // of the Delaunay    // Triangulation and the Delaunay Graph in Medium
  // Dimension," In Proceedings // of the 25th ACM Symposium on Computational
  // Geometry, 2009.                //
  //   It first randomly sorts the points into subgroups using the Biased
  //   Rand-//
  // omized Insertion Ordering (BRIO) of Amenta et al 2003, then sorts the //
  // points in each subgroup along the 3D Hilbert curve.  Inserting points in //
  // this order ensures a randomized "sprinkling" of the points over the //
  // domain, while sorting of each subset ensures locality. //
  //                                                                           //
 
  void transfernodes();
 
  // Point sorting.
  int transgc[8][3][8], tsb1mod3[8];
  void hilbert_init(int n);
  int hilbert_split(point *vertexarray, int arraysize, int gc0, int gc1, REAL,
                    REAL, REAL, REAL, REAL, REAL);
  void hilbert_sort3(point *vertexarray, int arraysize, int e, int d, REAL,
                     REAL, REAL, REAL, REAL, REAL, int depth);
  void brio_multiscale_sort(point *, int, int threshold, REAL ratio,
                            int *depth);
 
  // Point location.
  unsigned long randomnation(unsigned int choices);
  void randomsample(point searchpt, triface *searchtet);
  enum locateresult locate(point searchpt, triface *searchtet,
                           int chkencflag = 0);
 
  // Incremental flips.
  void flippush(badface *&, triface *);
  int incrementalflip(point newpt, int, flipconstraints *fc);
 
  // Incremental Delaunay construction.
  void initialdelaunay(point pa, point pb, point pc, point pd);
  void incrementaldelaunay(clock_t &);
 
  //                                                                           //
  // Surface triangulation //
  //                                                                           //
 
  void flipshpush(face *);
  void flip22(face *, int, int);
  void flip31(face *, int);
  long lawsonflip();
  int sinsertvertex(point newpt, face *, face *, int iloc, int bowywat, int);
  int sremovevertex(point delpt, face *, face *, int lawson);
 
  enum locateresult slocate(point, face *, int, int, int);
  enum interresult sscoutsegment(face *, point, int, int, int);
  void scarveholes(int, REAL *);
 
  void unifysegments();
  void identifyinputedges(point *);
  void mergefacets();
 
  enum interresult finddirection(triface *searchtet, point endpt);
 
  //                                                                           //
  // Constrained tetrahedralizations. //
  //                                                                           //
 
  int checkflipeligibility(int fliptype, point, point, point, point, point,
                           int level, int edgepivot, flipconstraints *fc);
 
  int removeedgebyflips(triface *, flipconstraints *);
  int removefacebyflips(triface *, flipconstraints *);
 
  int recoveredgebyflips(point, point, face *, triface *, int fullsearch);
  int add_steinerpt_in_schoenhardtpoly(triface *, int, int chkencflag);
  int add_steinerpt_in_segment(face *, int searchlevel);
  int addsteiner4recoversegment(face *, int);
  int recoversegments(arraypool *, int fullsearch, int steinerflag);
 
  int recoverfacebyflips(point, point, point, face *, triface *);
  int recoversubfaces(arraypool *, int steinerflag);
 
  int getvertexstar(int, point searchpt, arraypool *, arraypool *, arraypool *);
  int getedge(point, point, triface *);
  int reduceedgesatvertex(point startpt, arraypool *endptlist);
  int removevertexbyflips(point steinerpt);
 
  int suppressbdrysteinerpoint(point steinerpt);
  int suppresssteinerpoints();
 
  void recoverboundary(clock_t &);
 
  //                                                                           //
  // Mesh reconstruction //
  //                                                                           //
 
  void carveholes();
 
  // Comment: These three functions are implemented directly in:
  // src/mesh/meshGRegionBoundaryRecovery.cpp
  bool reconstructmesh(void *);
  void outsurfacemesh(const char *mfilename);
  void outmesh2medit(const char *mfilename);
 
  void enqueuesubface(memorypool *, face *);
  void enqueuetetrahedron(triface *);
 
  //                                                                           //
  // Mesh optimization //
  //                                                                           //
 
  long lawsonflip3d(flipconstraints *fc);
  void recoverdelaunay();
 
  int gettetrahedron(point, point, point, point, triface *);
  long improvequalitybyflips();
 
  int smoothpoint(point smtpt, arraypool *, int ccw, optparameters *opm);
  long improvequalitybysmoothing(optparameters *opm);
 
  int splitsliver(triface *, REAL, int);
  long removeslivers(int);
 
  void optimizemesh();
 
  void jettisonnodes();
 
  //                                                                           //
  // Constructor & destructor //
  //                                                                           //
 
  void initializetetgenmesh()
  {
    in = addin = NULL;
    b = NULL;
    bgm = NULL;
 
    tetrahedrons = subfaces = subsegs = points = NULL;
    badtetrahedrons = badsubfacs = badsubsegs = NULL;
    tet2segpool = tet2subpool = NULL;
    flippool = NULL;
 
    dummypoint = NULL;
    flipstack = NULL;
    unflipqueue = NULL;
 
    cavetetlist = cavebdrylist = caveoldtetlist = NULL;
    cavetetshlist = cavetetseglist = cavetetvertlist = NULL;
    caveencshlist = caveencseglist = NULL;
    caveshlist = caveshbdlist = cavesegshlist = NULL;
 
    subsegstack = subfacstack = subvertstack = NULL;
    encseglist = encshlist = NULL;
    idx2facetlist = NULL;
    facetverticeslist = NULL;
    segmentendpointslist = NULL;
 
    highordertable = NULL;
 
    numpointattrib = numelemattrib = 0;
    sizeoftensor = 0;
    pointmtrindex = 0;
    pointparamindex = 0;
    pointmarkindex = 0;
    point2simindex = 0;
    pointinsradiusindex = 0;
    elemattribindex = 0;
    volumeboundindex = 0;
    shmarkindex = 0;
    areaboundindex = 0;
    checksubsegflag = 0;
    checksubfaceflag = 0;
    checkconstraints = 0;
    nonconvex = 0;
    autofliplinklevel = 1;
    useinsertradius = 0;
    samples = 0l;
    randomseed = 1l;
    minfaceang = minfacetdihed = PI;
    tetprism_vol_sum = 0.0;
    longest = minedgelength = 0.0;
    xmax = xmin = ymax = ymin = zmax = zmin = 0.0;
 
    insegments = 0l;
    hullsize = 0l;
    meshedges = meshhulledges = 0l;
    steinerleft = -1;
    dupverts = 0l;
    unuverts = 0l;
    nonregularcount = 0l;
    st_segref_count = st_facref_count = st_volref_count = 0l;
    fillregioncount = cavitycount = cavityexpcount = 0l;
    flip14count = flip26count = flipn2ncount = 0l;
    flip23count = flip32count = flip44count = flip41count = 0l;
    flip22count = flip31count = 0l;
    totalworkmemory = 0l;
 
  } // tetgenmesh()
 
  void freememory()
  {
    if(bgm != NULL) { delete bgm; }
 
    if(points != (memorypool *)NULL) {
      delete points;
      delete[] dummypoint;
    }
    if(tetrahedrons != (memorypool *)NULL) { delete tetrahedrons; }
    if(subfaces != (memorypool *)NULL) {
      delete subfaces;
      delete subsegs;
    }
    if(tet2segpool != NULL) {
      delete tet2segpool;
      delete tet2subpool;
    }
 
    if(badtetrahedrons) { delete badtetrahedrons; }
    if(badsubfacs) { delete badsubfacs; }
    if(badsubsegs) { delete badsubsegs; }
    if(encseglist) { delete encseglist; }
    if(encshlist) { delete encshlist; }
 
    if(flippool != NULL) {
      delete flippool;
      delete unflipqueue;
    }
 
    if(cavetetlist != NULL) {
      delete cavetetlist;
      delete cavebdrylist;
      delete caveoldtetlist;
      delete cavetetvertlist;
    }
 
    if(caveshlist != NULL) {
      delete caveshlist;
      delete caveshbdlist;
      delete cavesegshlist;
      delete cavetetshlist;
      delete cavetetseglist;
      delete caveencshlist;
      delete caveencseglist;
    }
 
    if(subsegstack != NULL) {
      delete subsegstack;
      delete subfacstack;
      delete subvertstack;
    }
 
    if(idx2facetlist != NULL) {
      delete[] idx2facetlist;
      delete[] facetverticeslist;
    }
 
    if(segmentendpointslist != NULL) { delete[] segmentendpointslist; }
 
    if(highordertable != NULL) { delete[] highordertable; }
 
    initializetetgenmesh();
  }
 
  tetgenmesh() { initializetetgenmesh(); }
 
  ~tetgenmesh() { freememory(); } // ~tetgenmesh()
 
}; // End of class tetgenmesh.
 
//                                                                           //
// terminatetetgen()    Terminate TetGen with a given exit code.             //
//                                                                           //
 
// selfint_event, a structure to report self-intersections.
//
//   - e_type,  report the type of self-intersections,
//     it may be one of:
//     0, reserved.
//     1, two edges intersect,
//     2, an edge and a triangle intersect,
//     3, two triangles intersect,
//     4, two edges are overlapping,
//     5, an edge and a triangle are overlapping,
//     6, two triangles are overlapping,
//     7, a vertex lies in an edge,
//     8, a vertex lies in a facet,
 
class selfint_event {
public:
  int e_type;
  int f_marker1; // Tag of the 1st facet.
  int s_marker1; // Tag of the 1st segment.
  int f_vertices1[3];
  int f_marker2; // Tag of the 2nd facet.
  int s_marker2; // Tag of the 2nd segment.
  int f_vertices2[3];
  REAL int_point[3];
  selfint_event()
  {
    e_type = 0;
    f_marker1 = f_marker2 = 0;
    s_marker1 = s_marker2 = 0;
  }
};
 
static selfint_event sevent;
 
inline void terminatetetgen(tetgenmesh *m, int x)
{
#ifdef TETLIBRARY
  throw x;
#else
  switch(x) {
  case 1: // Out of memory.
    printf("Error:  Out of memory.\n");
    break;
  case 2: // Encounter an internal error.
    printf("Please report this bug to Hang.Si@wias-berlin.de. Include\n");
    printf("  the message above, your input data set, and the exact\n");
    printf("  command line you used to run this program, thank you.\n");
    break;
  case 3:
    printf("A self-intersection was detected. Program stopped.\n");
    printf("Hint: use -d option to detect all self-intersections.\n");
    break;
  case 4:
    printf("A very small input feature size was detected. Program stopped.\n");
    if(m) {
      printf("Hint: use -T option to set a smaller tolerance. Current is %g\n",
             m->b->epsilon);
    }
    break;
  case 5:
    printf("Two very close input facets were detected. Program stopped.\n");
    printf("Hint: use -Y option to avoid adding Steiner points in boundary.\n");
    break;
  case 10: printf("An input error was detected. Program stopped.\n"); break;
  } // switch (x)
  exit(x);
#endif // #ifdef TETLIBRARY
}
 
//                                                                           //
// Primitives for tetrahedra                                                 //
//                                                                           //
 
// encode()  compress a handle into a single pointer.  It relies on the
//   assumption that all addresses of tetrahedra are aligned to sixteen-
//   byte boundaries, so that the last four significant bits are zero.
 
inline tetgenmesh::tetrahedron tetgenmesh::encode(triface &t)
{
  return (tetrahedron)((uintptr_t)(t).tet | (uintptr_t)(t).ver);
}
 
inline tetgenmesh::tetrahedron tetgenmesh::encode2(tetrahedron *ptr, int ver)
{
  return (tetrahedron)((uintptr_t)(ptr) | (uintptr_t)(ver));
}
 
// decode()  converts a pointer to a handle. The version is extracted from
//   the four least significant bits of the pointer.
 
inline void tetgenmesh::decode(tetrahedron ptr, triface &t)
{
  (t).ver = (int)((uintptr_t)(ptr) & (uintptr_t)15);
  (t).tet = (tetrahedron *)((uintptr_t)(ptr) ^ (uintptr_t)(t).ver);
}
 
// bond()  connects two tetrahedra together. (t1,v1) and (t2,v2) must
//   refer to the same face and the same edge.
 
inline void tetgenmesh::bond(triface &t1, triface &t2)
{
  t1.tet[t1.ver & 3] = encode2(t2.tet, bondtbl[t1.ver][t2.ver]);
  t2.tet[t2.ver & 3] = encode2(t1.tet, bondtbl[t2.ver][t1.ver]);
}
 
// dissolve()  a bond (from one side).
 
inline void tetgenmesh::dissolve(triface &t) { t.tet[t.ver & 3] = NULL; }
 
// enext()  finds the next edge (counterclockwise) in the same face.
 
inline void tetgenmesh::enext(triface &t1, triface &t2)
{
  t2.tet = t1.tet;
  t2.ver = enexttbl[t1.ver];
}
 
inline void tetgenmesh::enextself(triface &t) { t.ver = enexttbl[t.ver]; }
 
// eprev()   finds the next edge (clockwise) in the same face.
 
inline void tetgenmesh::eprev(triface &t1, triface &t2)
{
  t2.tet = t1.tet;
  t2.ver = eprevtbl[t1.ver];
}
 
inline void tetgenmesh::eprevself(triface &t) { t.ver = eprevtbl[t.ver]; }
 
// esym()  finds the reversed edge.  It is in the other face of the
//   same tetrahedron.
 
inline void tetgenmesh::esym(triface &t1, triface &t2)
{
  (t2).tet = (t1).tet;
  (t2).ver = esymtbl[(t1).ver];
}
 
inline void tetgenmesh::esymself(triface &t) { (t).ver = esymtbl[(t).ver]; }
 
// enextesym()  finds the reversed edge of the next edge. It is in the other
//   face of the same tetrahedron. It is the combination esym() * enext().
 
inline void tetgenmesh::enextesym(triface &t1, triface &t2)
{
  t2.tet = t1.tet;
  t2.ver = enextesymtbl[t1.ver];
}
 
inline void tetgenmesh::enextesymself(triface &t)
{
  t.ver = enextesymtbl[t.ver];
}
 
// eprevesym()  finds the reversed edge of the previous edge.
 
inline void tetgenmesh::eprevesym(triface &t1, triface &t2)
{
  t2.tet = t1.tet;
  t2.ver = eprevesymtbl[t1.ver];
}
 
inline void tetgenmesh::eprevesymself(triface &t)
{
  t.ver = eprevesymtbl[t.ver];
}
 
// eorgoppo()    Finds the opposite face of the origin of the current edge.
//               Return the opposite edge of the current edge.
 
inline void tetgenmesh::eorgoppo(triface &t1, triface &t2)
{
  t2.tet = t1.tet;
  t2.ver = eorgoppotbl[t1.ver];
}
 
inline void tetgenmesh::eorgoppoself(triface &t) { t.ver = eorgoppotbl[t.ver]; }
 
// edestoppo()    Finds the opposite face of the destination of the current
//                edge. Return the opposite edge of the current edge.
 
inline void tetgenmesh::edestoppo(triface &t1, triface &t2)
{
  t2.tet = t1.tet;
  t2.ver = edestoppotbl[t1.ver];
}
 
inline void tetgenmesh::edestoppoself(triface &t)
{
  t.ver = edestoppotbl[t.ver];
}
 
// fsym()  finds the adjacent tetrahedron at the same face and the same edge.
 
inline void tetgenmesh::fsym(triface &t1, triface &t2)
{
  decode((t1).tet[(t1).ver & 3], t2);
  t2.ver = fsymtbl[t1.ver][t2.ver];
}
 
#define fsymself(t)                                                            \
  t1ver = (t).ver;                                                             \
  decode((t).tet[(t).ver & 3], (t));                                           \
  (t).ver = fsymtbl[t1ver][(t).ver]
 
// fnext()  finds the next face while rotating about an edge according to
//   a right-hand rule. The face is in the adjacent tetrahedron.  It is
//   the combination: fsym() * esym().
 
inline void tetgenmesh::fnext(triface &t1, triface &t2)
{
  decode(t1.tet[facepivot1[t1.ver]], t2);
  t2.ver = facepivot2[t1.ver][t2.ver];
}
 
#define fnextself(t)                                                           \
  t1ver = (t).ver;                                                             \
  decode((t).tet[facepivot1[(t).ver]], (t));                                   \
  (t).ver = facepivot2[t1ver][(t).ver]
 
// The following primtives get or set the origin, destination, face apex,
//   or face opposite of an ordered tetrahedron.
 
inline tetgenmesh::point tetgenmesh::org(triface &t)
{
  return (point)(t).tet[orgpivot[(t).ver]];
}
 
inline tetgenmesh::point tetgenmesh::dest(triface &t)
{
  return (point)(t).tet[destpivot[(t).ver]];
}
 
inline tetgenmesh::point tetgenmesh::apex(triface &t)
{
  return (point)(t).tet[apexpivot[(t).ver]];
}
 
inline tetgenmesh::point tetgenmesh::oppo(triface &t)
{
  return (point)(t).tet[oppopivot[(t).ver]];
}
 
inline void tetgenmesh::setorg(triface &t, point p)
{
  (t).tet[orgpivot[(t).ver]] = (tetrahedron)(p);
}
 
inline void tetgenmesh::setdest(triface &t, point p)
{
  (t).tet[destpivot[(t).ver]] = (tetrahedron)(p);
}
 
inline void tetgenmesh::setapex(triface &t, point p)
{
  (t).tet[apexpivot[(t).ver]] = (tetrahedron)(p);
}
 
inline void tetgenmesh::setoppo(triface &t, point p)
{
  (t).tet[oppopivot[(t).ver]] = (tetrahedron)(p);
}
 
#define setvertices(t, torg, tdest, tapex, toppo)                              \
  (t).tet[orgpivot[(t).ver]] = (tetrahedron)(torg);                            \
  (t).tet[destpivot[(t).ver]] = (tetrahedron)(tdest);                          \
  (t).tet[apexpivot[(t).ver]] = (tetrahedron)(tapex);                          \
  (t).tet[oppopivot[(t).ver]] = (tetrahedron)(toppo)
 
// Check or set a tetrahedron's attributes.
 
inline REAL tetgenmesh::elemattribute(tetrahedron *ptr, int attnum)
{
  return ((REAL *)(ptr))[elemattribindex + attnum];
}
 
inline void tetgenmesh::setelemattribute(tetrahedron *ptr, int attnum,
                                         REAL value)
{
  ((REAL *)(ptr))[elemattribindex + attnum] = value;
}
 
// Check or set a tetrahedron's maximum volume bound.
 
inline REAL tetgenmesh::volumebound(tetrahedron *ptr)
{
  return ((REAL *)(ptr))[volumeboundindex];
}
 
inline void tetgenmesh::setvolumebound(tetrahedron *ptr, REAL value)
{
  ((REAL *)(ptr))[volumeboundindex] = value;
}
 
// Get or set a tetrahedron's index (only used for output).
//    These two routines use the reserved slot ptr[10].
 
inline int tetgenmesh::elemindex(tetrahedron *ptr)
{
  int *iptr = (int *)&(ptr[10]);
  return iptr[0];
}
 
inline void tetgenmesh::setelemindex(tetrahedron *ptr, int value)
{
  int *iptr = (int *)&(ptr[10]);
  iptr[0] = value;
}
 
// Get or set a tetrahedron's marker.
//   Set 'value = 0' cleans all the face/edge flags.
 
inline int tetgenmesh::elemmarker(tetrahedron *ptr)
{
  return ((int *)(ptr))[elemmarkerindex];
}
 
inline void tetgenmesh::setelemmarker(tetrahedron *ptr, int value)
{
  ((int *)(ptr))[elemmarkerindex] = value;
}
 
// infect(), infected(), uninfect() -- primitives to flag or unflag a
//   tetrahedron. The last bit of the element marker is flagged (1)
//   or unflagged (0).
 
inline void tetgenmesh::infect(triface &t)
{
  ((int *)(t.tet))[elemmarkerindex] |= 1;
}
 
inline void tetgenmesh::uninfect(triface &t)
{
  ((int *)(t.tet))[elemmarkerindex] &= ~1;
}
 
inline bool tetgenmesh::infected(triface &t)
{
  return (((int *)(t.tet))[elemmarkerindex] & 1) != 0;
}
 
// marktest(), marktested(), unmarktest() -- primitives to flag or unflag a
//   tetrahedron.  Use the second lowerest bit of the element marker.
 
inline void tetgenmesh::marktest(triface &t)
{
  ((int *)(t.tet))[elemmarkerindex] |= 2;
}
 
inline void tetgenmesh::unmarktest(triface &t)
{
  ((int *)(t.tet))[elemmarkerindex] &= ~2;
}
 
inline bool tetgenmesh::marktested(triface &t)
{
  return (((int *)(t.tet))[elemmarkerindex] & 2) != 0;
}
 
// markface(), unmarkface(), facemarked() -- primitives to flag or unflag a
//   face of a tetrahedron.  From the last 3rd to 6th bits are used for
//   face markers, e.g., the last third bit corresponds to loc = 0.
 
inline void tetgenmesh::markface(triface &t)
{
  ((int *)(t.tet))[elemmarkerindex] |= (4 << (t.ver & 3));
}
 
inline void tetgenmesh::unmarkface(triface &t)
{
  ((int *)(t.tet))[elemmarkerindex] &= ~(4 << (t.ver & 3));
}
 
inline bool tetgenmesh::facemarked(triface &t)
{
  return (((int *)(t.tet))[elemmarkerindex] & (4 << (t.ver & 3))) != 0;
}
 
// markedge(), unmarkedge(), edgemarked() -- primitives to flag or unflag an
//   edge of a tetrahedron.  From the last 7th to 12th bits are used for
//   edge markers, e.g., the last 7th bit corresponds to the 0th edge, etc.
//   Remark: The last 7th bit is marked by 2^6 = 64.
 
inline void tetgenmesh::markedge(triface &t)
{
  ((int *)(t.tet))[elemmarkerindex] |= (int)(64 << ver2edge[(t).ver]);
}
 
inline void tetgenmesh::unmarkedge(triface &t)
{
  ((int *)(t.tet))[elemmarkerindex] &= ~(int)(64 << ver2edge[(t).ver]);
}
 
inline bool tetgenmesh::edgemarked(triface &t)
{
  return (((int *)(t.tet))[elemmarkerindex] & (int)(64 << ver2edge[(t).ver])) !=
         0;
}
 
// marktest2(), unmarktest2(), marktest2ed() -- primitives to flag and unflag
//   a tetrahedron. The 13th bit (2^12 = 4096) is used for this flag.
 
inline void tetgenmesh::marktest2(triface &t)
{
  ((int *)(t.tet))[elemmarkerindex] |= (int)(4096);
}
 
inline void tetgenmesh::unmarktest2(triface &t)
{
  ((int *)(t.tet))[elemmarkerindex] &= ~(int)(4096);
}
 
inline bool tetgenmesh::marktest2ed(triface &t)
{
  return (((int *)(t.tet))[elemmarkerindex] & (int)(4096)) != 0;
}
 
// elemcounter(), setelemcounter() -- primitives to read or ser a (small)
//   integer counter in this tet. It is saved from the 16th bit. On 32 bit
//   system, the range of the counter is [0, 2^15 = 32768].
 
inline int tetgenmesh::elemcounter(triface &t)
{
  return (((int *)(t.tet))[elemmarkerindex]) >> 16;
}
 
inline void tetgenmesh::setelemcounter(triface &t, int value)
{
  int c = ((int *)(t.tet))[elemmarkerindex];
  // Clear the old counter while keep the other flags.
  c &= 65535; // sum_{i=0^15} 2^i
  c |= (value << 16);
  ((int *)(t.tet))[elemmarkerindex] = c;
}
 
inline void tetgenmesh::increaseelemcounter(triface &t)
{
  int c = elemcounter(t);
  setelemcounter(t, c + 1);
}
 
inline void tetgenmesh::decreaseelemcounter(triface &t)
{
  int c = elemcounter(t);
  setelemcounter(t, c - 1);
}
 
// ishulltet()  tests if t is a hull tetrahedron.
 
inline bool tetgenmesh::ishulltet(triface &t)
{
  return (point)(t).tet[7] == dummypoint;
}
 
// isdeadtet()  tests if t is a tetrahedron is dead.
 
inline bool tetgenmesh::isdeadtet(triface &t)
{
  return ((t.tet == NULL) || (t.tet[4] == NULL));
}
 
//                                                                           //
// Primitives for subfaces and subsegments                                   //
//                                                                           //
 
// Each subface contains three pointers to its neighboring subfaces, with
//   edge versions.  To save memory, both information are kept in a single
//   pointer. To make this possible, all subfaces are aligned to eight-byte
//   boundaries, so that the last three bits of each pointer are zeros. An
//   edge version (in the range 0 to 5) is compressed into the last three
//   bits of each pointer by 'sencode()'.  'sdecode()' decodes a pointer,
//   extracting an edge version and a pointer to the beginning of a subface.
 
inline void tetgenmesh::sdecode(shellface sptr, face &s)
{
  s.shver = (int)((uintptr_t)(sptr) & (uintptr_t)7);
  s.sh = (shellface *)((uintptr_t)(sptr) ^ (uintptr_t)(s.shver));
}
 
inline tetgenmesh::shellface tetgenmesh::sencode(face &s)
{
  return (shellface)((uintptr_t)s.sh | (uintptr_t)s.shver);
}
 
inline tetgenmesh::shellface tetgenmesh::sencode2(shellface *sh, int shver)
{
  return (shellface)((uintptr_t)sh | (uintptr_t)shver);
}
 
// sbond() bonds two subfaces (s1) and (s2) together. s1 and s2 must refer
//   to the same edge. No requirement is needed on their orientations.
 
inline void tetgenmesh::sbond(face &s1, face &s2)
{
  s1.sh[s1.shver >> 1] = sencode(s2);
  s2.sh[s2.shver >> 1] = sencode(s1);
}
 
// sbond1() bonds s1 <== s2, i.e., after bonding, s1 is pointing to s2,
//   but s2 is not pointing to s1.  s1 and s2 must refer to the same edge.
//   No requirement is needed on their orientations.
 
inline void tetgenmesh::sbond1(face &s1, face &s2)
{
  s1.sh[s1.shver >> 1] = sencode(s2);
}
 
// Dissolve a subface bond (from one side).  Note that the other subface
//   will still think it's connected to this subface.
 
inline void tetgenmesh::sdissolve(face &s) { s.sh[s.shver >> 1] = NULL; }
 
// spivot() finds the adjacent subface (s2) for a given subface (s1).
//   s1 and s2 share at the same edge.
 
inline void tetgenmesh::spivot(face &s1, face &s2)
{
  shellface sptr = s1.sh[s1.shver >> 1];
  sdecode(sptr, s2);
}
 
inline void tetgenmesh::spivotself(face &s)
{
  shellface sptr = s.sh[s.shver >> 1];
  sdecode(sptr, s);
}
 
// These primitives determine or set the origin, destination, or apex
//   of a subface with respect to the edge version.
 
inline tetgenmesh::point tetgenmesh::sorg(face &s)
{
  return (point)s.sh[sorgpivot[s.shver]];
}
 
inline tetgenmesh::point tetgenmesh::sdest(face &s)
{
  return (point)s.sh[sdestpivot[s.shver]];
}
 
inline tetgenmesh::point tetgenmesh::sapex(face &s)
{
  return (point)s.sh[sapexpivot[s.shver]];
}
 
inline void tetgenmesh::setsorg(face &s, point pointptr)
{
  s.sh[sorgpivot[s.shver]] = (shellface)pointptr;
}
 
inline void tetgenmesh::setsdest(face &s, point pointptr)
{
  s.sh[sdestpivot[s.shver]] = (shellface)pointptr;
}
 
inline void tetgenmesh::setsapex(face &s, point pointptr)
{
  s.sh[sapexpivot[s.shver]] = (shellface)pointptr;
}
 
#define setshvertices(s, pa, pb, pc)                                           \
  setsorg(s, pa);                                                              \
  setsdest(s, pb);                                                             \
  setsapex(s, pc)
 
// sesym()  reserves the direction of the lead edge.
 
inline void tetgenmesh::sesym(face &s1, face &s2)
{
  s2.sh = s1.sh;
  s2.shver = (s1.shver ^ 1); // Inverse the last bit.
}
 
inline void tetgenmesh::sesymself(face &s) { s.shver ^= 1; }
 
// senext()  finds the next edge (counterclockwise) in the same orientation
//   of this face.
 
inline void tetgenmesh::senext(face &s1, face &s2)
{
  s2.sh = s1.sh;
  s2.shver = snextpivot[s1.shver];
}
 
inline void tetgenmesh::senextself(face &s) { s.shver = snextpivot[s.shver]; }
 
inline void tetgenmesh::senext2(face &s1, face &s2)
{
  s2.sh = s1.sh;
  s2.shver = snextpivot[snextpivot[s1.shver]];
}
 
inline void tetgenmesh::senext2self(face &s)
{
  s.shver = snextpivot[snextpivot[s.shver]];
}
 
// Check or set a subface's maximum area bound.
 
inline REAL tetgenmesh::areabound(face &s)
{
  return ((REAL *)(s.sh))[areaboundindex];
}
 
inline void tetgenmesh::setareabound(face &s, REAL value)
{
  ((REAL *)(s.sh))[areaboundindex] = value;
}
 
// These two primitives read or set a shell marker.  Shell markers are used
//   to hold user boundary information.
 
inline int tetgenmesh::shellmark(face &s)
{
  return ((int *)(s.sh))[shmarkindex];
}
 
inline void tetgenmesh::setshellmark(face &s, int value)
{
  ((int *)(s.sh))[shmarkindex] = value;
}
 
// sinfect(), sinfected(), suninfect() -- primitives to flag or unflag a
//   subface. The last bit of ((int *) ((s).sh))[shmarkindex+1] is flagged.
 
inline void tetgenmesh::sinfect(face &s)
{
  ((int *)((s).sh))[shmarkindex + 1] =
    (((int *)((s).sh))[shmarkindex + 1] | (int)1);
}
 
inline void tetgenmesh::suninfect(face &s)
{
  ((int *)((s).sh))[shmarkindex + 1] =
    (((int *)((s).sh))[shmarkindex + 1] & ~(int)1);
}
 
// Test a subface for viral infection.
 
inline bool tetgenmesh::sinfected(face &s)
{
  return (((int *)((s).sh))[shmarkindex + 1] & (int)1) != 0;
}
 
// smarktest(), smarktested(), sunmarktest() -- primitives to flag or unflag
//   a subface. The last 2nd bit of the integer is flagged.
 
inline void tetgenmesh::smarktest(face &s)
{
  ((int *)((s).sh))[shmarkindex + 1] =
    (((int *)((s).sh))[shmarkindex + 1] | (int)2);
}
 
inline void tetgenmesh::sunmarktest(face &s)
{
  ((int *)((s).sh))[shmarkindex + 1] =
    (((int *)((s).sh))[shmarkindex + 1] & ~(int)2);
}
 
inline bool tetgenmesh::smarktested(face &s)
{
  return ((((int *)((s).sh))[shmarkindex + 1] & (int)2) != 0);
}
 
// smarktest2(), smarktest2ed(), sunmarktest2() -- primitives to flag or
//   unflag a subface. The last 3rd bit of the integer is flagged.
 
inline void tetgenmesh::smarktest2(face &s)
{
  ((int *)((s).sh))[shmarkindex + 1] =
    (((int *)((s).sh))[shmarkindex + 1] | (int)4);
}
 
inline void tetgenmesh::sunmarktest2(face &s)
{
  ((int *)((s).sh))[shmarkindex + 1] =
    (((int *)((s).sh))[shmarkindex + 1] & ~(int)4);
}
 
inline bool tetgenmesh::smarktest2ed(face &s)
{
  return ((((int *)((s).sh))[shmarkindex + 1] & (int)4) != 0);
}
 
// The last 4th bit of ((int *) ((s).sh))[shmarkindex+1] is flagged.
 
inline void tetgenmesh::smarktest3(face &s)
{
  ((int *)((s).sh))[shmarkindex + 1] =
    (((int *)((s).sh))[shmarkindex + 1] | (int)8);
}
 
inline void tetgenmesh::sunmarktest3(face &s)
{
  ((int *)((s).sh))[shmarkindex + 1] =
    (((int *)((s).sh))[shmarkindex + 1] & ~(int)8);
}
 
inline bool tetgenmesh::smarktest3ed(face &s)
{
  return ((((int *)((s).sh))[shmarkindex + 1] & (int)8) != 0);
}
 
// Each facet has a unique index (automatically indexed). Starting from '0'.
// We save this index in the same field of the shell type.
 
inline void tetgenmesh::setfacetindex(face &s, int value)
{
  ((int *)(s.sh))[shmarkindex + 2] = value;
}
 
inline int tetgenmesh::getfacetindex(face &s)
{
  return ((int *)(s.sh))[shmarkindex + 2];
}
 
//                                                                           //
// Primitives for interacting between tetrahedra and subfaces                //
//                                                                           //
 
// tsbond() bond a tetrahedron (t) and a subface (s) together.
// Note that t and s must be the same face and the same edge. Moreover,
//   t and s have the same orientation.
// Since the edge number in t and in s can be any number in {0,1,2}. We bond
//   the edge in s which corresponds to t's 0th edge, and vice versa.
 
inline void tetgenmesh::tsbond(triface &t, face &s)
{
  if((t).tet[9] == NULL) {
    // Allocate space for this tet.
    (t).tet[9] = (tetrahedron)tet2subpool->alloc();
    // Initialize.
    for(int i = 0; i < 4; i++) { ((shellface *)(t).tet[9])[i] = NULL; }
  }
  // Bond t <== s.
  ((shellface *)(t).tet[9])[(t).ver & 3] =
    sencode2((s).sh, tsbondtbl[t.ver][s.shver]);
  // Bond s <== t.
  s.sh[9 + ((s).shver & 1)] =
    (shellface)encode2((t).tet, stbondtbl[t.ver][s.shver]);
}
 
// tspivot() finds a subface (s) abutting on the given tetrahdera (t).
//   Return s.sh = NULL if there is no subface at t. Otherwise, return
//   the subface s, and s and t must be at the same edge wth the same
//   orientation.
 
inline void tetgenmesh::tspivot(triface &t, face &s)
{
  if((t).tet[9] == NULL) {
    (s).sh = NULL;
    return;
  }
  // Get the attached subface s.
  sdecode(((shellface *)(t).tet[9])[(t).ver & 3], (s));
  (s).shver = tspivottbl[t.ver][s.shver];
}
 
// Quickly check if the handle (t, v) is a subface.
#define issubface(t) ((t).tet[9] && ((t).tet[9])[(t).ver & 3])
 
// stpivot() finds a tetrahedron (t) abutting a given subface (s).
//   Return the t (if it exists) with the same edge and the same
//   orientation of s.
 
inline void tetgenmesh::stpivot(face &s, triface &t)
{
  decode((tetrahedron)s.sh[9 + (s.shver & 1)], t);
  if((t).tet == NULL) { return; }
  (t).ver = stpivottbl[t.ver][s.shver];
}
 
// Quickly check if this subface is attached to a tetrahedron.
 
#define isshtet(s) ((s).sh[9 + ((s).shver & 1)])
 
// tsdissolve() dissolve a bond (from the tetrahedron side).
 
inline void tetgenmesh::tsdissolve(triface &t)
{
  if((t).tet[9] != NULL) { ((shellface *)(t).tet[9])[(t).ver & 3] = NULL; }
}
 
// stdissolve() dissolve a bond (from the subface side).
 
inline void tetgenmesh::stdissolve(face &s)
{
  (s).sh[9] = NULL;
  (s).sh[10] = NULL;
}
 
//                                                                           //
// Primitives for interacting between subfaces and segments                  //
//                                                                           //
 
// ssbond() bond a subface to a subsegment.
 
inline void tetgenmesh::ssbond(face &s, face &edge)
{
  s.sh[6 + (s.shver >> 1)] = sencode(edge);
  edge.sh[0] = sencode(s);
}
 
inline void tetgenmesh::ssbond1(face &s, face &edge)
{
  s.sh[6 + (s.shver >> 1)] = sencode(edge);
  // edge.sh[0] = sencode(s);
}
 
// ssdisolve() dissolve a bond (from the subface side)
 
inline void tetgenmesh::ssdissolve(face &s) { s.sh[6 + (s.shver >> 1)] = NULL; }
 
// sspivot() finds a subsegment abutting a subface.
 
inline void tetgenmesh::sspivot(face &s, face &edge)
{
  sdecode((shellface)s.sh[6 + (s.shver >> 1)], edge);
}
 
// Quickly check if the edge is a subsegment.
 
#define isshsubseg(s) ((s).sh[6 + ((s).shver >> 1)])
 
//                                                                           //
// Primitives for interacting between tetrahedra and segments                //
//                                                                           //
 
inline void tetgenmesh::tssbond1(triface &t, face &s)
{
  if((t).tet[8] == NULL) {
    // Allocate space for this tet.
    (t).tet[8] = (tetrahedron)tet2segpool->alloc();
    // Initialization.
    for(int i = 0; i < 6; i++) { ((shellface *)(t).tet[8])[i] = NULL; }
  }
  ((shellface *)(t).tet[8])[ver2edge[(t).ver]] = sencode((s));
}
 
inline void tetgenmesh::sstbond1(face &s, triface &t)
{
  ((tetrahedron *)(s).sh)[9] = encode(t);
}
 
inline void tetgenmesh::tssdissolve1(triface &t)
{
  if((t).tet[8] != NULL) {
    ((shellface *)(t).tet[8])[ver2edge[(t).ver]] = NULL;
  }
}
 
inline void tetgenmesh::sstdissolve1(face &s)
{
  ((tetrahedron *)(s).sh)[9] = NULL;
}
 
inline void tetgenmesh::tsspivot1(triface &t, face &s)
{
  if((t).tet[8] != NULL) {
    sdecode(((shellface *)(t).tet[8])[ver2edge[(t).ver]], s);
  }
  else {
    (s).sh = NULL;
  }
}
 
// Quickly check whether 't' is a segment or not.
 
#define issubseg(t) ((t).tet[8] && ((t).tet[8])[ver2edge[(t).ver]])
 
inline void tetgenmesh::sstpivot1(face &s, triface &t)
{
  decode((tetrahedron)s.sh[9], t);
}
 
//                                                                           //
// Primitives for points                                                     //
//                                                                           //
 
inline int tetgenmesh::pointmark(point pt)
{
  return ((int *)(pt))[pointmarkindex];
}
 
inline void tetgenmesh::setpointmark(point pt, int value)
{
  ((int *)(pt))[pointmarkindex] = value;
}
 
// These two primitives set and read the type of the point.
 
inline enum tetgenmesh::verttype tetgenmesh::pointtype(point pt)
{
  return (enum verttype)(((int *)(pt))[pointmarkindex + 1] >> (int)8);
}
 
inline void tetgenmesh::setpointtype(point pt, enum verttype value)
{
  ((int *)(pt))[pointmarkindex + 1] =
    ((int)value << 8) + (((int *)(pt))[pointmarkindex + 1] & (int)255);
}
 
// Read and set the geometry tag of the point (used by -s option).
 
inline int tetgenmesh::pointgeomtag(point pt)
{
  return ((int *)(pt))[pointmarkindex + 2];
}
 
inline void tetgenmesh::setpointgeomtag(point pt, int value)
{
  ((int *)(pt))[pointmarkindex + 2] = value;
}
 
// Read and set the u,v coordinates of the point (used by -s option).
 
inline REAL tetgenmesh::pointgeomuv(point pt, int i)
{
  return pt[pointparamindex + i];
}
 
inline void tetgenmesh::setpointgeomuv(point pt, int i, REAL value)
{
  pt[pointparamindex + i] = value;
}
 
// pinfect(), puninfect(), pinfected() -- primitives to flag or unflag
//   a point. The last bit of the integer '[pointindex+1]' is flagged.
 
inline void tetgenmesh::pinfect(point pt)
{
  ((int *)(pt))[pointmarkindex + 1] |= (int)1;
}
 
inline void tetgenmesh::puninfect(point pt)
{
  ((int *)(pt))[pointmarkindex + 1] &= ~(int)1;
}
 
inline bool tetgenmesh::pinfected(point pt)
{
  return (((int *)(pt))[pointmarkindex + 1] & (int)1) != 0;
}
 
// pmarktest(), punmarktest(), pmarktested() -- more primitives to
//   flag or unflag a point.
 
inline void tetgenmesh::pmarktest(point pt)
{
  ((int *)(pt))[pointmarkindex + 1] |= (int)2;
}
 
inline void tetgenmesh::punmarktest(point pt)
{
  ((int *)(pt))[pointmarkindex + 1] &= ~(int)2;
}
 
inline bool tetgenmesh::pmarktested(point pt)
{
  return (((int *)(pt))[pointmarkindex + 1] & (int)2) != 0;
}
 
inline void tetgenmesh::pmarktest2(point pt)
{
  ((int *)(pt))[pointmarkindex + 1] |= (int)4;
}
 
inline void tetgenmesh::punmarktest2(point pt)
{
  ((int *)(pt))[pointmarkindex + 1] &= ~(int)4;
}
 
inline bool tetgenmesh::pmarktest2ed(point pt)
{
  return (((int *)(pt))[pointmarkindex + 1] & (int)4) != 0;
}
 
inline void tetgenmesh::pmarktest3(point pt)
{
  ((int *)(pt))[pointmarkindex + 1] |= (int)8;
}
 
inline void tetgenmesh::punmarktest3(point pt)
{
  ((int *)(pt))[pointmarkindex + 1] &= ~(int)8;
}
 
inline bool tetgenmesh::pmarktest3ed(point pt)
{
  return (((int *)(pt))[pointmarkindex + 1] & (int)8) != 0;
}
 
// These following primitives set and read a pointer to a tetrahedron
//   a subface/subsegment, a point, or a tet of background mesh.
 
inline tetgenmesh::tetrahedron tetgenmesh::point2tet(point pt)
{
  return ((tetrahedron *)(pt))[point2simindex];
}
 
inline void tetgenmesh::setpoint2tet(point pt, tetrahedron value)
{
  ((tetrahedron *)(pt))[point2simindex] = value;
}
 
inline tetgenmesh::point tetgenmesh::point2ppt(point pt)
{
  return (point)((tetrahedron *)(pt))[point2simindex + 1];
}
 
inline void tetgenmesh::setpoint2ppt(point pt, point value)
{
  ((tetrahedron *)(pt))[point2simindex + 1] = (tetrahedron)value;
}
 
inline tetgenmesh::shellface tetgenmesh::point2sh(point pt)
{
  return (shellface)((tetrahedron *)(pt))[point2simindex + 2];
}
 
inline void tetgenmesh::setpoint2sh(point pt, shellface value)
{
  ((tetrahedron *)(pt))[point2simindex + 2] = (tetrahedron)value;
}
 
inline tetgenmesh::tetrahedron tetgenmesh::point2bgmtet(point pt)
{
  return ((tetrahedron *)(pt))[point2simindex + 3];
}
 
inline void tetgenmesh::setpoint2bgmtet(point pt, tetrahedron value)
{
  ((tetrahedron *)(pt))[point2simindex + 3] = value;
}
 
// The primitives for saving and getting the insertion radius.
inline void tetgenmesh::setpointinsradius(point pt, REAL value)
{
  pt[pointinsradiusindex] = value;
}
 
inline REAL tetgenmesh::getpointinsradius(point pt)
{
  return pt[pointinsradiusindex];
}
 
inline bool tetgenmesh::issteinerpoint(point pt)
{
  return (pointtype(pt) == FREESEGVERTEX) ||
         (pointtype(pt) == FREEFACETVERTEX) || (pointtype(pt) == FREEVOLVERTEX);
}
 
// point2tetorg()    Get the tetrahedron whose origin is the point.
 
inline void tetgenmesh::point2tetorg(point pa, triface &searchtet)
{
  decode(point2tet(pa), searchtet);
  if((point)searchtet.tet[4] == pa) { searchtet.ver = 11; }
  else if((point)searchtet.tet[5] == pa) {
    searchtet.ver = 3;
  }
  else if((point)searchtet.tet[6] == pa) {
    searchtet.ver = 7;
  }
  else {
    searchtet.ver = 0;
  }
}
 
// point2shorg()    Get the subface/segment whose origin is the point.
 
inline void tetgenmesh::point2shorg(point pa, face &searchsh)
{
  sdecode(point2sh(pa), searchsh);
  if((point)searchsh.sh[3] == pa) { searchsh.shver = 0; }
  else if((point)searchsh.sh[4] == pa) {
    searchsh.shver = (searchsh.sh[5] != NULL ? 2 : 1);
  }
  else {
    searchsh.shver = 4;
  }
}
 
// farsorg()    Return the origin of the subsegment.
// farsdest()   Return the destination of the subsegment.
 
inline tetgenmesh::point tetgenmesh::farsorg(face &s)
{
  face travesh, neighsh;
 
  travesh = s;
  while(1) {
    senext2(travesh, neighsh);
    spivotself(neighsh);
    if(neighsh.sh == NULL) break;
    if(sorg(neighsh) != sorg(travesh)) sesymself(neighsh);
    senext2(neighsh, travesh);
  }
  return sorg(travesh);
}
 
inline tetgenmesh::point tetgenmesh::farsdest(face &s)
{
  face travesh, neighsh;
 
  travesh = s;
  while(1) {
    senext(travesh, neighsh);
    spivotself(neighsh);
    if(neighsh.sh == NULL) break;
    if(sdest(neighsh) != sdest(travesh)) sesymself(neighsh);
    senext(neighsh, travesh);
  }
  return sdest(travesh);
}
 
//                                                                           //
// Linear algebra operators.                                                 //
//                                                                           //
 
// dot() returns the dot product: v1 dot v2.
inline REAL tetgenmesh::dot(REAL *v1, REAL *v2)
{
  return v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2];
}
 
// cross() computes the cross product: n = v1 cross v2.
inline void tetgenmesh::cross(REAL *v1, REAL *v2, REAL *n)
{
  n[0] = v1[1] * v2[2] - v2[1] * v1[2];
  n[1] = -(v1[0] * v2[2] - v2[0] * v1[2]);
  n[2] = v1[0] * v2[1] - v2[0] * v1[1];
}
 
// distance() computes the Euclidean distance between two points.
inline REAL tetgenmesh::distance(REAL *p1, REAL *p2)
{
  return sqrt((p2[0] - p1[0]) * (p2[0] - p1[0]) +
              (p2[1] - p1[1]) * (p2[1] - p1[1]) +
              (p2[2] - p1[2]) * (p2[2] - p1[2]));
}
 
inline REAL tetgenmesh::norm2(REAL x, REAL y, REAL z)
{
  return (x) * (x) + (y) * (y) + (z) * (z);
}
 
#endif // tetgenBRH