glukoz/blender
0
0
Fork 0
forked from blender/blender
Code
softbody-stable-v1
blender/extern/tetgen/tetgen.h
over0219 b86af3ecb7 added tetgen
2020-06-08 14:38:27 -05:00

3341 lines
132 KiB
C++
Raw Permalink Blame History

///////////////////////////////////////////////////////////////////////////////
// //
// TetGen //
// //
// A Quality Tetrahedral Mesh Generator and A 3D Delaunay Triangulator //
// //
// Version 1.5 //
// November 4, 2013 //
// //
// 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. //
// //
///////////////////////////////////////////////////////////////////////////////
#ifndef tetgenH
#define tetgenH
#ifndef TETGEN_UNUSED
#define TETGEN_UNUSED(x) (void)(x)
#endif
#define TETLIBRARY
// To compile TetGen as a library instead of an executable program, define
// the TETLIBRARY symbol.
// #define TETLIBRARY
// Uncomment the following line to disable assert macros. These macros were
// inserted in the code where I hoped to catch bugs. They may slow down the
// speed of TetGen.
// #define NDEBUG
// TetGen default uses the double precision (64 bit) for a real number.
// Alternatively, one can use the single precision (32 bit) 'float' if the
// memory is limited.
#define REAL double // #define REAL float
// Maximum number of characters in a file name (including the null).
#define FILENAMESIZE 1024
// Maximum number of chars in a line read from a file (including the null).
#define INPUTLINESIZE 2048
// TetGen only uses the C standard library.
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <assert.h>
// The types 'intptr_t' and 'uintptr_t' are signed and unsigned integer types,
// respectively. They are guaranteed to be the same width as a pointer.
// They are defined in <stdint.h> by the C99 Standard. However, Microsoft
// Visual C++ 2003 -- 2008 (Visual C++ 7.1 - 9) doesn't ship with this header
// file. In such case, we can define them by ourself.
// Update (learned from Stack Overflow): Visual Studio 2010 and Visual C++ 2010
// Express both have stdint.h
// The following piece of code was provided by Steven Johnson (MIT). Define the
// symbol _MSC_VER if you are using Microsoft Visual C++. Moreover, define
// the _WIN64 symbol if you are running TetGen on Win64 systems.
#ifdef _MSC_VER // Microsoft Visual C++
# ifdef _WIN64
typedef __int64 intptr_t;
typedef unsigned __int64 uintptr_t;
# else // not _WIN64
typedef int intptr_t;
typedef unsigned int uintptr_t;
# endif
#else // not Visual C++
# include <stdint.h>
#endif
///////////////////////////////////////////////////////////////////////////////
// //
// tetgenio //
// //
// A structure for transferring data into and out of TetGen's mesh structure,//
// 'tetgenmesh' (declared below). //
// //
// The input of TetGen is either a 3D point set, or a 3D piecewise linear //
// complex (PLC), or a tetrahedral mesh. Depending on the input object and //
// the specified options, the output of TetGen is either a Delaunay (or wei- //
// ghted Delaunay) tetrahedralization, or a constrained (Delaunay) tetrahed- //
// ralization, or a quality tetrahedral mesh. //
// //
// A piecewise linear complex (PLC) represents a 3D polyhedral domain with //
// possibly internal boundaries(subdomains). It is introduced in [Miller et //
// al, 1996]. Basically it is a set of "cells", i.e., vertices, edges, poly- //
// gons, and polyhedra, and the intersection of any two of its cells is the //
// union of other cells of it. //
// //
// TetGen uses a set of files to describe the inputs and outputs. Each file //
// is identified from its file extension (.node, .ele, .face, .edge, etc). //
// //
// The 'tetgenio' structure is a collection of arrays of data, i.e., points, //
// facets, tetrahedra, and so forth. It contains functions to read and write //
// (input and output) files of TetGen as well as other supported mesh files. //
// //
// Once an object of tetgenio is declared, no array is created. One has to //
// allocate enough memory for them. On deletion of this object, the memory //
// occupied by these arrays needs to be freed. The routine deinitialize() //
// will be automatically called. It frees the memory for an array if it is //
// not a NULL. Note that it assumes that the memory is allocated by the C++ //
// "new" operator. Otherwise, the user is responsible to free them and all //
// pointers must be NULL before the call of the destructor. //
// //
///////////////////////////////////////////////////////////////////////////////
class tetgenio {
public:
// A "polygon" describes a simple polygon (no holes). It is not necessarily
// convex. Each polygon contains a number of corners (points) and the same
// number of sides (edges). The points of the polygon must be given in
// either counterclockwise or clockwise order and they form a ring, so
// every two consecutive points forms an edge of the polygon.
typedef struct {
int *vertexlist;
int numberofvertices;
} polygon;
// A "facet" describes a polygonal region possibly with holes, edges, and
// points floating in it. Each facet consists of a list of polygons and
// a list of hole points (which lie strictly inside holes).
typedef struct {
polygon *polygonlist;
int numberofpolygons;
REAL *holelist;
int numberofholes;
} facet;
// A "voroedge" is an edge of the Voronoi diagram. It corresponds to a
// Delaunay face. Each voroedge is either a line segment connecting
// two Voronoi vertices or a ray starting from a Voronoi vertex to an
// "infinite vertex". 'v1' and 'v2' are two indices pointing to the
// list of Voronoi vertices. 'v1' must be non-negative, while 'v2' may
// be -1 if it is a ray, in this case, the unit normal of this ray is
// given in 'vnormal'.
typedef struct {
int v1, v2;
REAL vnormal[3];
} voroedge;
// A "vorofacet" is an facet of the Voronoi diagram. It corresponds to a
// Delaunay edge. Each Voronoi facet is a convex polygon formed by a
// list of Voronoi edges, it may not be closed. 'c1' and 'c2' are two
// indices pointing into the list of Voronoi cells, i.e., the two cells
// share this facet. 'elist' is an array of indices pointing into the
// list of Voronoi edges, 'elist[0]' saves the number of Voronoi edges
// (including rays) of this facet.
typedef struct {
int c1, c2;
int *elist;
} vorofacet;
// Additional parameters associated with an input (or mesh) vertex.
// These informations are provided by CAD libraries.
typedef struct {
REAL uv[2];
int tag;
int type; // 0, 1, or 2.
} pointparam;
// Callback functions for meshing PSCs.
typedef REAL (* GetVertexParamOnEdge)(void*, int, int);
typedef void (* GetSteinerOnEdge)(void*, int, REAL, REAL*);
typedef void (* GetVertexParamOnFace)(void*, int, int, REAL*);
typedef void (* GetEdgeSteinerParamOnFace)(void*, int, REAL, int, REAL*);
typedef void (* GetSteinerOnFace)(void*, int, REAL*, REAL*);
// A callback function for mesh refinement.
typedef bool (* TetSizeFunc)(REAL*, REAL*, REAL*, REAL*, REAL*, REAL);
// Items are numbered starting from 'firstnumber' (0 or 1), default is 0.
int firstnumber;
// Dimension of the mesh (2 or 3), default is 3.
int mesh_dim;
// Does the lines in .node file contain index or not, default is 1.
int useindex;
// 'pointlist': An array of point coordinates. The first point's x
// coordinate is at index [0] and its y coordinate at index [1], its
// z coordinate is at index [2], followed by the coordinates of the
// remaining points. Each point occupies three REALs.
// 'pointattributelist': An array of point attributes. Each point's
// attributes occupy 'numberofpointattributes' REALs.
// 'pointmtrlist': An array of metric tensors at points. Each point's
// tensor occupies 'numberofpointmtr' REALs.
// 'pointmarkerlist': An array of point markers; one integer per point.
REAL *pointlist;
REAL *pointattributelist;
REAL *pointmtrlist;
int *pointmarkerlist;
pointparam *pointparamlist;
int numberofpoints;
int numberofpointattributes;
int numberofpointmtrs;
// 'tetrahedronlist': An array of tetrahedron corners. The first
// tetrahedron's first corner is at index [0], followed by its other
// corners, followed by six nodes on the edges of the tetrahedron if the
// second order option (-o2) is applied. Each tetrahedron occupies
// 'numberofcorners' ints. The second order nodes are ouput only.
// 'tetrahedronattributelist': An array of tetrahedron attributes. Each
// tetrahedron's attributes occupy 'numberoftetrahedronattributes' REALs.
// 'tetrahedronvolumelist': An array of constraints, i.e. tetrahedron's
// volume; one REAL per element. Input only.
// 'neighborlist': An array of tetrahedron neighbors; 4 ints per element.
// Output only.
int *tetrahedronlist;
REAL *tetrahedronattributelist;
REAL *tetrahedronvolumelist;
int *neighborlist;
int numberoftetrahedra;
int numberofcorners;
int numberoftetrahedronattributes;
// 'facetlist': An array of facets. Each entry is a structure of facet.
// 'facetmarkerlist': An array of facet markers; one int per facet.
facet *facetlist;
int *facetmarkerlist;
int numberoffacets;
// 'holelist': An array of holes (in volume). Each hole is given by a
// seed (point) which lies strictly inside it. The first seed's x, y and z
// coordinates are at indices [0], [1] and [2], followed by the
// remaining seeds. Three REALs per hole.
REAL *holelist;
int numberofholes;
// 'regionlist': An array of regions (subdomains). Each region is given by
// a seed (point) which lies strictly inside it. The first seed's x, y and
// z coordinates are at indices [0], [1] and [2], followed by the regional
// attribute at index [3], followed by the maximum volume at index [4].
// Five REALs per region.
// Note that each regional attribute is used only if you select the 'A'
// switch, and each volume constraint is used only if you select the
// 'a' switch (with no number following).
REAL *regionlist;
int numberofregions;
// 'facetconstraintlist': An array of facet constraints. Each constraint
// specifies a maximum area bound on the subfaces of that facet. The
// first facet constraint is given by a facet marker at index [0] and its
// maximum area bound at index [1], followed by the remaining facet con-
// straints. Two REALs per facet constraint. Note: the facet marker is
// actually an integer.
REAL *facetconstraintlist;
int numberoffacetconstraints;
// 'segmentconstraintlist': An array of segment constraints. Each constraint
// specifies a maximum length bound on the subsegments of that segment.
// The first constraint is given by the two endpoints of the segment at
// index [0] and [1], and the maximum length bound at index [2], followed
// by the remaining segment constraints. Three REALs per constraint.
// Note the segment endpoints are actually integers.
REAL *segmentconstraintlist;
int numberofsegmentconstraints;
// 'trifacelist': An array of face (triangle) corners. The first face's
// three corners are at indices [0], [1] and [2], followed by the remaining
// faces. Three ints per face.
// 'trifacemarkerlist': An array of face markers; one int per face.
// 'o2facelist': An array of second order nodes (on the edges) of the face.
// It is output only if the second order option (-o2) is applied. The
// first face's three second order nodes are at [0], [1], and [2],
// followed by the remaining faces. Three ints per face.
// 'adjtetlist': An array of adjacent tetrahedra to the faces. The first
// face's two adjacent tetrahedra are at indices [0] and [1], followed by
// the remaining faces. A '-1' indicates outside (no adj. tet). This list
// is output when '-nn' switch is used. Output only.
int *trifacelist;
int *trifacemarkerlist;
int *o2facelist;
int *adjtetlist;
int numberoftrifaces;
// 'edgelist': An array of edge endpoints. The first edge's endpoints
// are at indices [0] and [1], followed by the remaining edges.
// Two ints per edge.
// 'edgemarkerlist': An array of edge markers; one int per edge.
// 'o2edgelist': An array of midpoints of edges. It is output only if the
// second order option (-o2) is applied. One int per edge.
// 'edgeadjtetlist': An array of adjacent tetrahedra to the edges. One
// tetrahedron (an integer) per edge.
int *edgelist;
int *edgemarkerlist;
int *o2edgelist;
int *edgeadjtetlist;
int numberofedges;
// 'vpointlist': An array of Voronoi vertex coordinates (like pointlist).
// 'vedgelist': An array of Voronoi edges. Each entry is a 'voroedge'.
// 'vfacetlist': An array of Voronoi facets. Each entry is a 'vorofacet'.
// 'vcelllist': An array of Voronoi cells. Each entry is an array of
// indices pointing into 'vfacetlist'. The 0th entry is used to store
// the length of this array.
REAL *vpointlist;
voroedge *vedgelist;
vorofacet *vfacetlist;
int **vcelllist;
int numberofvpoints;
int numberofvedges;
int numberofvfacets;
int numberofvcells;
// Variable (and callback functions) for meshing PSCs.
void *geomhandle;
GetVertexParamOnEdge getvertexparamonedge;
GetSteinerOnEdge getsteineronedge;
GetVertexParamOnFace getvertexparamonface;
GetEdgeSteinerParamOnFace getedgesteinerparamonface;
GetSteinerOnFace getsteineronface;
// A callback function.
TetSizeFunc tetunsuitable;
// Input & output routines.
bool load_node_call(FILE* infile, int markers, int uvflag, char*);
bool load_node(char*);
bool load_edge(char*);
bool load_face(char*);
bool load_tet(char*);
bool load_vol(char*);
bool load_var(char*);
bool load_mtr(char*);
bool load_pbc(char*);
bool load_poly(char*);
bool load_off(char*);
bool load_ply(char*);
bool load_stl(char*);
bool load_vtk(char*);
bool load_medit(char*, int);
bool load_plc(char*, int);
bool load_tetmesh(char*, int);
void save_nodes(char*);
void save_elements(char*);
void save_faces(char*);
void save_edges(char*);
void save_neighbors(char*);
void save_poly(char*);
void save_faces2smesh(char*);
// Read line and parse string functions.
char *readline(char* string, FILE* infile, int *linenumber);
char *findnextfield(char* string);
char *readnumberline(char* string, FILE* infile, char* infilename);
char *findnextnumber(char* string);
static void init(polygon* p) {
p->vertexlist = (int *) NULL;
p->numberofvertices = 0;
}
static void init(facet* f) {
f->polygonlist = (polygon *) NULL;
f->numberofpolygons = 0;
f->holelist = (REAL *) NULL;
f->numberofholes = 0;
}
// Initialize routine.
void initialize()
{
firstnumber = 0;
mesh_dim = 3;
useindex = 1;
pointlist = (REAL *) NULL;
pointattributelist = (REAL *) NULL;
pointmtrlist = (REAL *) NULL;
pointmarkerlist = (int *) NULL;
pointparamlist = (pointparam *) NULL;
numberofpoints = 0;
numberofpointattributes = 0;
numberofpointmtrs = 0;
tetrahedronlist = (int *) NULL;
tetrahedronattributelist = (REAL *) NULL;
tetrahedronvolumelist = (REAL *) NULL;
neighborlist = (int *) NULL;
numberoftetrahedra = 0;
numberofcorners = 4;
numberoftetrahedronattributes = 0;
trifacelist = (int *) NULL;
trifacemarkerlist = (int *) NULL;
o2facelist = (int *) NULL;
adjtetlist = (int *) NULL;
numberoftrifaces = 0;
edgelist = (int *) NULL;
edgemarkerlist = (int *) NULL;
o2edgelist = (int *) NULL;
edgeadjtetlist = (int *) NULL;
numberofedges = 0;
facetlist = (facet *) NULL;
facetmarkerlist = (int *) NULL;
numberoffacets = 0;
holelist = (REAL *) NULL;
numberofholes = 0;
regionlist = (REAL *) NULL;
numberofregions = 0;
facetconstraintlist = (REAL *) NULL;
numberoffacetconstraints = 0;
segmentconstraintlist = (REAL *) NULL;
numberofsegmentconstraints = 0;
vpointlist = (REAL *) NULL;
vedgelist = (voroedge *) NULL;
vfacetlist = (vorofacet *) NULL;
vcelllist = (int **) NULL;
numberofvpoints = 0;
numberofvedges = 0;
numberofvfacets = 0;
numberofvcells = 0;
tetunsuitable = NULL;
geomhandle = NULL;
getvertexparamonedge = NULL;
getsteineronedge = NULL;
getvertexparamonface = NULL;
getedgesteinerparamonface = NULL;
getsteineronface = NULL;
}
// Free the memory allocated in 'tetgenio'. Note that it assumes that the
// memory was allocated by the "new" operator (C++).
void deinitialize()
{
int i, j;
if (pointlist != (REAL *) NULL) {
delete [] pointlist;
}
if (pointattributelist != (REAL *) NULL) {
delete [] pointattributelist;
}
if (pointmtrlist != (REAL *) NULL) {
delete [] pointmtrlist;
}
if (pointmarkerlist != (int *) NULL) {
delete [] pointmarkerlist;
}
if (pointparamlist != (pointparam *) NULL) {
delete [] pointparamlist;
}
if (tetrahedronlist != (int *) NULL) {
delete [] tetrahedronlist;
}
if (tetrahedronattributelist != (REAL *) NULL) {
delete [] tetrahedronattributelist;
}
if (tetrahedronvolumelist != (REAL *) NULL) {
delete [] tetrahedronvolumelist;
}
if (neighborlist != (int *) NULL) {
delete [] neighborlist;
}
if (trifacelist != (int *) NULL) {
delete [] trifacelist;
}
if (trifacemarkerlist != (int *) NULL) {
delete [] trifacemarkerlist;
}
if (o2facelist != (int *) NULL) {
delete [] o2facelist;
}
if (adjtetlist != (int *) NULL) {
delete [] adjtetlist;
}
if (edgelist != (int *) NULL) {
delete [] edgelist;
}
if (edgemarkerlist != (int *) NULL) {
delete [] edgemarkerlist;
}
if (o2edgelist != (int *) NULL) {
delete [] o2edgelist;
}
if (edgeadjtetlist != (int *) NULL) {
delete [] edgeadjtetlist;
}
if (facetlist != (facet *) NULL) {
facet *f;
polygon *p;
for (i = 0; i < numberoffacets; i++) {
f = &facetlist[i];
for (j = 0; j < f->numberofpolygons; j++) {
p = &f->polygonlist[j];
delete [] p->vertexlist;
}
delete [] f->polygonlist;
if (f->holelist != (REAL *) NULL) {
delete [] f->holelist;
}
}
delete [] facetlist;
}
if (facetmarkerlist != (int *) NULL) {
delete [] facetmarkerlist;
}
if (holelist != (REAL *) NULL) {
delete [] holelist;
}
if (regionlist != (REAL *) NULL) {
delete [] regionlist;
}
if (facetconstraintlist != (REAL *) NULL) {
delete [] facetconstraintlist;
}
if (segmentconstraintlist != (REAL *) NULL) {
delete [] segmentconstraintlist;
}
if (vpointlist != (REAL *) NULL) {
delete [] vpointlist;
}
if (vedgelist != (voroedge *) NULL) {
delete [] vedgelist;
}
if (vfacetlist != (vorofacet *) NULL) {
for (i = 0; i < numberofvfacets; i++) {
delete [] vfacetlist[i].elist;
}
delete [] vfacetlist;
}
if (vcelllist != (int **) NULL) {
for (i = 0; i < numberofvcells; i++) {
delete [] vcelllist[i];
}
delete [] vcelllist;
}
}
// Constructor & destructor.
tetgenio() {initialize();}
~tetgenio() {deinitialize();}
}; // class tetgenio
///////////////////////////////////////////////////////////////////////////////
// //
// 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 conforming; // '-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 reversetetori; // '-R', 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_param; // '-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 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_ang_tol; // '-p', 179.9.
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 minedgelength; // 0.0.
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;
conforming = 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;
reversetetori = 0;
docheck = 0;
quiet = 0;
verbose = 0;
vertexperblock = 4092;
tetrahedraperblock = 8188;
shellfaceperblock = 4092;
nobisect_param = 2;
addsteiner_algo = 1;
coarsen_param = 0;
weighted_param = 0;
fliplinklevel = -1; // No limit on linklevel.
flipstarsize = -1; // No limit on flip star size.
fliplinklevelinc = 1;
reflevel = 3;
optscheme = 7; // 1 & 2 & 4, // min_max_dihedral.
optlevel = 2;
delmaxfliplevel = 1;
order = 1;
steinerleft = -1;
no_sort = 0;
hilbert_order = 52; //-1;
hilbert_limit = 8;
brio_threshold = 64;
brio_ratio = 0.125;
facet_ang_tol = 179.9;
maxvolume = -1.0;
minratio = 2.0;
mindihedral = 0.0; // 5.0;
optmaxdihedral = 165.00; // without -q, default is 179.0
optminsmtdihed = 179.00; // without -q, default is 179.999
optminslidihed = 179.00; // without -q, default is 179.999
epsilon = 1.0e-8;
minedgelength = 0.0;
coarsen_percent = 1.0;
object = NODES;
commandline[0] = '\0';
infilename[0] = '\0';
outfilename[0] = '\0';
addinfilename[0] = '\0';
bgmeshfilename[0] = '\0';
}
}; // class tetgenbehavior
///////////////////////////////////////////////////////////////////////////////
// //
// Robust Geometric predicates //
// //
// Geometric predicates are simple tests of spatial relations of a set of d- //
// dimensional points, such as the orientation test and the point-in-sphere //
// test. Each of these tests is performed by evaluating the sign of a deter- //
// minant of a matrix whose entries are the coordinates of these points. If //
// the computation is performed by using the floating-point numbers, e.g., //
// the single or double precision numbers in C/C++, roundoff error may cause //
// an incorrect result. This may either lead to a wrong result or eventually //
// lead to a failure of the program. Computing the predicates exactly will //
// avoid the error and make the program robust. //
// //
// The following routines are the robust geometric predicates for 3D orient- //
// ation test and point-in-sphere test. They were implemented by Shewchuk. //
// The source code are generously provided by him in the public domain, //
// http://www.cs.cmu.edu/~quake/robust.html. predicates.cxx is a C++ version //
// of the original C code. //
// //
// The original predicates of Shewchuk only use "dynamic filters", i.e., it //
// computes the error at run time step by step. TetGen first adds a "static //
// filter" in each predicate. It estimates the maximal possible error in all //
// cases. So it can safely and quickly answer many easy cases. //
// //
///////////////////////////////////////////////////////////////////////////////
void exactinit(int, int, int, REAL, REAL, REAL);
REAL orient3d(REAL *pa, REAL *pb, REAL *pc, REAL *pd);
REAL insphere(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL *pe);
REAL orient4d(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL *pe,
REAL ah, REAL bh, REAL ch, REAL dh, REAL eh);
///////////////////////////////////////////////////////////////////////////////
// //
// 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 max_min_aspectratio; // Maximize the minimum 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;
max_min_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,
COLLISIONFACE, ACROSSSEG, ACROSSSUB};
// 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 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 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);
// 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);
REAL facedihedral(REAL* pa, REAL* pb, REAL* pc1, REAL* pc2);
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);
///////////////////////////////////////////////////////////////////////////////
// //
// 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);
// 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);
void scarveholes(int, REAL*);
void triangulate(int, arraypool*, arraypool*, int, REAL*);
void unifysubfaces(face*, face*);
void unifysegments();
void mergefacets();
void identifypscedges(point*);
void meshsurface();
void interecursive(shellface** subfacearray, int arraysize, int axis,
REAL, REAL, REAL, REAL, REAL, REAL, int* internum);
void detectinterfaces();
///////////////////////////////////////////////////////////////////////////////
// //
// Constrained Delaunay tetrahedralization //
// //
// A constrained Delaunay tetrahedralization (CDT) is a variation of a Dela- //
// unay tetrahedralization (DT) that is constrained to respect the boundary //
// of a 3D PLC (domain). In a CDT of a 3D PLC, every vertex or edge of the //
// PLC is also a vertex or an edge of the CDT, every polygon of the PLC is a //
// union of triangles of the CDT. A crucial difference between a CDT and a //
// DT is that triangles in the PLC's polygons are not required to be locally //
// Delaunay, which frees the CDT to better respect the PLC's polygons. CDTs //
// have optimal properties similar to those of DTs. //
// //
// Steiner Points and Steiner CDTs. It is known that even a simple 3D polyh- //
// edron may not have a tetrahedralization which only uses its own vertices. //
// Some extra points, so-called "Steiner points" are needed in order to form //
// a tetrahedralization of such polyhedron. It is true for tetrahedralizing //
// a 3D PLC as well. A Steiner CDT of a 3D PLC is a CDT containing Steiner //
// points. The CDT algorithms of TetGen in general create Steiner CDTs. //
// Almost all of the Steiner points are added in the edges of the PLC. They //
// guarantee the existence of a CDT of the modified PLC. //
// //
// The routine constraineddelaunay() starts from a DT of the vertices of a //
// PLC and creates a (Steiner) CDT of the PLC (including Steiner points). It //
// is constructed by two steps, (1) segment recovery and (2) facet (polygon) //
// recovery. Each step is accomplished by its own algorithm. //
// //
// The routine delaunizesegments() implements the segment recovery algorithm //
// of Si, H. and Gaertner, K. "Meshing Piecewise Linear Complexes by Constr- //
// ained Delaunay Tetrahedralizations," In Proceedings of the 14th Internat- //
// ional Meshing Roundtable, 147--163, 2005. It adds Steiner points into //
// non-Delaunay segments until all subsegments appear together in a DT. The //
// running time of this algorithm is proportional to the number of added //
// Steiner points. //
// //
// There are two incremental facet recovery algorithms: the cavity re-trian- //
// gulation algorithm of Si, H. and Gaertner, K. "3D Boundary Recovery by //
// Constrained Delaunay Tetrahedralization," International Journal for Numer-//
// ical Methods in Engineering, 85:1341-1364, 2011, and the flip algorithm //
// of Shewchuk, J. "Updating and Constructing Constrained Delaunay and //
// Constrained Regular Triangulations by Flips." In Proceedings of the 19th //
// ACM Symposium on Computational Geometry, 86-95, 2003. //
// //
// It is guaranteed in theory, no Steiner point is needed in both algorithms //
// However, a facet with non-coplanar vertices might cause the additions of //
// Steiner points. It is discussed in the paper of Si, H., and Shewchuk, J.,//
// "Incrementally Constructing and Updating Constrained Delaunay //
// Tetrahedralizations with Finite Precision Coordinates." In Proceedings of //
// the 21th International Meshing Roundtable, 2012. //
// //
// Our implementation of the facet recovery algorithms recover a "missing //
// region" at a time. Each missing region is a subset of connected interiors //
// of a polygon. The routine formcavity() creates the cavity of crossing //
// tetrahedra of the missing region. //
// //
// The cavity re-triangulation algorithm is implemented by three subroutines,//
// delaunizecavity(), fillcavity(), and carvecavity(). Since it may fail due //
// to non-coplanar vertices, the subroutine restorecavity() is used to rest- //
// ore the original cavity. //
// //
// The routine flipinsertfacet() implements the flip algorithm. The subrout- //
// ine flipcertify() is used to maintain the priority queue of flips. //
// //
// The routine refineregion() is called when the facet recovery algorithm //
// fail to recover a missing region. It inserts Steiner points to refine the //
// missing region. In order to avoid inserting Steiner points very close to //
// existing segments. The classical encroachment rules of the Delaunay //
// refinement algorithm are used to choose the Steiner points. //
// //
// The routine constrainedfacets() does the facet recovery by using either //
// the cavity re-triangulation algorithm (default) or the flip algorithm. It //
// results a CDT of the (modified) PLC (including Steiner points). //
// //
///////////////////////////////////////////////////////////////////////////////
void makesegmentendpointsmap();
enum interresult finddirection(triface* searchtet, point endpt);
enum interresult scoutsegment(point, point, triface*, point*, arraypool*);
int getsteinerptonsegment(face* seg, point refpt, point steinpt);
void delaunizesegments();
enum interresult scoutsubface(face* searchsh, triface* searchtet);
void formregion(face*, arraypool*, arraypool*, arraypool*);
int scoutcrossedge(triface& crosstet, arraypool*, arraypool*);
bool formcavity(triface*, arraypool*, arraypool*, arraypool*, arraypool*,
arraypool*, arraypool*);
// Facet recovery by cavity re-triangulation [Si and Gaertner 2011].
void delaunizecavity(arraypool*, arraypool*, arraypool*, arraypool*,
arraypool*, arraypool*);
bool fillcavity(arraypool*, arraypool*, arraypool*, arraypool*,
arraypool*, arraypool*, triface* crossedge);
void carvecavity(arraypool*, arraypool*, arraypool*);
void restorecavity(arraypool*, arraypool*, arraypool*, arraypool*);
// Facet recovery by flips [Shewchuk 2003].
void flipcertify(triface *chkface, badface **pqueue, point, point, point);
void flipinsertfacet(arraypool*, arraypool*, arraypool*, arraypool*);
bool fillregion(arraypool* missingshs, arraypool*, arraypool* newshs);
int insertpoint_cdt(point, triface*, face*, face*, insertvertexflags*,
arraypool*, arraypool*, arraypool*, arraypool*,
arraypool*, arraypool*);
void refineregion(face&, arraypool*, arraypool*, arraypool*, arraypool*,
arraypool*, arraypool*);
void constrainedfacets();
void constraineddelaunay(clock_t&);
///////////////////////////////////////////////////////////////////////////////
// //
// 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, 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();
void reconstructmesh();
int scoutpoint(point, triface*, int randflag);
REAL getpointmeshsize(point, triface*, int iloc);
void interpolatemeshsize();
void insertconstrainedpoints(point *insertarray, int arylen, int rejflag);
void insertconstrainedpoints(tetgenio *addio);
void collectremovepoints(arraypool *remptlist);
void meshcoarsening();
///////////////////////////////////////////////////////////////////////////////
// //
// Mesh refinement //
// //
// The purpose of mesh refinement is to obtain a tetrahedral mesh with well- //
// -shaped tetrahedra and appropriate mesh size. It is necessary to insert //
// new Steiner points to achieve this property. The questions are (1) how to //
// choose the Steiner points? and (2) how to insert them? //
// //
// Delaunay refinement is a technique first developed by Chew [1989] and //
// Ruppert [1993, 1995] to generate quality triangular meshes in the plane. //
// It provides guarantee on the smallest angle of the triangles. Rupper's //
// algorithm guarantees that the mesh is size-optimal (to within a constant //
// factor) among all meshes with the same quality. //
// Shewchuk generalized Ruppert's algorithm into 3D in his PhD thesis //
// [Shewchuk 1997]. A short version of his algorithm appears in "Tetrahedral //
// Mesh Generation by Delaunay Refinement," In Proceedings of the 14th ACM //
// Symposium on Computational Geometry, 86-95, 1998. It guarantees that all //
// tetrahedra of the output mesh have a "radius-edge ratio" (equivalent to //
// the minimal face angle) bounded. However, it does not remove slivers, a //
// type of very flat tetrahedra which can have no small face angles but have //
// very small (and large) dihedral angles. Moreover, it may not terminate if //
// the input PLC contains "sharp features", e.g., two edges (or two facets) //
// meet at an acute angle (or dihedral angle). //
// //
// TetGen uses the basic Delaunay refinement scheme to insert Steiner points.//
// While it always maintains a constrained Delaunay mesh. The algorithm is //
// described in Si, H., "Adaptive Constrained Delaunay Mesh Generation," //
// International Journal for Numerical Methods in Engineering, 75:856-880. //
// This algorithm always terminates and sharp features are easily preserved. //
// The mesh has good quality (same as Shewchuk's Delaunay refinement algori- //
// thm) in the bulk of the mesh domain. Moreover, it supports the generation //
// of adaptive mesh according to a (isotropic) mesh sizing function. //
// //
///////////////////////////////////////////////////////////////////////////////
void makefacetverticesmap();
int segsegadjacent(face *, face *);
int segfacetadjacent(face *checkseg, face *checksh);
int facetfacetadjacent(face *, face *);
int checkseg4encroach(point pa, point pb, point checkpt);
int checkseg4split(face *chkseg, point&, int&);
int splitsegment(face *splitseg, point encpt, REAL, point, point, int, int);
void repairencsegs(int chkencflag);
void enqueuesubface(memorypool*, face*);
int checkfac4encroach(point, point, point, point checkpt, REAL*, REAL*);
int checkfac4split(face *chkfac, point& encpt, int& qflag, REAL *ccent);
int splitsubface(face *splitfac, point, point, int qflag, REAL *ccent, int);
void repairencfacs(int chkencflag);
void enqueuetetrahedron(triface*);
int checktet4split(triface *chktet, int& qflag, REAL *ccent);
int splittetrahedron(triface* splittet,int qflag,REAL *ccent, int);
void repairbadtets(int chkencflag);
void delaunayrefinement();
///////////////////////////////////////////////////////////////////////////////
// //
// 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();
///////////////////////////////////////////////////////////////////////////////
// //
// Mesh check and statistics //
// //
///////////////////////////////////////////////////////////////////////////////
// Mesh validations.
int checkmesh(int topoflag);
int checkshells();
int checksegments();
int checkdelaunay();
int checkregular(int);
int checkconforming(int);
// Mesh statistics.
void printfcomma(unsigned long n);
void qualitystatistics();
void memorystatistics();
void statistics();
///////////////////////////////////////////////////////////////////////////////
// //
// Mesh output //
// //
///////////////////////////////////////////////////////////////////////////////
void jettisonnodes();
void highorder();
void numberedges();
void outnodes(tetgenio*);
void outmetrics(tetgenio*);
void outelements(tetgenio*);
void outfaces(tetgenio*);
void outhullfaces(tetgenio*);
void outsubfaces(tetgenio*);
void outedges(tetgenio*);
void outsubsegments(tetgenio*);
void outneighbors(tetgenio*);
void outvoronoi(tetgenio*);
void outsmesh(char*);
void outmesh2medit(char*);
void outmesh2vtk(char*);
///////////////////////////////////////////////////////////////////////////////
// //
// Constructor & destructor //
// //
///////////////////////////////////////////////////////////////////////////////
tetgenmesh()
{
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;
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 = 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 (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;
}
}
~tetgenmesh()
{
freememory();
} // ~tetgenmesh()
}; // End of class tetgenmesh.
///////////////////////////////////////////////////////////////////////////////
// //
// tetrahedralize() Interface for using TetGen's library to generate //
// Delaunay tetrahedralizations, constrained Delaunay //
// tetrahedralizations, quality tetrahedral meshes. //
// //
// 'in' is an object of 'tetgenio' which contains a PLC you want to tetrahed-//
// ralize or a previously generated tetrahedral mesh you want to refine. It //
// must not be a NULL. 'out' is another object of 'tetgenio' for storing the //
// generated tetrahedral mesh. It can be a NULL. If so, the output will be //
// saved to file(s). If 'bgmin' != NULL, it contains a background mesh which //
// defines a mesh size function. //
// //
///////////////////////////////////////////////////////////////////////////////
void tetrahedralize(tetgenbehavior *b, tetgenio *in, tetgenio *out,
tetgenio *addin = NULL, tetgenio *bgmin = NULL);
#ifdef TETLIBRARY
void tetrahedralize(char *switches, tetgenio *in, tetgenio *out,
tetgenio *addin = NULL, tetgenio *bgmin = NULL);
#endif // #ifdef TETLIBRARY
///////////////////////////////////////////////////////////////////////////////
// //
// terminatetetgen() Terminate TetGen with a given exit code. //
// //
///////////////////////////////////////////////////////////////////////////////
inline void terminatetetgen(tetgenmesh *m, int x)
{
// Release the allocated memory.
if (m) {
m->freememory();
}
#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");
printf("Hint: use -T option to set a smaller tolerance.\n");
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[pointmtrindex + sizeoftensor - 1] = value;
}
inline REAL tetgenmesh::getpointinsradius(point pt)
{
return pt[pointmtrindex + sizeoftensor - 1];
}
// 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 {
assert((point) searchtet.tet[7] == pa); // SELF_CHECK
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 {
assert((point) searchsh.sh[5] == pa); // SELF_CHECK
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 // #ifndef tetgenH
View Git Blame Copy Permalink