1///////////////////////////////////////////////////////////////////////////////
2// //
3// TetGen //
4// //
5// A Quality Tetrahedral Mesh Generator and A 3D Delaunay Triangulator //
6// //
7// Version 1.5 //
8// November 4, 2013 //
9// //
10// TetGen is freely available through the website: http://www.tetgen.org. //
11// It may be copied, modified, and redistributed for non-commercial use. //
12// Please consult the file LICENSE for the detailed copyright notices. //
13// //
14///////////////////////////////////////////////////////////////////////////////
15
16
17#ifndef tetgenH
18#define tetgenH
19
20#define _CRT_SECURE_NO_WARNINGS
21
22// To compile TetGen as a library instead of an executable program, define
23// the TETLIBRARY symbol.
24
25#define TETLIBRARY
26
27// Uncomment the following line to disable assert macros. These macros were
28// inserted in the code where I hoped to catch bugs. They may slow down the
29// speed of TetGen.
30
31// #define NDEBUG
32
33// TetGen default uses the double precision (64 bit) for a real number.
34// Alternatively, one can use the single precision (32 bit) 'float' if the
35// memory is limited.
36
37#define REAL double // #define REAL float
38
39// Maximum number of characters in a file name (including the null).
40
41#define FILENAMESIZE 1024
42
43// Maximum number of chars in a line read from a file (including the null).
44
45#define INPUTLINESIZE 2048
46
47// TetGen only uses the C standard library.
48
49#include <stdio.h>
50#include <stdlib.h>
51#include <string.h>
52#include <math.h>
53#include <time.h>
54#include <assert.h>
55
56// The types 'intptr_t' and 'uintptr_t' are signed and unsigned integer types,
57// respectively. They are guaranteed to be the same width as a pointer.
58// They are defined in <stdint.h> by the C99 Standard. However, Microsoft
59// Visual C++ 2003 -- 2008 (Visual C++ 7.1 - 9) doesn't ship with this header
60// file. In such case, we can define them by ourself.
61// Update (learned from Stack Overflow): Visual Studio 2010 and Visual C++ 2010
62// Express both have stdint.h
63
64// The following piece of code was provided by Steven Johnson (MIT). Define the
65// symbol _MSC_VER if you are using Microsoft Visual C++. Moreover, define
66// the _WIN64 symbol if you are running TetGen on Win64 systems.
67
68#ifdef _MSC_VER // Microsoft Visual C++
69# ifdef _WIN64
70 typedef __int64 intptr_t;
71 typedef unsigned __int64 uintptr_t;
72# else // not _WIN64
73 typedef int intptr_t;
74 typedef unsigned int uintptr_t;
75# endif
76#else // not Visual C++
77# include <stdint.h>
78#endif
79
80///////////////////////////////////////////////////////////////////////////////
81// //
82// tetgenio //
83// //
84// A structure for transferring data into and out of TetGen's mesh structure,//
85// 'tetgenmesh' (declared below). //
86// //
87// The input of TetGen is either a 3D point set, or a 3D piecewise linear //
88// complex (PLC), or a tetrahedral mesh. Depending on the input object and //
89// the specified options, the output of TetGen is either a Delaunay (or wei- //
90// ghted Delaunay) tetrahedralization, or a constrained (Delaunay) tetrahed- //
91// ralization, or a quality tetrahedral mesh. //
92// //
93// A piecewise linear complex (PLC) represents a 3D polyhedral domain with //
94// possibly internal boundaries(subdomains). It is introduced in [Miller et //
95// al, 1996]. Basically it is a set of "cells", i.e., vertices, edges, poly- //
96// gons, and polyhedra, and the intersection of any two of its cells is the //
97// union of other cells of it. //
98// //
99// TetGen uses a set of files to describe the inputs and outputs. Each file //
100// is identified from its file extension (.node, .ele, .face, .edge, etc). //
101// //
102// The 'tetgenio' structure is a collection of arrays of data, i.e., points, //
103// facets, tetrahedra, and so forth. It contains functions to read and write //
104// (input and output) files of TetGen as well as other supported mesh files. //
105// //
106// Once an object of tetgenio is declared, no array is created. One has to //
107// allocate enough memory for them. On deletion of this object, the memory //
108// occupied by these arrays needs to be freed. The routine deinitialize() //
109// will be automatically called. It frees the memory for an array if it is //
110// not a NULL. Note that it assumes that the memory is allocated by the C++ //
111// "new" operator. Otherwise, the user is responsible to free them and all //
112// pointers must be NULL before the call of the destructor. //
113// //
114///////////////////////////////////////////////////////////////////////////////
115
116class tetgenio {
117
118public:
119
120 // A "polygon" describes a simple polygon (no holes). It is not necessarily
121 // convex. Each polygon contains a number of corners (points) and the same
122 // number of sides (edges). The points of the polygon must be given in
123 // either counterclockwise or clockwise order and they form a ring, so
124 // every two consecutive points forms an edge of the polygon.
125 typedef struct {
126 int *vertexlist;
127 int numberofvertices;
128 } polygon;
129
130 // A "facet" describes a polygonal region possibly with holes, edges, and
131 // points floating in it. Each facet consists of a list of polygons and
132 // a list of hole points (which lie strictly inside holes).
133 typedef struct {
134 polygon *polygonlist;
135 int numberofpolygons;
136 REAL *holelist;
137 int numberofholes;
138 } facet;
139
140 // A "voroedge" is an edge of the Voronoi diagram. It corresponds to a
141 // Delaunay face. Each voroedge is either a line segment connecting
142 // two Voronoi vertices or a ray starting from a Voronoi vertex to an
143 // "infinite vertex". 'v1' and 'v2' are two indices pointing to the
144 // list of Voronoi vertices. 'v1' must be non-negative, while 'v2' may
145 // be -1 if it is a ray, in this case, the unit normal of this ray is
146 // given in 'vnormal'.
147 typedef struct {
148 int v1, v2;
149 REAL vnormal[3];
150 } voroedge;
151
152 // A "vorofacet" is an facet of the Voronoi diagram. It corresponds to a
153 // Delaunay edge. Each Voronoi facet is a convex polygon formed by a
154 // list of Voronoi edges, it may not be closed. 'c1' and 'c2' are two
155 // indices pointing into the list of Voronoi cells, i.e., the two cells
156 // share this facet. 'elist' is an array of indices pointing into the
157 // list of Voronoi edges, 'elist[0]' saves the number of Voronoi edges
158 // (including rays) of this facet.
159 typedef struct {
160 int c1, c2;
161 int *elist;
162 } vorofacet;
163
164
165 // Additional parameters associated with an input (or mesh) vertex.
166 // These informations are provided by CAD libraries.
167 typedef struct {
168 REAL uv[2];
169 int tag;
170 int type; // 0, 1, or 2.
171 } pointparam;
172
173 // Callback functions for meshing PSCs.
174 typedef REAL (* GetVertexParamOnEdge)(void*, int, int);
175 typedef void (* GetSteinerOnEdge)(void*, int, REAL, REAL*);
176 typedef void (* GetVertexParamOnFace)(void*, int, int, REAL*);
177 typedef void (* GetEdgeSteinerParamOnFace)(void*, int, REAL, int, REAL*);
178 typedef void (* GetSteinerOnFace)(void*, int, REAL*, REAL*);
179
180 // A callback function for mesh refinement.
181 typedef bool (* TetSizeFunc)(REAL*, REAL*, REAL*, REAL*, REAL*, REAL);
182
183 // Items are numbered starting from 'firstnumber' (0 or 1), default is 0.
184 int firstnumber;
185
186 // Dimension of the mesh (2 or 3), default is 3.
187 int mesh_dim;
188
189 // Does the lines in .node file contain index or not, default is 1.
190 int useindex;
191
192 // 'pointlist': An array of point coordinates. The first point's x
193 // coordinate is at index [0] and its y coordinate at index [1], its
194 // z coordinate is at index [2], followed by the coordinates of the
195 // remaining points. Each point occupies three REALs.
196 // 'pointattributelist': An array of point attributes. Each point's
197 // attributes occupy 'numberofpointattributes' REALs.
198 // 'pointmtrlist': An array of metric tensors at points. Each point's
199 // tensor occupies 'numberofpointmtr' REALs.
200 // 'pointmarkerlist': An array of point markers; one integer per point.
201 REAL *pointlist;
202 REAL *pointattributelist;
203 REAL *pointmtrlist;
204 int *pointmarkerlist;
205 pointparam *pointparamlist;
206 int numberofpoints;
207 int numberofpointattributes;
208 int numberofpointmtrs;
209
210 // 'tetrahedronlist': An array of tetrahedron corners. The first
211 // tetrahedron's first corner is at index [0], followed by its other
212 // corners, followed by six nodes on the edges of the tetrahedron if the
213 // second order option (-o2) is applied. Each tetrahedron occupies
214 // 'numberofcorners' ints. The second order nodes are ouput only.
215 // 'tetrahedronattributelist': An array of tetrahedron attributes. Each
216 // tetrahedron's attributes occupy 'numberoftetrahedronattributes' REALs.
217 // 'tetrahedronvolumelist': An array of constraints, i.e. tetrahedron's
218 // volume; one REAL per element. Input only.
219 // 'neighborlist': An array of tetrahedron neighbors; 4 ints per element.
220 // Output only.
221 int *tetrahedronlist;
222 REAL *tetrahedronattributelist;
223 REAL *tetrahedronvolumelist;
224 int *neighborlist;
225 int numberoftetrahedra;
226 int numberofcorners;
227 int numberoftetrahedronattributes;
228
229 // 'facetlist': An array of facets. Each entry is a structure of facet.
230 // 'facetmarkerlist': An array of facet markers; one int per facet.
231 facet *facetlist;
232 int *facetmarkerlist;
233 int numberoffacets;
234
235 // 'holelist': An array of holes (in volume). Each hole is given by a
236 // seed (point) which lies strictly inside it. The first seed's x, y and z
237 // coordinates are at indices [0], [1] and [2], followed by the
238 // remaining seeds. Three REALs per hole.
239 REAL *holelist;
240 int numberofholes;
241
242 // 'regionlist': An array of regions (subdomains). Each region is given by
243 // a seed (point) which lies strictly inside it. The first seed's x, y and
244 // z coordinates are at indices [0], [1] and [2], followed by the regional
245 // attribute at index [3], followed by the maximum volume at index [4].
246 // Five REALs per region.
247 // Note that each regional attribute is used only if you select the 'A'
248 // switch, and each volume constraint is used only if you select the
249 // 'a' switch (with no number following).
250 REAL *regionlist;
251 int numberofregions;
252
253 // 'facetconstraintlist': An array of facet constraints. Each constraint
254 // specifies a maximum area bound on the subfaces of that facet. The
255 // first facet constraint is given by a facet marker at index [0] and its
256 // maximum area bound at index [1], followed by the remaining facet con-
257 // straints. Two REALs per facet constraint. Note: the facet marker is
258 // actually an integer.
259 REAL *facetconstraintlist;
260 int numberoffacetconstraints;
261
262 // 'segmentconstraintlist': An array of segment constraints. Each constraint
263 // specifies a maximum length bound on the subsegments of that segment.
264 // The first constraint is given by the two endpoints of the segment at
265 // index [0] and [1], and the maximum length bound at index [2], followed
266 // by the remaining segment constraints. Three REALs per constraint.
267 // Note the segment endpoints are actually integers.
268 REAL *segmentconstraintlist;
269 int numberofsegmentconstraints;
270
271
272 // 'trifacelist': An array of face (triangle) corners. The first face's
273 // three corners are at indices [0], [1] and [2], followed by the remaining
274 // faces. Three ints per face.
275 // 'trifacemarkerlist': An array of face markers; one int per face.
276 // 'o2facelist': An array of second order nodes (on the edges) of the face.
277 // It is output only if the second order option (-o2) is applied. The
278 // first face's three second order nodes are at [0], [1], and [2],
279 // followed by the remaining faces. Three ints per face.
280 // 'adjtetlist': An array of adjacent tetrahedra to the faces. The first
281 // face's two adjacent tetrahedra are at indices [0] and [1], followed by
282 // the remaining faces. A '-1' indicates outside (no adj. tet). This list
283 // is output when '-nn' switch is used. Output only.
284 int *trifacelist;
285 int *trifacemarkerlist;
286 int *o2facelist;
287 int *adjtetlist;
288 int numberoftrifaces;
289
290 // 'edgelist': An array of edge endpoints. The first edge's endpoints
291 // are at indices [0] and [1], followed by the remaining edges.
292 // Two ints per edge.
293 // 'edgemarkerlist': An array of edge markers; one int per edge.
294 // 'o2edgelist': An array of midpoints of edges. It is output only if the
295 // second order option (-o2) is applied. One int per edge.
296 // 'edgeadjtetlist': An array of adjacent tetrahedra to the edges. One
297 // tetrahedron (an integer) per edge.
298 int *edgelist;
299 int *edgemarkerlist;
300 int *o2edgelist;
301 int *edgeadjtetlist;
302 int numberofedges;
303
304 // 'vpointlist': An array of Voronoi vertex coordinates (like pointlist).
305 // 'vedgelist': An array of Voronoi edges. Each entry is a 'voroedge'.
306 // 'vfacetlist': An array of Voronoi facets. Each entry is a 'vorofacet'.
307 // 'vcelllist': An array of Voronoi cells. Each entry is an array of
308 // indices pointing into 'vfacetlist'. The 0th entry is used to store
309 // the length of this array.
310 REAL *vpointlist;
311 voroedge *vedgelist;
312 vorofacet *vfacetlist;
313 int **vcelllist;
314 int numberofvpoints;
315 int numberofvedges;
316 int numberofvfacets;
317 int numberofvcells;
318
319 // Variable (and callback functions) for meshing PSCs.
320 void *geomhandle;
321 GetVertexParamOnEdge getvertexparamonedge;
322 GetSteinerOnEdge getsteineronedge;
323 GetVertexParamOnFace getvertexparamonface;
324 GetEdgeSteinerParamOnFace getedgesteinerparamonface;
325 GetSteinerOnFace getsteineronface;
326
327 // A callback function.
328 TetSizeFunc tetunsuitable;
329
330 // Input & output routines.
331 bool load_node_call(FILE* infile, int markers, int uvflag, char*);
332 bool load_node(char*);
333 bool load_edge(char*);
334 bool load_face(char*);
335 bool load_tet(char*);
336 bool load_vol(char*);
337 bool load_var(char*);
338 bool load_mtr(char*);
339 bool load_pbc(char*);
340 bool load_poly(char*);
341 bool load_off(char*);
342 bool load_ply(char*);
343 bool load_stl(char*);
344 bool load_vtk(char*);
345 bool load_medit(char*, int);
346 bool load_plc(char*, int);
347 bool load_tetmesh(char*, int);
348 void save_nodes(char*);
349 void save_elements(char*);
350 void save_faces(char*);
351 void save_edges(char*);
352 void save_neighbors(char*);
353 void save_poly(char*);
354 void save_faces2smesh(char*);
355
356 // Read line and parse string functions.
357 char *readline(char* string, FILE* infile, int *linenumber);
358 char *findnextfield(char* string);
359 char *readnumberline(char* string, FILE* infile, char* infilename);
360 char *findnextnumber(char* string);
361
362 static void init(polygon* p) {
363 p->vertexlist = (int *) NULL;
364 p->numberofvertices = 0;
365 }
366
367 static void init(facet* f) {
368 f->polygonlist = (polygon *) NULL;
369 f->numberofpolygons = 0;
370 f->holelist = (REAL *) NULL;
371 f->numberofholes = 0;
372 }
373
374 // Initialize routine.
375 void initialize()
376 {
377 firstnumber = 0;
378 mesh_dim = 3;
379 useindex = 1;
380
381 pointlist = (REAL *) NULL;
382 pointattributelist = (REAL *) NULL;
383 pointmtrlist = (REAL *) NULL;
384 pointmarkerlist = (int *) NULL;
385 pointparamlist = (pointparam *) NULL;
386 numberofpoints = 0;
387 numberofpointattributes = 0;
388 numberofpointmtrs = 0;
389
390 tetrahedronlist = (int *) NULL;
391 tetrahedronattributelist = (REAL *) NULL;
392 tetrahedronvolumelist = (REAL *) NULL;
393 neighborlist = (int *) NULL;
394 numberoftetrahedra = 0;
395 numberofcorners = 4;
396 numberoftetrahedronattributes = 0;
397
398 trifacelist = (int *) NULL;
399 trifacemarkerlist = (int *) NULL;
400 o2facelist = (int *) NULL;
401 adjtetlist = (int *) NULL;
402 numberoftrifaces = 0;
403
404 edgelist = (int *) NULL;
405 edgemarkerlist = (int *) NULL;
406 o2edgelist = (int *) NULL;
407 edgeadjtetlist = (int *) NULL;
408 numberofedges = 0;
409
410 facetlist = (facet *) NULL;
411 facetmarkerlist = (int *) NULL;
412 numberoffacets = 0;
413
414 holelist = (REAL *) NULL;
415 numberofholes = 0;
416
417 regionlist = (REAL *) NULL;
418 numberofregions = 0;
419
420 facetconstraintlist = (REAL *) NULL;
421 numberoffacetconstraints = 0;
422 segmentconstraintlist = (REAL *) NULL;
423 numberofsegmentconstraints = 0;
424
425
426 vpointlist = (REAL *) NULL;
427 vedgelist = (voroedge *) NULL;
428 vfacetlist = (vorofacet *) NULL;
429 vcelllist = (int **) NULL;
430 numberofvpoints = 0;
431 numberofvedges = 0;
432 numberofvfacets = 0;
433 numberofvcells = 0;
434
435 tetunsuitable = NULL;
436
437 geomhandle = NULL;
438 getvertexparamonedge = NULL;
439 getsteineronedge = NULL;
440 getvertexparamonface = NULL;
441 getedgesteinerparamonface = NULL;
442 getsteineronface = NULL;
443 }
444
445 // Free the memory allocated in 'tetgenio'. Note that it assumes that the
446 // memory was allocated by the "new" operator (C++).
447 void deinitialize()
448 {
449 int i, j;
450
451 if (pointlist != (REAL *) NULL) {
452 delete [] pointlist;
453 }
454 if (pointattributelist != (REAL *) NULL) {
455 delete [] pointattributelist;
456 }
457 if (pointmtrlist != (REAL *) NULL) {
458 delete [] pointmtrlist;
459 }
460 if (pointmarkerlist != (int *) NULL) {
461 delete [] pointmarkerlist;
462 }
463 if (pointparamlist != (pointparam *) NULL) {
464 delete [] pointparamlist;
465 }
466
467 if (tetrahedronlist != (int *) NULL) {
468 delete [] tetrahedronlist;
469 }
470 if (tetrahedronattributelist != (REAL *) NULL) {
471 delete [] tetrahedronattributelist;
472 }
473 if (tetrahedronvolumelist != (REAL *) NULL) {
474 delete [] tetrahedronvolumelist;
475 }
476 if (neighborlist != (int *) NULL) {
477 delete [] neighborlist;
478 }
479
480 if (trifacelist != (int *) NULL) {
481 delete [] trifacelist;
482 }
483 if (trifacemarkerlist != (int *) NULL) {
484 delete [] trifacemarkerlist;
485 }
486 if (o2facelist != (int *) NULL) {
487 delete [] o2facelist;
488 }
489 if (adjtetlist != (int *) NULL) {
490 delete [] adjtetlist;
491 }
492
493 if (edgelist != (int *) NULL) {
494 delete [] edgelist;
495 }
496 if (edgemarkerlist != (int *) NULL) {
497 delete [] edgemarkerlist;
498 }
499 if (o2edgelist != (int *) NULL) {
500 delete [] o2edgelist;
501 }
502 if (edgeadjtetlist != (int *) NULL) {
503 delete [] edgeadjtetlist;
504 }
505
506 if (facetlist != (facet *) NULL) {
507 facet *f;
508 polygon *p;
509 for (i = 0; i < numberoffacets; i++) {
510 f = &facetlist[i];
511 for (j = 0; j < f->numberofpolygons; j++) {
512 p = &f->polygonlist[j];
513 delete [] p->vertexlist;
514 }
515 delete [] f->polygonlist;
516 if (f->holelist != (REAL *) NULL) {
517 delete [] f->holelist;
518 }
519 }
520 delete [] facetlist;
521 }
522 if (facetmarkerlist != (int *) NULL) {
523 delete [] facetmarkerlist;
524 }
525
526 if (holelist != (REAL *) NULL) {
527 delete [] holelist;
528 }
529 if (regionlist != (REAL *) NULL) {
530 delete [] regionlist;
531 }
532 if (facetconstraintlist != (REAL *) NULL) {
533 delete [] facetconstraintlist;
534 }
535 if (segmentconstraintlist != (REAL *) NULL) {
536 delete [] segmentconstraintlist;
537 }
538 if (vpointlist != (REAL *) NULL) {
539 delete [] vpointlist;
540 }
541 if (vedgelist != (voroedge *) NULL) {
542 delete [] vedgelist;
543 }
544 if (vfacetlist != (vorofacet *) NULL) {
545 for (i = 0; i < numberofvfacets; i++) {
546 delete [] vfacetlist[i].elist;
547 }
548 delete [] vfacetlist;
549 }
550 if (vcelllist != (int **) NULL) {
551 for (i = 0; i < numberofvcells; i++) {
552 delete [] vcelllist[i];
553 }
554 delete [] vcelllist;
555 }
556 }
557
558 // Constructor & destructor.
559 tetgenio() {initialize();}
560 ~tetgenio() {deinitialize();}
561
562}; // class tetgenio
563
564///////////////////////////////////////////////////////////////////////////////
565// //
566// tetgenbehavior //
567// //
568// A structure for maintaining the switches and parameters used by TetGen's //
569// mesh data structure and algorithms. //
570// //
571// All switches and parameters are initialized with default values. They can //
572// be set by the command line arguments (a list of strings) of TetGen. //
573// //
574// NOTE: Some of the switches are incompatible. While some may depend on //
575// other switches. The routine parse_commandline() sets the switches from //
576// the command line (a list of strings) and checks the consistency of the //
577// applied switches. //
578// //
579///////////////////////////////////////////////////////////////////////////////
580
581class tetgenbehavior {
582
583public:
584
585 // Switches of TetGen.
586 int plc; // '-p', 0.
587 int psc; // '-s', 0.
588 int refine; // '-r', 0.
589 int quality; // '-q', 0.
590 int nobisect; // '-Y', 0.
591 int coarsen; // '-R', 0.
592 int weighted; // '-w', 0.
593 int brio_hilbert; // '-b', 1.
594 int incrflip; // '-l', 0.
595 int flipinsert; // '-L', 0.
596 int metric; // '-m', 0.
597 int varvolume; // '-a', 0.
598 int fixedvolume; // '-a', 0.
599 int regionattrib; // '-A', 0.
600 int conforming; // '-D', 0.
601 int insertaddpoints; // '-i', 0.
602 int diagnose; // '-d', 0.
603 int convex; // '-c', 0.
604 int nomergefacet; // '-M', 0.
605 int nomergevertex; // '-M', 0.
606 int noexact; // '-X', 0.
607 int nostaticfilter; // '-X', 0.
608 int zeroindex; // '-z', 0.
609 int facesout; // '-f', 0.
610 int edgesout; // '-e', 0.
611 int neighout; // '-n', 0.
612 int voroout; // '-v', 0.
613 int meditview; // '-g', 0.
614 int vtkview; // '-k', 0.
615 int nobound; // '-B', 0.
616 int nonodewritten; // '-N', 0.
617 int noelewritten; // '-E', 0.
618 int nofacewritten; // '-F', 0.
619 int noiterationnum; // '-I', 0.
620 int nojettison; // '-J', 0.
621 int reversetetori; // '-R', 0.
622 int docheck; // '-C', 0.
623 int quiet; // '-Q', 0.
624 int verbose; // '-V', 0.
625
626 // Parameters of TetGen.
627 int vertexperblock; // '-x', 4092.
628 int tetrahedraperblock; // '-x', 8188.
629 int shellfaceperblock; // '-x', 2044.
630 int nobisect_param; // '-Y', 2.
631 int addsteiner_algo; // '-Y/', 1.
632 int coarsen_param; // '-R', 0.
633 int weighted_param; // '-w', 0.
634 int fliplinklevel; // -1.
635 int flipstarsize; // -1.
636 int fliplinklevelinc; // 1.
637 int reflevel; // '-D', 3.
638 int optlevel; // '-O', 2.
639 int optscheme; // '-O', 7.
640 int delmaxfliplevel; // 1.
641 int order; // '-o', 1.
642 int steinerleft; // '-S', 0.
643 int no_sort; // 0.
644 int hilbert_order; // '-b///', 52.
645 int hilbert_limit; // '-b//' 8.
646 int brio_threshold; // '-b' 64.
647 REAL brio_ratio; // '-b/' 0.125.
648 REAL facet_ang_tol; // '-p', 179.9.
649 REAL maxvolume; // '-a', -1.0.
650 REAL minratio; // '-q', 0.0.
651 REAL mindihedral; // '-q', 5.0.
652 REAL optmaxdihedral; // 165.0.
653 REAL optminsmtdihed; // 179.0.
654 REAL optminslidihed; // 179.0.
655 REAL epsilon; // '-T', 1.0e-8.
656 REAL minedgelength; // 0.0.
657 REAL coarsen_percent; // -R1/#, 1.0.
658
659 // Strings of command line arguments and input/output file names.
660 char commandline[1024];
661 char infilename[1024];
662 char outfilename[1024];
663 char addinfilename[1024];
664 char bgmeshfilename[1024];
665
666 // The input object of TetGen. They are recognized by either the input
667 // file extensions or by the specified options.
668 // Currently the following objects are supported:
669 // - NODES, a list of nodes (.node);
670 // - POLY, a piecewise linear complex (.poly or .smesh);
671 // - OFF, a polyhedron (.off, Geomview's file format);
672 // - PLY, a polyhedron (.ply, file format from gatech, only ASCII);
673 // - STL, a surface mesh (.stl, stereolithography format);
674 // - MEDIT, a surface mesh (.mesh, Medit's file format);
675 // - MESH, a tetrahedral mesh (.ele).
676 // If no extension is available, the imposed command line switch
677 // (-p or -r) implies the object.
678 enum objecttype {NODES, POLY, OFF, PLY, STL, MEDIT, VTK, MESH} object;
679
680
681 void syntax();
682 void usage();
683
684 // Command line parse routine.
685 bool parse_commandline(int argc, char **argv);
686 bool parse_commandline(char *switches) {
687 return parse_commandline(0, &switches);
688 }
689
690 // Initialize all variables.
691 tetgenbehavior()
692 {
693 plc = 0;
694 psc = 0;
695 refine = 0;
696 quality = 0;
697 nobisect = 0;
698 coarsen = 0;
699 metric = 0;
700 weighted = 0;
701 brio_hilbert = 1;
702 incrflip = 0;
703 flipinsert = 0;
704 varvolume = 0;
705 fixedvolume = 0;
706 noexact = 0;
707 nostaticfilter = 0;
708 insertaddpoints = 0;
709 regionattrib = 0;
710 conforming = 0;
711 diagnose = 0;
712 convex = 0;
713 zeroindex = 0;
714 facesout = 0;
715 edgesout = 0;
716 neighout = 0;
717 voroout = 0;
718 meditview = 0;
719 vtkview = 0;
720 nobound = 0;
721 nonodewritten = 0;
722 noelewritten = 0;
723 nofacewritten = 0;
724 noiterationnum = 0;
725 nomergefacet = 0;
726 nomergevertex = 0;
727 nojettison = 0;
728 reversetetori = 0;
729 docheck = 0;
730 quiet = 0;
731 verbose = 0;
732
733 vertexperblock = 4092;
734 tetrahedraperblock = 8188;
735 shellfaceperblock = 4092;
736 nobisect_param = 2;
737 addsteiner_algo = 1;
738 coarsen_param = 0;
739 weighted_param = 0;
740 fliplinklevel = -1; // No limit on linklevel.
741 flipstarsize = -1; // No limit on flip star size.
742 fliplinklevelinc = 1;
743 reflevel = 3;
744 optscheme = 7; // 1 & 2 & 4, // min_max_dihedral.
745 optlevel = 2;
746 delmaxfliplevel = 1;
747 order = 1;
748 steinerleft = -1;
749 no_sort = 0;
750 hilbert_order = 52; //-1;
751 hilbert_limit = 8;
752 brio_threshold = 64;
753 brio_ratio = 0.125;
754 facet_ang_tol = 179.9;
755 maxvolume = -1.0;
756 minratio = 2.0;
757 mindihedral = 0.0; // 5.0;
758 optmaxdihedral = 165.00; // without -q, default is 179.0
759 optminsmtdihed = 179.00; // without -q, default is 179.999
760 optminslidihed = 179.00; // without -q, default is 179.999
761 epsilon = 1.0e-8;
762 minedgelength = 0.0;
763 coarsen_percent = 1.0;
764 object = NODES;
765
766 commandline[0] = '\0';
767 infilename[0] = '\0';
768 outfilename[0] = '\0';
769 addinfilename[0] = '\0';
770 bgmeshfilename[0] = '\0';
771
772 }
773
774}; // class tetgenbehavior
775
776///////////////////////////////////////////////////////////////////////////////
777// //
778// Robust Geometric predicates //
779// //
780// Geometric predicates are simple tests of spatial relations of a set of d- //
781// dimensional points, such as the orientation test and the point-in-sphere //
782// test. Each of these tests is performed by evaluating the sign of a deter- //
783// minant of a matrix whose entries are the coordinates of these points. If //
784// the computation is performed by using the floating-point numbers, e.g., //
785// the single or double precision numbers in C/C++, roundoff error may cause //
786// an incorrect result. This may either lead to a wrong result or eventually //
787// lead to a failure of the program. Computing the predicates exactly will //
788// avoid the error and make the program robust. //
789// //
790// The following routines are the robust geometric predicates for 3D orient- //
791// ation test and point-in-sphere test. They were implemented by Shewchuk. //
792// The source code are generously provided by him in the public domain, //
793// http://www.cs.cmu.edu/~quake/robust.html. predicates.cxx is a C++ version //
794// of the original C code. //
795// //
796// The original predicates of Shewchuk only use "dynamic filters", i.e., it //
797// computes the error at run time step by step. TetGen first adds a "static //
798// filter" in each predicate. It estimates the maximal possible error in all //
799// cases. So it can safely and quickly answer many easy cases. //
800// //
801///////////////////////////////////////////////////////////////////////////////
802
803void exactinit(int, int, int, REAL, REAL, REAL);
804REAL orient3d(REAL *pa, REAL *pb, REAL *pc, REAL *pd);
805REAL insphere(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL *pe);
806REAL orient4d(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL *pe,
807 REAL ah, REAL bh, REAL ch, REAL dh, REAL eh);
808
809///////////////////////////////////////////////////////////////////////////////
810// //
811// tetgenmesh //
812// //
813// A structure for creating and updating tetrahedral meshes. //
814// //
815///////////////////////////////////////////////////////////////////////////////
816
817class tetgenmesh {
818
819public:
820
821///////////////////////////////////////////////////////////////////////////////
822// //
823// Mesh data structure //
824// //
825// A tetrahedral mesh T of a 3D piecewise linear complex (PLC) X is a 3D //
826// simplicial complex whose underlying space is equal to the space of X. T //
827// contains a 2D subcomplex S which is a triangular mesh of the boundary of //
828// X. S contains a 1D subcomplex L which is a linear mesh of the boundary of //
829// S. Faces and edges in S and L are respectively called subfaces and segme- //
830// nts to distinguish them from others in T. //
831// //
832// TetGen stores the tetrahedra and vertices of T. The basic structure of a //
833// tetrahedron contains pointers to its vertices and adjacent tetrahedra. A //
834// vertex stores its x-, y-, and z-coordinates, and a pointer to a tetrahed- //
835// ron containing it. Both tetrahedra and vertices may contain user data. //
836// //
837// Each face of T belongs to either two tetrahedra or one tetrahedron. In //
838// the latter case, the face is an exterior boundary face of T. TetGen adds //
839// fictitious tetrahedra (one-to-one) at such faces, and connects them to an //
840// "infinite vertex" (which has no geometric coordinates). One can imagine //
841// such a vertex lies in 4D space and is visible by all exterior boundary //
842// faces. The extended set of tetrahedra (including the infinite vertex) is //
843// a tetrahedralization of a 3-pseudomanifold without boundary. It has the //
844// property that every face is shared by exactly two tetrahedra. //
845// //
846// The current version of TetGen stores explicitly the subfaces and segments //
847// (which are in surface mesh S and the linear mesh L), respectively. Extra //
848// pointers are allocated in tetrahedra and subfaces to point each others. //
849// //
850///////////////////////////////////////////////////////////////////////////////
851
852 // The tetrahedron data structure. It includes the following fields:
853 // - a list of four adjoining tetrahedra;
854 // - a list of four vertices;
855 // - a pointer to a list of four subfaces (optional, for -p switch);
856 // - a pointer to a list of six segments (optional, for -p switch);
857 // - a list of user-defined floating-point attributes (optional);
858 // - a volume constraint (optional, for -a switch);
859 // - an integer of element marker (and flags);
860 // The structure of a tetrahedron is an array of pointers. Its actual size
861 // (the length of the array) is determined at runtime.
862
863 typedef REAL **tetrahedron;
864
865 // The subface data structure. It includes the following fields:
866 // - a list of three adjoining subfaces;
867 // - a list of three vertices;
868 // - a list of three adjoining segments;
869 // - two adjoining tetrahedra;
870 // - an area constraint (optional, for -q switch);
871 // - an integer for boundary marker;
872 // - an integer for type, flags, etc.
873
874 typedef REAL **shellface;
875
876 // The point data structure. It includes the following fields:
877 // - x, y and z coordinates;
878 // - a list of user-defined point attributes (optional);
879 // - u, v coordinates (optional, for -s switch);
880 // - a metric tensor (optional, for -q or -m switch);
881 // - a pointer to an adjacent tetrahedron;
882 // - a pointer to a parent (or a duplicate) point;
883 // - a pointer to an adjacent subface or segment (optional, -p switch);
884 // - a pointer to a tet in background mesh (optional, for -m switch);
885 // - an integer for boundary marker (point index);
886 // - an integer for point type (and flags).
887 // - an integer for geometry tag (optional, for -s switch).
888 // The structure of a point is an array of REALs. Its acutal size is
889 // determined at the runtime.
890
891 typedef REAL *point;
892
893///////////////////////////////////////////////////////////////////////////////
894// //
895// Handles //
896// //
897// Navigation and manipulation in a tetrahedralization are accomplished by //
898// operating on structures referred as ``handles". A handle is a pair (t,v), //
899// where t is a pointer to a tetrahedron, and v is a 4-bit integer, in the //
900// range from 0 to 11. v is called the ``version'' of a tetrahedron, it rep- //
901// resents a directed edge of a specific face of the tetrahedron. //
902// //
903// There are 12 even permutations of the four vertices, each of them corres- //
904// ponds to a directed edge (a version) of the tetrahedron. The 12 versions //
905// can be grouped into 4 distinct ``edge rings'' in 4 ``oriented faces'' of //
906// this tetrahedron. One can encode each version (a directed edge) into a //
907// 4-bit integer such that the two upper bits encode the index (from 0 to 2) //
908// of this edge in the edge ring, and the two lower bits encode the index ( //
909// from 0 to 3) of the oriented face which contains this edge. //
910// //
911// The four vertices of a tetrahedron are indexed from 0 to 3 (according to //
912// their storage in the data structure). Give each face the same index as //
913// the node opposite it in the tetrahedron. Denote the edge connecting face //
914// i to face j as i/j. We number the twelve versions as follows: //
915// //
916// | edge 0 edge 1 edge 2 //
917// --------|-------------------------------- //
918// face 0 | 0 (0/1) 4 (0/3) 8 (0/2) //
919// face 1 | 1 (1/2) 5 (1/3) 9 (1/0) //
920// face 2 | 2 (2/3) 6 (2/1) 10 (2/0) //
921// face 3 | 3 (3/0) 7 (3/1) 11 (3/2) //
922// //
923// Similarly, navigation and manipulation in a (boundary) triangulation are //
924// done by using handles of triangles. Each handle is a pair (s, v), where s //
925// is a pointer to a triangle, and v is a version in the range from 0 to 5. //
926// Each version corresponds to a directed edge of this triangle. //
927// //
928// Number the three vertices of a triangle from 0 to 2 (according to their //
929// storage in the data structure). Give each edge the same index as the node //
930// opposite it in the triangle. The six versions of a triangle are: //
931// //
932// | edge 0 edge 1 edge 2 //
933// ---------------|-------------------------- //
934// ccw orieation | 0 2 4 //
935// cw orieation | 1 3 5 //
936// //
937// In the following, a 'triface' is a handle of tetrahedron, and a 'face' is //
938// a handle of a triangle. //
939// //
940///////////////////////////////////////////////////////////////////////////////
941
942 class triface {
943 public:
944 tetrahedron *tet;
945 int ver; // Range from 0 to 11.
946 triface() : tet(0), ver(0) {}
947 triface& operator=(const triface& t) {
948 tet = t.tet; ver = t.ver;
949 return *this;
950 }
951 };
952
953 class face {
954 public:
955 shellface *sh;
956 int shver; // Range from 0 to 5.
957 face() : sh(0), shver(0) {}
958 face& operator=(const face& s) {
959 sh = s.sh; shver = s.shver;
960 return *this;
961 }
962 };
963
964///////////////////////////////////////////////////////////////////////////////
965// //
966// Arraypool //
967// //
968// A dynamic linear array. (It is written by J. Shewchuk) //
969// //
970// Each arraypool contains an array of pointers to a number of blocks. Each //
971// block contains the same fixed number of objects. Each index of the array //
972// addresses a particular object in the pool. The most significant bits add- //
973// ress the index of the block containing the object. The less significant //
974// bits address this object within the block. //
975// //
976// 'objectbytes' is the size of one object in blocks; 'log2objectsperblock' //
977// is the base-2 logarithm of 'objectsperblock'; 'objects' counts the number //
978// of allocated objects; 'totalmemory' is the total memory in bytes. //
979// //
980///////////////////////////////////////////////////////////////////////////////
981
982 class arraypool {
983
984 public:
985
986 int objectbytes;
987 int objectsperblock;
988 int log2objectsperblock;
989 int objectsperblockmark;
990 int toparraylen;
991 char **toparray;
992 long objects;
993 unsigned long totalmemory;
994
995 void restart();
996 void poolinit(int sizeofobject, int log2objperblk);
997 char* getblock(int objectindex);
998 void* lookup(int objectindex);
999 int newindex(void **newptr);
1000
1001 arraypool(int sizeofobject, int log2objperblk);
1002 ~arraypool();
1003 };
1004
1005// fastlookup() -- A fast, unsafe operation. Return the pointer to the object
1006// with a given index. Note: The object's block must have been allocated,
1007// i.e., by the function newindex().
1008
1009#define fastlookup(pool, index) \
1010 (void *) ((pool)->toparray[(index) >> (pool)->log2objectsperblock] + \
1011 ((index) & (pool)->objectsperblockmark) * (pool)->objectbytes)
1012
1013///////////////////////////////////////////////////////////////////////////////
1014// //
1015// Memorypool //
1016// //
1017// A structure for memory allocation. (It is written by J. Shewchuk) //
1018// //
1019// firstblock is the first block of items. nowblock is the block from which //
1020// items are currently being allocated. nextitem points to the next slab //
1021// of free memory for an item. deaditemstack is the head of a linked list //
1022// (stack) of deallocated items that can be recycled. unallocateditems is //
1023// the number of items that remain to be allocated from nowblock. //
1024// //
1025// Traversal is the process of walking through the entire list of items, and //
1026// is separate from allocation. Note that a traversal will visit items on //
1027// the "deaditemstack" stack as well as live items. pathblock points to //
1028// the block currently being traversed. pathitem points to the next item //
1029// to be traversed. pathitemsleft is the number of items that remain to //
1030// be traversed in pathblock. //
1031// //
1032///////////////////////////////////////////////////////////////////////////////
1033
1034 class memorypool {
1035
1036 public:
1037
1038 void **firstblock, **nowblock;
1039 void *nextitem;
1040 void *deaditemstack;
1041 void **pathblock;
1042 void *pathitem;
1043 int alignbytes;
1044 int itembytes, itemwords;
1045 int itemsperblock;
1046 long items, maxitems;
1047 int unallocateditems;
1048 int pathitemsleft;
1049
1050 memorypool();
1051 memorypool(int, int, int, int);
1052 ~memorypool();
1053
1054 void poolinit(int, int, int, int);
1055 void restart();
1056 void *alloc();
1057 void dealloc(void*);
1058 void traversalinit();
1059 void *traverse();
1060 };
1061
1062///////////////////////////////////////////////////////////////////////////////
1063// //
1064// badface //
1065// //
1066// Despite of its name, a 'badface' can be used to represent one of the //
1067// following objects: //
1068// - a face of a tetrahedron which is (possibly) non-Delaunay; //
1069// - an encroached subsegment or subface; //
1070// - a bad-quality tetrahedron, i.e, has too large radius-edge ratio; //
1071// - a sliver, i.e., has good radius-edge ratio but nearly zero volume; //
1072// - a recently flipped face (saved for undoing the flip later). //
1073// //
1074///////////////////////////////////////////////////////////////////////////////
1075
1076 class badface {
1077 public:
1078 triface tt;
1079 face ss;
1080 REAL key, cent[6]; // circumcenter or cos(dihedral angles) at 6 edges.
1081 point forg, fdest, fapex, foppo, noppo;
1082 badface *nextitem;
1083 badface() : key(0), forg(0), fdest(0), fapex(0), foppo(0), noppo(0),
1084 nextitem(0) {}
1085 };
1086
1087///////////////////////////////////////////////////////////////////////////////
1088// //
1089// insertvertexflags //
1090// //
1091// A collection of flags that pass to the routine insertvertex(). //
1092// //
1093///////////////////////////////////////////////////////////////////////////////
1094
1095 class insertvertexflags {
1096
1097 public:
1098
1099 int iloc; // input/output.
1100 int bowywat, lawson;
1101 int splitbdflag, validflag, respectbdflag;
1102 int rejflag, chkencflag, cdtflag;
1103 int assignmeshsize;
1104 int sloc, sbowywat;
1105
1106 // Used by Delaunay refinement.
1107 int refineflag; // 0, 1, 2, 3
1108 triface refinetet;
1109 face refinesh;
1110 int smlenflag; // for useinsertradius.
1111 REAL smlen; // for useinsertradius.
1112 point parentpt;
1113
1114 insertvertexflags() {
1115 iloc = bowywat = lawson = 0;
1116 splitbdflag = validflag = respectbdflag = 0;
1117 rejflag = chkencflag = cdtflag = 0;
1118 assignmeshsize = 0;
1119 sloc = sbowywat = 0;
1120
1121 refineflag = 0;
1122 refinetet.tet = NULL;
1123 refinesh.sh = NULL;
1124 smlenflag = 0;
1125 smlen = 0.0;
1126 }
1127 };
1128
1129///////////////////////////////////////////////////////////////////////////////
1130// //
1131// flipconstraints //
1132// //
1133// A structure of a collection of data (options and parameters) which pass //
1134// to the edge flip function flipnm(). //
1135// //
1136///////////////////////////////////////////////////////////////////////////////
1137
1138 class flipconstraints {
1139
1140 public:
1141
1142 // Elementary flip flags.
1143 int enqflag; // (= flipflag)
1144 int chkencflag;
1145
1146 // Control flags
1147 int unflip; // Undo the performed flips.
1148 int collectnewtets; // Collect the new tets created by flips.
1149 int collectencsegflag;
1150
1151 // Optimization flags.
1152 int remove_ndelaunay_edge; // Remove a non-Delaunay edge.
1153 REAL bak_tetprism_vol; // The value to be minimized.
1154 REAL tetprism_vol_sum;
1155 int remove_large_angle; // Remove a large dihedral angle at edge.
1156 REAL cosdihed_in; // The input cosine of the dihedral angle (> 0).
1157 REAL cosdihed_out; // The improved cosine of the dihedral angle.
1158
1159 // Boundary recovery flags.
1160 int checkflipeligibility;
1161 point seg[2]; // A constraining edge to be recovered.
1162 point fac[3]; // A constraining face to be recovered.
1163 point remvert; // A vertex to be removed.
1164
1165
1166 flipconstraints() {
1167 enqflag = 0;
1168 chkencflag = 0;
1169
1170 unflip = 0;
1171 collectnewtets = 0;
1172 collectencsegflag = 0;
1173
1174 remove_ndelaunay_edge = 0;
1175 bak_tetprism_vol = 0.0;
1176 tetprism_vol_sum = 0.0;
1177 remove_large_angle = 0;
1178 cosdihed_in = 0.0;
1179 cosdihed_out = 0.0;
1180
1181 checkflipeligibility = 0;
1182 seg[0] = NULL;
1183 fac[0] = NULL;
1184 remvert = NULL;
1185 }
1186 };
1187
1188///////////////////////////////////////////////////////////////////////////////
1189// //
1190// optparameters //
1191// //
1192// Optimization options and parameters. //
1193// //
1194///////////////////////////////////////////////////////////////////////////////
1195
1196 class optparameters {
1197
1198 public:
1199
1200 // The one of goals of optimization.
1201 int max_min_volume; // Maximize the minimum volume.
1202 int max_min_aspectratio; // Maximize the minimum aspect ratio.
1203 int min_max_dihedangle; // Minimize the maximum dihedral angle.
1204
1205 // The initial and improved value.
1206 REAL initval, imprval;
1207
1208 int numofsearchdirs;
1209 REAL searchstep;
1210 int maxiter; // Maximum smoothing iterations (disabled by -1).
1211 int smthiter; // Performed iterations.
1212
1213
1214 optparameters() {
1215 max_min_volume = 0;
1216 max_min_aspectratio = 0;
1217 min_max_dihedangle = 0;
1218
1219 initval = imprval = 0.0;
1220
1221 numofsearchdirs = 10;
1222 searchstep = 0.01;
1223 maxiter = -1; // Unlimited smoothing iterations.
1224 smthiter = 0;
1225
1226 }
1227 };
1228
1229
1230///////////////////////////////////////////////////////////////////////////////
1231// //
1232// Labels (enumeration declarations) used by TetGen. //
1233// //
1234///////////////////////////////////////////////////////////////////////////////
1235
1236 // Labels that signify the type of a vertex.
1237 enum verttype {UNUSEDVERTEX, DUPLICATEDVERTEX, RIDGEVERTEX, ACUTEVERTEX,
1238 FACETVERTEX, VOLVERTEX, FREESEGVERTEX, FREEFACETVERTEX,
1239 FREEVOLVERTEX, NREGULARVERTEX, DEADVERTEX};
1240
1241 // Labels that signify the result of triangle-triangle intersection test.
1242 enum interresult {DISJOINT, INTERSECT, SHAREVERT, SHAREEDGE, SHAREFACE,
1243 TOUCHEDGE, TOUCHFACE, ACROSSVERT, ACROSSEDGE, ACROSSFACE,
1244 COLLISIONFACE, ACROSSSEG, ACROSSSUB};
1245
1246 // Labels that signify the result of point location.
1247 enum locateresult {UNKNOWN, OUTSIDE, INTETRAHEDRON, ONFACE, ONEDGE, ONVERTEX,
1248 ENCVERTEX, ENCSEGMENT, ENCSUBFACE, NEARVERTEX, NONREGULAR,
1249 INSTAR, BADELEMENT};
1250
1251///////////////////////////////////////////////////////////////////////////////
1252// //
1253// Variables of TetGen //
1254// //
1255///////////////////////////////////////////////////////////////////////////////
1256
1257 // Pointer to the input data (a set of nodes, a PLC, or a mesh).
1258 tetgenio *in, *addin;
1259
1260 // Pointer to the switches and parameters.
1261 tetgenbehavior *b;
1262
1263 // Pointer to a background mesh (contains size specification map).
1264 tetgenmesh *bgm;
1265
1266 // Memorypools to store mesh elements (points, tetrahedra, subfaces, and
1267 // segments) and extra pointers between tetrahedra, subfaces, and segments.
1268 memorypool *tetrahedrons, *subfaces, *subsegs, *points;
1269 memorypool *tet2subpool, *tet2segpool;
1270
1271 // Memorypools to store bad-quality (or encroached) elements.
1272 memorypool *badtetrahedrons, *badsubfacs, *badsubsegs;
1273
1274 // A memorypool to store faces to be flipped.
1275 memorypool *flippool;
1276 arraypool *unflipqueue;
1277 badface *flipstack;
1278
1279 // Arrays used for point insertion (the Bowyer-Watson algorithm).
1280 arraypool *cavetetlist, *cavebdrylist, *caveoldtetlist;
1281 arraypool *cavetetshlist, *cavetetseglist, *cavetetvertlist;
1282 arraypool *caveencshlist, *caveencseglist;
1283 arraypool *caveshlist, *caveshbdlist, *cavesegshlist;
1284
1285 // Stacks used for CDT construction and boundary recovery.
1286 arraypool *subsegstack, *subfacstack, *subvertstack;
1287
1288 // Arrays of encroached segments and subfaces (for mesh refinement).
1289 arraypool *encseglist, *encshlist;
1290
1291 // The map between facets to their vertices (for mesh refinement).
1292 int *idx2facetlist;
1293 point *facetverticeslist;
1294
1295 // The map between segments to their endpoints (for mesh refinement).
1296 point *segmentendpointslist;
1297
1298 // The infinite vertex.
1299 point dummypoint;
1300 // The recently visited tetrahedron, subface.
1301 triface recenttet;
1302 face recentsh;
1303
1304 // PI is the ratio of a circle's circumference to its diameter.
1305 static REAL PI;
1306
1307 // Array (size = numberoftetrahedra * 6) for storing high-order nodes of
1308 // tetrahedra (only used when -o2 switch is selected).
1309 point *highordertable;
1310
1311 // Various variables.
1312 int numpointattrib; // Number of point attributes.
1313 int numelemattrib; // Number of tetrahedron attributes.
1314 int sizeoftensor; // Number of REALs per metric tensor.
1315 int pointmtrindex; // Index to find the metric tensor of a point.
1316 int pointparamindex; // Index to find the u,v coordinates of a point.
1317 int point2simindex; // Index to find a simplex adjacent to a point.
1318 int pointmarkindex; // Index to find boundary marker of a point.
1319 int elemattribindex; // Index to find attributes of a tetrahedron.
1320 int volumeboundindex; // Index to find volume bound of a tetrahedron.
1321 int elemmarkerindex; // Index to find marker of a tetrahedron.
1322 int shmarkindex; // Index to find boundary marker of a subface.
1323 int areaboundindex; // Index to find area bound of a subface.
1324 int checksubsegflag; // Are there segments in the tetrahedralization yet?
1325 int checksubfaceflag; // Are there subfaces in the tetrahedralization yet?
1326 int checkconstraints; // Are there variant (node, seg, facet) constraints?
1327 int nonconvex; // Is current mesh non-convex?
1328 int autofliplinklevel; // The increase of link levels, default is 1.
1329 int useinsertradius; // Save the insertion radius for Steiner points.
1330 long samples; // Number of random samples for point location.
1331 unsigned long randomseed; // Current random number seed.
1332 REAL cosmaxdihed, cosmindihed; // The cosine values of max/min dihedral.
1333 REAL cossmtdihed; // The cosine value of a bad dihedral to be smoothed.
1334 REAL cosslidihed; // The cosine value of the max dihedral of a sliver.
1335 REAL minfaceang, minfacetdihed; // The minimum input (dihedral) angles.
1336 REAL tetprism_vol_sum; // The total volume of tetrahedral-prisms (in 4D).
1337 REAL longest; // The longest possible edge length.
1338 REAL xmax, xmin, ymax, ymin, zmax, zmin; // Bounding box of points.
1339
1340 // Counters.
1341 long insegments; // Number of input segments.
1342 long hullsize; // Number of exterior boundary faces.
1343 long meshedges; // Number of mesh edges.
1344 long meshhulledges; // Number of boundary mesh edges.
1345 long steinerleft; // Number of Steiner points not yet used.
1346 long dupverts; // Are there duplicated vertices?
1347 long unuverts; // Are there unused vertices?
1348 long nonregularcount; // Are there non-regular vertices?
1349 long st_segref_count, st_facref_count, st_volref_count; // Steiner points.
1350 long fillregioncount, cavitycount, cavityexpcount;
1351 long flip14count, flip26count, flipn2ncount;
1352 long flip23count, flip32count, flip44count, flip41count;
1353 long flip31count, flip22count;
1354 unsigned long totalworkmemory; // Total memory used by working arrays.
1355
1356
1357///////////////////////////////////////////////////////////////////////////////
1358// //
1359// Mesh manipulation primitives //
1360// //
1361///////////////////////////////////////////////////////////////////////////////
1362
1363 // Fast lookup tables for mesh manipulation primitives.
1364 static int bondtbl[12][12], fsymtbl[12][12];
1365 static int esymtbl[12], enexttbl[12], eprevtbl[12];
1366 static int enextesymtbl[12], eprevesymtbl[12];
1367 static int eorgoppotbl[12], edestoppotbl[12];
1368 static int facepivot1[12], facepivot2[12][12];
1369 static int orgpivot[12], destpivot[12], apexpivot[12], oppopivot[12];
1370 static int tsbondtbl[12][6], stbondtbl[12][6];
1371 static int tspivottbl[12][6], stpivottbl[12][6];
1372 static int ver2edge[12], edge2ver[6], epivot[12];
1373 static int sorgpivot [6], sdestpivot[6], sapexpivot[6];
1374 static int snextpivot[6];
1375
1376 void inittables();
1377
1378 // Primitives for tetrahedra.
1379 inline tetrahedron encode(triface& t);
1380 inline tetrahedron encode2(tetrahedron* ptr, int ver);
1381 inline void decode(tetrahedron ptr, triface& t);
1382 inline void bond(triface& t1, triface& t2);
1383 inline void dissolve(triface& t);
1384 inline void esym(triface& t1, triface& t2);
1385 inline void esymself(triface& t);
1386 inline void enext(triface& t1, triface& t2);
1387 inline void enextself(triface& t);
1388 inline void eprev(triface& t1, triface& t2);
1389 inline void eprevself(triface& t);
1390 inline void enextesym(triface& t1, triface& t2);
1391 inline void enextesymself(triface& t);
1392 inline void eprevesym(triface& t1, triface& t2);
1393 inline void eprevesymself(triface& t);
1394 inline void eorgoppo(triface& t1, triface& t2);
1395 inline void eorgoppoself(triface& t);
1396 inline void edestoppo(triface& t1, triface& t2);
1397 inline void edestoppoself(triface& t);
1398 inline void fsym(triface& t1, triface& t2);
1399 inline void fsymself(triface& t);
1400 inline void fnext(triface& t1, triface& t2);
1401 inline void fnextself(triface& t);
1402 inline point org (triface& t);
1403 inline point dest(triface& t);
1404 inline point apex(triface& t);
1405 inline point oppo(triface& t);
1406 inline void setorg (triface& t, point p);
1407 inline void setdest(triface& t, point p);
1408 inline void setapex(triface& t, point p);
1409 inline void setoppo(triface& t, point p);
1410 inline REAL elemattribute(tetrahedron* ptr, int attnum);
1411 inline void setelemattribute(tetrahedron* ptr, int attnum, REAL value);
1412 inline REAL volumebound(tetrahedron* ptr);
1413 inline void setvolumebound(tetrahedron* ptr, REAL value);
1414 inline int elemindex(tetrahedron* ptr);
1415 inline void setelemindex(tetrahedron* ptr, int value);
1416 inline int elemmarker(tetrahedron* ptr);
1417 inline void setelemmarker(tetrahedron* ptr, int value);
1418 inline void infect(triface& t);
1419 inline void uninfect(triface& t);
1420 inline bool infected(triface& t);
1421 inline void marktest(triface& t);
1422 inline void unmarktest(triface& t);
1423 inline bool marktested(triface& t);
1424 inline void markface(triface& t);
1425 inline void unmarkface(triface& t);
1426 inline bool facemarked(triface& t);
1427 inline void markedge(triface& t);
1428 inline void unmarkedge(triface& t);
1429 inline bool edgemarked(triface& t);
1430 inline void marktest2(triface& t);
1431 inline void unmarktest2(triface& t);
1432 inline bool marktest2ed(triface& t);
1433 inline int elemcounter(triface& t);
1434 inline void setelemcounter(triface& t, int value);
1435 inline void increaseelemcounter(triface& t);
1436 inline void decreaseelemcounter(triface& t);
1437 inline bool ishulltet(triface& t);
1438 inline bool isdeadtet(triface& t);
1439
1440 // Primitives for subfaces and subsegments.
1441 inline void sdecode(shellface sptr, face& s);
1442 inline shellface sencode(face& s);
1443 inline shellface sencode2(shellface *sh, int shver);
1444 inline void spivot(face& s1, face& s2);
1445 inline void spivotself(face& s);
1446 inline void sbond(face& s1, face& s2);
1447 inline void sbond1(face& s1, face& s2);
1448 inline void sdissolve(face& s);
1449 inline point sorg(face& s);
1450 inline point sdest(face& s);
1451 inline point sapex(face& s);
1452 inline void setsorg(face& s, point pointptr);
1453 inline void setsdest(face& s, point pointptr);
1454 inline void setsapex(face& s, point pointptr);
1455 inline void sesym(face& s1, face& s2);
1456 inline void sesymself(face& s);
1457 inline void senext(face& s1, face& s2);
1458 inline void senextself(face& s);
1459 inline void senext2(face& s1, face& s2);
1460 inline void senext2self(face& s);
1461 inline REAL areabound(face& s);
1462 inline void setareabound(face& s, REAL value);
1463 inline int shellmark(face& s);
1464 inline void setshellmark(face& s, int value);
1465 inline void sinfect(face& s);
1466 inline void suninfect(face& s);
1467 inline bool sinfected(face& s);
1468 inline void smarktest(face& s);
1469 inline void sunmarktest(face& s);
1470 inline bool smarktested(face& s);
1471 inline void smarktest2(face& s);
1472 inline void sunmarktest2(face& s);
1473 inline bool smarktest2ed(face& s);
1474 inline void smarktest3(face& s);
1475 inline void sunmarktest3(face& s);
1476 inline bool smarktest3ed(face& s);
1477 inline void setfacetindex(face& f, int value);
1478 inline int getfacetindex(face& f);
1479
1480 // Primitives for interacting tetrahedra and subfaces.
1481 inline void tsbond(triface& t, face& s);
1482 inline void tsdissolve(triface& t);
1483 inline void stdissolve(face& s);
1484 inline void tspivot(triface& t, face& s);
1485 inline void stpivot(face& s, triface& t);
1486
1487 // Primitives for interacting tetrahedra and segments.
1488 inline void tssbond1(triface& t, face& seg);
1489 inline void sstbond1(face& s, triface& t);
1490 inline void tssdissolve1(triface& t);
1491 inline void sstdissolve1(face& s);
1492 inline void tsspivot1(triface& t, face& s);
1493 inline void sstpivot1(face& s, triface& t);
1494
1495 // Primitives for interacting subfaces and segments.
1496 inline void ssbond(face& s, face& edge);
1497 inline void ssbond1(face& s, face& edge);
1498 inline void ssdissolve(face& s);
1499 inline void sspivot(face& s, face& edge);
1500
1501 // Primitives for points.
1502 inline int pointmark(point pt);
1503 inline void setpointmark(point pt, int value);
1504 inline enum verttype pointtype(point pt);
1505 inline void setpointtype(point pt, enum verttype value);
1506 inline int pointgeomtag(point pt);
1507 inline void setpointgeomtag(point pt, int value);
1508 inline REAL pointgeomuv(point pt, int i);
1509 inline void setpointgeomuv(point pt, int i, REAL value);
1510 inline void pinfect(point pt);
1511 inline void puninfect(point pt);
1512 inline bool pinfected(point pt);
1513 inline void pmarktest(point pt);
1514 inline void punmarktest(point pt);
1515 inline bool pmarktested(point pt);
1516 inline void pmarktest2(point pt);
1517 inline void punmarktest2(point pt);
1518 inline bool pmarktest2ed(point pt);
1519 inline void pmarktest3(point pt);
1520 inline void punmarktest3(point pt);
1521 inline bool pmarktest3ed(point pt);
1522 inline tetrahedron point2tet(point pt);
1523 inline void setpoint2tet(point pt, tetrahedron value);
1524 inline shellface point2sh(point pt);
1525 inline void setpoint2sh(point pt, shellface value);
1526 inline point point2ppt(point pt);
1527 inline void setpoint2ppt(point pt, point value);
1528 inline tetrahedron point2bgmtet(point pt);
1529 inline void setpoint2bgmtet(point pt, tetrahedron value);
1530 inline void setpointinsradius(point pt, REAL value);
1531 inline REAL getpointinsradius(point pt);
1532
1533 // Advanced primitives.
1534 inline void point2tetorg(point pt, triface& t);
1535 inline void point2shorg(point pa, face& s);
1536 inline point farsorg(face& seg);
1537 inline point farsdest(face& seg);
1538
1539///////////////////////////////////////////////////////////////////////////////
1540// //
1541// Memory managment //
1542// //
1543///////////////////////////////////////////////////////////////////////////////
1544
1545 void tetrahedrondealloc(tetrahedron*);
1546 tetrahedron *tetrahedrontraverse();
1547 tetrahedron *alltetrahedrontraverse();
1548 void shellfacedealloc(memorypool*, shellface*);
1549 shellface *shellfacetraverse(memorypool*);
1550 void pointdealloc(point);
1551 point pointtraverse();
1552
1553 void makeindex2pointmap(point*&);
1554 void makepoint2submap(memorypool*, int*&, face*&);
1555 void maketetrahedron(triface*);
1556 void makeshellface(memorypool*, face*);
1557 void makepoint(point*, enum verttype);
1558
1559 void initializepools();
1560
1561///////////////////////////////////////////////////////////////////////////////
1562// //
1563// Advanced geometric predicates and calculations //
1564// //
1565// TetGen uses a simplified symbolic perturbation scheme from Edelsbrunner, //
1566// et al [*]. Hence the point-in-sphere test never returns a zero. The idea //
1567// is to perturb the weights of vertices in the fourth dimension. TetGen //
1568// uses the indices of the vertices decide the amount of perturbation. It is //
1569// implemented in the routine insphere_s().
1570// //
1571// The routine tri_edge_test() determines whether or not a triangle and an //
1572// edge intersect in 3D. If they intersect, their intersection type is also //
1573// reported. This test is a combination of n 3D orientation tests (n is bet- //
1574// ween 3 and 9). It uses the robust orient3d() test to make the branch dec- //
1575// isions. The routine tri_tri_test() determines whether or not two triang- //
1576// les intersect in 3D. It also uses the robust orient3d() test. //
1577// //
1578// There are a number of routines to calculate geometrical quantities, e.g., //
1579// circumcenters, angles, dihedral angles, face normals, face areas, etc. //
1580// They are so far done by the default floating-point arithmetics which are //
1581// non-robust. They should be improved in the future. //
1582// //
1583///////////////////////////////////////////////////////////////////////////////
1584
1585 // Symbolic perturbations (robust)
1586 REAL insphere_s(REAL*, REAL*, REAL*, REAL*, REAL*);
1587 REAL orient4d_s(REAL*, REAL*, REAL*, REAL*, REAL*,
1588 REAL, REAL, REAL, REAL, REAL);
1589
1590 // Triangle-edge intersection test (robust)
1591 int tri_edge_2d(point, point, point, point, point, point, int, int*, int*);
1592 int tri_edge_tail(point, point, point, point, point, point, REAL, REAL, int,
1593 int*, int*);
1594 int tri_edge_test(point, point, point, point, point, point, int, int*, int*);
1595
1596 // Triangle-triangle intersection test (robust)
1597 int tri_edge_inter_tail(point, point, point, point, point, REAL, REAL);
1598 int tri_tri_inter(point, point, point, point, point, point);
1599
1600 // Linear algebra functions
1601 inline REAL dot(REAL* v1, REAL* v2);
1602 inline void cross(REAL* v1, REAL* v2, REAL* n);
1603 bool lu_decmp(REAL lu[4][4], int n, int* ps, REAL* d, int N);
1604 void lu_solve(REAL lu[4][4], int n, int* ps, REAL* b, int N);
1605
1606 // An embedded 2-dimensional geometric predicate (non-robust)
1607 REAL incircle3d(point pa, point pb, point pc, point pd);
1608
1609 // Geometric calculations (non-robust)
1610 REAL orient3dfast(REAL *pa, REAL *pb, REAL *pc, REAL *pd);
1611 inline REAL norm2(REAL x, REAL y, REAL z);
1612 inline REAL distance(REAL* p1, REAL* p2);
1613 void facenormal(point pa, point pb, point pc, REAL *n, int pivot, REAL *lav);
1614 REAL shortdistance(REAL* p, REAL* e1, REAL* e2);
1615 REAL triarea(REAL* pa, REAL* pb, REAL* pc);
1616 REAL interiorangle(REAL* o, REAL* p1, REAL* p2, REAL* n);
1617 void projpt2edge(REAL* p, REAL* e1, REAL* e2, REAL* prj);
1618 void projpt2face(REAL* p, REAL* f1, REAL* f2, REAL* f3, REAL* prj);
1619 REAL facedihedral(REAL* pa, REAL* pb, REAL* pc1, REAL* pc2);
1620 bool tetalldihedral(point, point, point, point, REAL*, REAL*, REAL*);
1621 void tetallnormal(point, point, point, point, REAL N[4][3], REAL* volume);
1622 REAL tetaspectratio(point, point, point, point);
1623 bool circumsphere(REAL*, REAL*, REAL*, REAL*, REAL* cent, REAL* radius);
1624 bool orthosphere(REAL*,REAL*,REAL*,REAL*,REAL,REAL,REAL,REAL,REAL*,REAL*);
1625 void planelineint(REAL*, REAL*, REAL*, REAL*, REAL*, REAL*, REAL*);
1626 int linelineint(REAL*, REAL*, REAL*, REAL*, REAL*, REAL*, REAL*, REAL*);
1627 REAL tetprismvol(REAL* pa, REAL* pb, REAL* pc, REAL* pd);
1628 bool calculateabovepoint(arraypool*, point*, point*, point*);
1629 void calculateabovepoint4(point, point, point, point);
1630
1631///////////////////////////////////////////////////////////////////////////////
1632// //
1633// Local mesh transformations //
1634// //
1635// A local transformation replaces a small set of tetrahedra with another //
1636// set of tetrahedra which fills the same space and the same boundaries. //
1637// In 3D, the most simplest local transformations are the elementary flips //
1638// performed within the convex hull of five vertices: 2-to-3, 3-to-2, 1-to-4,//
1639// and 4-to-1 flips, where the numbers indicate the number of tetrahedra //
1640// before and after each flip. The 1-to-4 and 4-to-1 flip involve inserting //
1641// or deleting a vertex, respectively. //
1642// There are complex local transformations which can be decomposed as a //
1643// combination of elementary flips. For example,a 4-to-4 flip which replaces //
1644// two coplanar edges can be regarded by a 2-to-3 flip and a 3-to-2 flip. //
1645// Note that the first 2-to-3 flip will temporarily create a degenerate tet- //
1646// rahedron which is removed immediately by the followed 3-to-2 flip. More //
1647// generally, a n-to-m flip, where n > 3, m = (n - 2) * 2, which removes an //
1648// edge can be done by first performing a sequence of (n - 3) 2-to-3 flips //
1649// followed by a 3-to-2 flip. //
1650// //
1651// The routines flip23(), flip32(), and flip41() perform the three element- //
1652// ray flips. The flip14() is available inside the routine insertpoint(). //
1653// //
1654// The routines flipnm() and flipnm_post() implement a generalized edge flip //
1655// algorithm which uses a combination of elementary flips. //
1656// //
1657// The routine insertpoint() implements a variant of Bowyer-Watson's cavity //
1658// algorithm to insert a vertex. It works for arbitrary tetrahedralization, //
1659// either Delaunay, or constrained Delaunay, or non-Delaunay. //
1660// //
1661///////////////////////////////////////////////////////////////////////////////
1662
1663 // The elementary flips.
1664 void flip23(triface*, int, flipconstraints* fc);
1665 void flip32(triface*, int, flipconstraints* fc);
1666 void flip41(triface*, int, flipconstraints* fc);
1667
1668 // A generalized edge flip.
1669 int flipnm(triface*, int n, int level, int, flipconstraints* fc);
1670 int flipnm_post(triface*, int n, int nn, int, flipconstraints* fc);
1671
1672 // Point insertion.
1673 int insertpoint(point, triface*, face*, face*, insertvertexflags*);
1674 void insertpoint_abort(face*, insertvertexflags*);
1675
1676///////////////////////////////////////////////////////////////////////////////
1677// //
1678// Delaunay tetrahedralization //
1679// //
1680// The routine incrementaldelaunay() implemented two incremental algorithms //
1681// for constructing Delaunay tetrahedralizations (DTs): the Bowyer-Watson //
1682// (B-W) algorithm and the incremental flip algorithm of Edelsbrunner and //
1683// Shah, "Incremental topological flipping works for regular triangulation," //
1684// Algorithmica, 15:233-241, 1996. //
1685// //
1686// The routine incrementalflip() implements the flip algorithm of [Edelsbru- //
1687// nner and Shah, 1996]. It flips a queue of locally non-Delaunay faces (in //
1688// an arbitrary order). The success is guaranteed when the Delaunay tetrah- //
1689// edralization is constructed incrementally by adding one vertex at a time. //
1690// //
1691// The routine locate() finds a tetrahedron contains a new point in current //
1692// DT. It uses a simple stochastic walk algorithm: starting from an arbitr- //
1693// ary tetrahedron in DT, it finds the destination by visit one tetrahedron //
1694// at a time, randomly chooses a tetrahedron if there are more than one //
1695// choices. This algorithm terminates due to Edelsbrunner's acyclic theorem. //
1696// Choose a good starting tetrahedron is crucial to the speed of the walk. //
1697// TetGen originally uses the "jump-and-walk" algorithm of Muecke, E.P., //
1698// Saias, I., and Zhu, B. "Fast Randomized Point Location Without Preproces- //
1699// sing." In Proceedings of the 12th ACM Symposium on Computational Geometry,//
1700// 274-283, 1996. It first randomly samples several tetrahedra in the DT //
1701// and then choosing the closet one to start walking. //
1702// The above algorithm slows download dramatically as the number of points //
1703// grows -- reported in Amenta, N., Choi, S. and Rote, G., "Incremental //
1704// construction con {BRIO}," In Proceedings of 19th ACM Symposium on //
1705// Computational Geometry, 211-219, 2003. On the other hand, Liu and //
1706// Snoeyink showed that the point location can be made in constant time if //
1707// the points are pre-sorted so that the nearby points in space have nearby //
1708// indices, then adding the points in this order. They sorted the points //
1709// along the 3D Hilbert curve. //
1710// //
1711// The routine hilbert_sort3() sorts a set of 3D points along the 3D Hilbert //
1712// curve. It recursively splits a point set according to the Hilbert indices //
1713// mapped to the subboxes of the bounding box of the point set. //
1714// The Hilbert indices is calculated by Butz's algorithm in 1971. A nice //
1715// exposition of this algorithm can be found in the paper of Hamilton, C., //
1716// "Compact Hilbert Indices", Technical Report CS-2006-07, Computer Science, //
1717// Dalhousie University, 2006 (the Section 2). My implementation also refer- //
1718// enced Steven Witham's implementation of "Hilbert walk" (hopefully, it is //
1719// still available at: http://www.tiac.net/~sw/2008/10/Hilbert/). //
1720// //
1721// TetGen sorts the points using the method in the paper of Boissonnat,J.-D.,//
1722// Devillers, O. and Hornus, S. "Incremental Construction of the Delaunay //
1723// Triangulation and the Delaunay Graph in Medium Dimension," In Proceedings //
1724// of the 25th ACM Symposium on Computational Geometry, 2009. //
1725// It first randomly sorts the points into subgroups using the Biased Rand-//
1726// omized Insertion Ordering (BRIO) of Amenta et al 2003, then sorts the //
1727// points in each subgroup along the 3D Hilbert curve. Inserting points in //
1728// this order ensures a randomized "sprinkling" of the points over the //
1729// domain, while sorting of each subset ensures locality. //
1730// //
1731///////////////////////////////////////////////////////////////////////////////
1732
1733 void transfernodes();
1734
1735 // Point sorting.
1736 int transgc[8][3][8], tsb1mod3[8];
1737 void hilbert_init(int n);
1738 int hilbert_split(point* vertexarray, int arraysize, int gc0, int gc1,
1739 REAL, REAL, REAL, REAL, REAL, REAL);
1740 void hilbert_sort3(point* vertexarray, int arraysize, int e, int d,
1741 REAL, REAL, REAL, REAL, REAL, REAL, int depth);
1742 void brio_multiscale_sort(point*,int,int threshold,REAL ratio,int* depth);
1743
1744 // Point location.
1745 unsigned long randomnation(unsigned int choices);
1746 void randomsample(point searchpt, triface *searchtet);
1747 enum locateresult locate(point searchpt, triface *searchtet);
1748
1749 // Incremental flips.
1750 void flippush(badface*&, triface*);
1751 int incrementalflip(point newpt, int, flipconstraints *fc);
1752
1753 // Incremental Delaunay construction.
1754 void initialdelaunay(point pa, point pb, point pc, point pd);
1755 void incrementaldelaunay(clock_t&);
1756
1757///////////////////////////////////////////////////////////////////////////////
1758// //
1759// Surface triangulation //
1760// //
1761///////////////////////////////////////////////////////////////////////////////
1762
1763 void flipshpush(face*);
1764 void flip22(face*, int, int);
1765 void flip31(face*, int);
1766 long lawsonflip();
1767 int sinsertvertex(point newpt, face*, face*, int iloc, int bowywat, int);
1768 int sremovevertex(point delpt, face*, face*, int lawson);
1769
1770 enum locateresult slocate(point, face*, int, int, int);
1771 enum interresult sscoutsegment(face*, point);
1772 void scarveholes(int, REAL*);
1773 void triangulate(int, arraypool*, arraypool*, int, REAL*);
1774
1775 void unifysubfaces(face*, face*);
1776 void unifysegments();
1777 void mergefacets();
1778 void identifypscedges(point*);
1779 void meshsurface();
1780
1781 void interecursive(shellface** subfacearray, int arraysize, int axis,
1782 REAL, REAL, REAL, REAL, REAL, REAL, int* internum);
1783 void detectinterfaces();
1784
1785///////////////////////////////////////////////////////////////////////////////
1786// //
1787// Constrained Delaunay tetrahedralization //
1788// //
1789// A constrained Delaunay tetrahedralization (CDT) is a variation of a Dela- //
1790// unay tetrahedralization (DT) that is constrained to respect the boundary //
1791// of a 3D PLC (domain). In a CDT of a 3D PLC, every vertex or edge of the //
1792// PLC is also a vertex or an edge of the CDT, every polygon of the PLC is a //
1793// union of triangles of the CDT. A crucial difference between a CDT and a //
1794// DT is that triangles in the PLC's polygons are not required to be locally //
1795// Delaunay, which frees the CDT to better respect the PLC's polygons. CDTs //
1796// have optimal properties similar to those of DTs. //
1797// //
1798// Steiner Points and Steiner CDTs. It is known that even a simple 3D polyh- //
1799// edron may not have a tetrahedralization which only uses its own vertices. //
1800// Some extra points, so-called "Steiner points" are needed in order to form //
1801// a tetrahedralization of such polyhedron. It is true for tetrahedralizing //
1802// a 3D PLC as well. A Steiner CDT of a 3D PLC is a CDT containing Steiner //
1803// points. The CDT algorithms of TetGen in general create Steiner CDTs. //
1804// Almost all of the Steiner points are added in the edges of the PLC. They //
1805// guarantee the existence of a CDT of the modified PLC. //
1806// //
1807// The routine constraineddelaunay() starts from a DT of the vertices of a //
1808// PLC and creates a (Steiner) CDT of the PLC (including Steiner points). It //
1809// is constructed by two steps, (1) segment recovery and (2) facet (polygon) //
1810// recovery. Each step is accomplished by its own algorithm. //
1811// //
1812// The routine delaunizesegments() implements the segment recovery algorithm //
1813// of Si, H. and Gaertner, K. "Meshing Piecewise Linear Complexes by Constr- //
1814// ained Delaunay Tetrahedralizations," In Proceedings of the 14th Internat- //
1815// ional Meshing Roundtable, 147--163, 2005. It adds Steiner points into //
1816// non-Delaunay segments until all subsegments appear together in a DT. The //
1817// running time of this algorithm is proportional to the number of added //
1818// Steiner points. //
1819// //
1820// There are two incremental facet recovery algorithms: the cavity re-trian- //
1821// gulation algorithm of Si, H. and Gaertner, K. "3D Boundary Recovery by //
1822// Constrained Delaunay Tetrahedralization," International Journal for Numer-//
1823// ical Methods in Engineering, 85:1341-1364, 2011, and the flip algorithm //
1824// of Shewchuk, J. "Updating and Constructing Constrained Delaunay and //
1825// Constrained Regular Triangulations by Flips." In Proceedings of the 19th //
1826// ACM Symposium on Computational Geometry, 86-95, 2003. //
1827// //
1828// It is guaranteed in theory, no Steiner point is needed in both algorithms //
1829// However, a facet with non-coplanar vertices might cause the additions of //
1830// Steiner points. It is discussed in the paper of Si, H., and Shewchuk, J.,//
1831// "Incrementally Constructing and Updating Constrained Delaunay //
1832// Tetrahedralizations with Finite Precision Coordinates." In Proceedings of //
1833// the 21th International Meshing Roundtable, 2012. //
1834// //
1835// Our implementation of the facet recovery algorithms recover a "missing //
1836// region" at a time. Each missing region is a subset of connected interiors //
1837// of a polygon. The routine formcavity() creates the cavity of crossing //
1838// tetrahedra of the missing region. //
1839// //
1840// The cavity re-triangulation algorithm is implemented by three subroutines,//
1841// delaunizecavity(), fillcavity(), and carvecavity(). Since it may fail due //
1842// to non-coplanar vertices, the subroutine restorecavity() is used to rest- //
1843// ore the original cavity. //
1844// //
1845// The routine flipinsertfacet() implements the flip algorithm. The subrout- //
1846// ine flipcertify() is used to maintain the priority queue of flips. //
1847// //
1848// The routine refineregion() is called when the facet recovery algorithm //
1849// fail to recover a missing region. It inserts Steiner points to refine the //
1850// missing region. In order to avoid inserting Steiner points very close to //
1851// existing segments. The classical encroachment rules of the Delaunay //
1852// refinement algorithm are used to choose the Steiner points. //
1853// //
1854// The routine constrainedfacets() does the facet recovery by using either //
1855// the cavity re-triangulation algorithm (default) or the flip algorithm. It //
1856// results a CDT of the (modified) PLC (including Steiner points). //
1857// //
1858///////////////////////////////////////////////////////////////////////////////
1859
1860 void makesegmentendpointsmap();
1861
1862 enum interresult finddirection(triface* searchtet, point endpt);
1863 enum interresult scoutsegment(point, point, triface*, point*, arraypool*);
1864 int getsteinerptonsegment(face* seg, point refpt, point steinpt);
1865 void delaunizesegments();
1866
1867 enum interresult scoutsubface(face* searchsh, triface* searchtet);
1868 void formregion(face*, arraypool*, arraypool*, arraypool*);
1869 int scoutcrossedge(triface& crosstet, arraypool*, arraypool*);
1870 bool formcavity(triface*, arraypool*, arraypool*, arraypool*, arraypool*,
1871 arraypool*, arraypool*);
1872
1873 // Facet recovery by cavity re-triangulation [Si and Gaertner 2011].
1874 void delaunizecavity(arraypool*, arraypool*, arraypool*, arraypool*,
1875 arraypool*, arraypool*);
1876 bool fillcavity(arraypool*, arraypool*, arraypool*, arraypool*,
1877 arraypool*, arraypool*, triface* crossedge);
1878 void carvecavity(arraypool*, arraypool*, arraypool*);
1879 void restorecavity(arraypool*, arraypool*, arraypool*, arraypool*);
1880
1881 // Facet recovery by flips [Shewchuk 2003].
1882 void flipcertify(triface *chkface, badface **pqueue, point, point, point);
1883 void flipinsertfacet(arraypool*, arraypool*, arraypool*, arraypool*);
1884
1885 bool fillregion(arraypool* missingshs, arraypool*, arraypool* newshs);
1886
1887 int insertpoint_cdt(point, triface*, face*, face*, insertvertexflags*,
1888 arraypool*, arraypool*, arraypool*, arraypool*,
1889 arraypool*, arraypool*);
1890 void refineregion(face&, arraypool*, arraypool*, arraypool*, arraypool*,
1891 arraypool*, arraypool*);
1892
1893 void constrainedfacets();
1894
1895 void constraineddelaunay(clock_t&);
1896
1897///////////////////////////////////////////////////////////////////////////////
1898// //
1899// Constrained tetrahedralizations. //
1900// //
1901///////////////////////////////////////////////////////////////////////////////
1902
1903 int checkflipeligibility(int fliptype, point, point, point, point, point,
1904 int level, int edgepivot, flipconstraints* fc);
1905
1906 int removeedgebyflips(triface*, flipconstraints*);
1907 int removefacebyflips(triface*, flipconstraints*);
1908
1909 int recoveredgebyflips(point, point, triface*, int fullsearch);
1910 int add_steinerpt_in_schoenhardtpoly(triface*, int, int chkencflag);
1911 int add_steinerpt_in_segment(face*, int searchlevel);
1912 int addsteiner4recoversegment(face*, int);
1913 int recoversegments(arraypool*, int fullsearch, int steinerflag);
1914
1915 int recoverfacebyflips(point, point, point, face*, triface*);
1916 int recoversubfaces(arraypool*, int steinerflag);
1917
1918 int getvertexstar(int, point searchpt, arraypool*, arraypool*, arraypool*);
1919 int getedge(point, point, triface*);
1920 int reduceedgesatvertex(point startpt, arraypool* endptlist);
1921 int removevertexbyflips(point steinerpt);
1922
1923 int suppressbdrysteinerpoint(point steinerpt);
1924 int suppresssteinerpoints();
1925
1926 void recoverboundary(clock_t&);
1927
1928///////////////////////////////////////////////////////////////////////////////
1929// //
1930// Mesh reconstruction //
1931// //
1932///////////////////////////////////////////////////////////////////////////////
1933
1934 void carveholes();
1935
1936 void reconstructmesh();
1937
1938 int scoutpoint(point, triface*, int randflag);
1939 REAL getpointmeshsize(point, triface*, int iloc);
1940 void interpolatemeshsize();
1941
1942 void insertconstrainedpoints(point *insertarray, int arylen, int rejflag);
1943 void insertconstrainedpoints(tetgenio *addio);
1944
1945 void collectremovepoints(arraypool *remptlist);
1946 void meshcoarsening();
1947
1948///////////////////////////////////////////////////////////////////////////////
1949// //
1950// Mesh refinement //
1951// //
1952// The purpose of mesh refinement is to obtain a tetrahedral mesh with well- //
1953// -shaped tetrahedra and appropriate mesh size. It is necessary to insert //
1954// new Steiner points to achieve this property. The questions are (1) how to //
1955// choose the Steiner points? and (2) how to insert them? //
1956// //
1957// Delaunay refinement is a technique first developed by Chew [1989] and //
1958// Ruppert [1993, 1995] to generate quality triangular meshes in the plane. //
1959// It provides guarantee on the smallest angle of the triangles. Rupper's //
1960// algorithm guarantees that the mesh is size-optimal (to within a constant //
1961// factor) among all meshes with the same quality. //
1962// Shewchuk generalized Ruppert's algorithm into 3D in his PhD thesis //
1963// [Shewchuk 1997]. A short version of his algorithm appears in "Tetrahedral //
1964// Mesh Generation by Delaunay Refinement," In Proceedings of the 14th ACM //
1965// Symposium on Computational Geometry, 86-95, 1998. It guarantees that all //
1966// tetrahedra of the output mesh have a "radius-edge ratio" (equivalent to //
1967// the minimal face angle) bounded. However, it does not remove slivers, a //
1968// type of very flat tetrahedra which can have no small face angles but have //
1969// very small (and large) dihedral angles. Moreover, it may not terminate if //
1970// the input PLC contains "sharp features", e.g., two edges (or two facets) //
1971// meet at an acute angle (or dihedral angle). //
1972// //
1973// TetGen uses the basic Delaunay refinement scheme to insert Steiner points.//
1974// While it always maintains a constrained Delaunay mesh. The algorithm is //
1975// described in Si, H., "Adaptive Constrained Delaunay Mesh Generation," //
1976// International Journal for Numerical Methods in Engineering, 75:856-880. //
1977// This algorithm always terminates and sharp features are easily preserved. //
1978// The mesh has good quality (same as Shewchuk's Delaunay refinement algori- //
1979// thm) in the bulk of the mesh domain. Moreover, it supports the generation //
1980// of adaptive mesh according to a (isotropic) mesh sizing function. //
1981// //
1982///////////////////////////////////////////////////////////////////////////////
1983
1984 void makefacetverticesmap();
1985 int segsegadjacent(face *, face *);
1986 int segfacetadjacent(face *checkseg, face *checksh);
1987 int facetfacetadjacent(face *, face *);
1988
1989 int checkseg4encroach(point pa, point pb, point checkpt);
1990 int checkseg4split(face *chkseg, point&, int&);
1991 int splitsegment(face *splitseg, point encpt, REAL, point, point, int, int);
1992 void repairencsegs(int chkencflag);
1993
1994 void enqueuesubface(memorypool*, face*);
1995 int checkfac4encroach(point, point, point, point checkpt, REAL*, REAL*);
1996 int checkfac4split(face *chkfac, point& encpt, int& qflag, REAL *ccent);
1997 int splitsubface(face *splitfac, point, point, int qflag, REAL *ccent, int);
1998 void repairencfacs(int chkencflag);
1999
2000 void enqueuetetrahedron(triface*);
2001 int checktet4split(triface *chktet, int& qflag, REAL *ccent);
2002 int splittetrahedron(triface* splittet,int qflag,REAL *ccent, int);
2003 void repairbadtets(int chkencflag);
2004
2005 void delaunayrefinement();
2006
2007///////////////////////////////////////////////////////////////////////////////
2008// //
2009// Mesh optimization //
2010// //
2011///////////////////////////////////////////////////////////////////////////////
2012
2013 long lawsonflip3d(flipconstraints *fc);
2014 void recoverdelaunay();
2015
2016 int gettetrahedron(point, point, point, point, triface *);
2017 long improvequalitybyflips();
2018
2019 int smoothpoint(point smtpt, arraypool*, int ccw, optparameters *opm);
2020 long improvequalitybysmoothing(optparameters *opm);
2021
2022 int splitsliver(triface *, REAL, int);
2023 long removeslivers(int);
2024
2025 void optimizemesh();
2026
2027///////////////////////////////////////////////////////////////////////////////
2028// //
2029// Mesh check and statistics //
2030// //
2031///////////////////////////////////////////////////////////////////////////////
2032
2033 // Mesh validations.
2034 int checkmesh(int topoflag);
2035 int checkshells();
2036 int checksegments();
2037 int checkdelaunay();
2038 int checkregular(int);
2039 int checkconforming(int);
2040
2041 // Mesh statistics.
2042 void printfcomma(unsigned long n);
2043 void qualitystatistics();
2044 void memorystatistics();
2045 void statistics();
2046
2047///////////////////////////////////////////////////////////////////////////////
2048// //
2049// Mesh output //
2050// //
2051///////////////////////////////////////////////////////////////////////////////
2052
2053 void jettisonnodes();
2054 void highorder();
2055 void numberedges();
2056 void outnodes(tetgenio*);
2057 void outmetrics(tetgenio*);
2058 void outelements(tetgenio*);
2059 void outfaces(tetgenio*);
2060 void outhullfaces(tetgenio*);
2061 void outsubfaces(tetgenio*);
2062 void outedges(tetgenio*);
2063 void outsubsegments(tetgenio*);
2064 void outneighbors(tetgenio*);
2065 void outvoronoi(tetgenio*);
2066 void outsmesh(char*);
2067 void outmesh2medit(char*);
2068 void outmesh2vtk(char*);
2069
2070
2071
2072///////////////////////////////////////////////////////////////////////////////
2073// //
2074// Constructor & destructor //
2075// //
2076///////////////////////////////////////////////////////////////////////////////
2077
2078 tetgenmesh()
2079 {
2080 in = addin = NULL;
2081 b = NULL;
2082 bgm = NULL;
2083
2084 tetrahedrons = subfaces = subsegs = points = NULL;
2085 badtetrahedrons = badsubfacs = badsubsegs = NULL;
2086 tet2segpool = tet2subpool = NULL;
2087 flippool = NULL;
2088
2089 dummypoint = NULL;
2090 flipstack = NULL;
2091 unflipqueue = NULL;
2092
2093 cavetetlist = cavebdrylist = caveoldtetlist = NULL;
2094 cavetetshlist = cavetetseglist = cavetetvertlist = NULL;
2095 caveencshlist = caveencseglist = NULL;
2096 caveshlist = caveshbdlist = cavesegshlist = NULL;
2097
2098 subsegstack = subfacstack = subvertstack = NULL;
2099 encseglist = encshlist = NULL;
2100 idx2facetlist = NULL;
2101 facetverticeslist = NULL;
2102 segmentendpointslist = NULL;
2103
2104 highordertable = NULL;
2105
2106 numpointattrib = numelemattrib = 0;
2107 sizeoftensor = 0;
2108 pointmtrindex = 0;
2109 pointparamindex = 0;
2110 pointmarkindex = 0;
2111 point2simindex = 0;
2112 elemattribindex = 0;
2113 volumeboundindex = 0;
2114 shmarkindex = 0;
2115 areaboundindex = 0;
2116 checksubsegflag = 0;
2117 checksubfaceflag = 0;
2118 checkconstraints = 0;
2119 nonconvex = 0;
2120 autofliplinklevel = 1;
2121 useinsertradius = 0;
2122 samples = 0l;
2123 randomseed = 1l;
2124 minfaceang = minfacetdihed = PI;
2125 tetprism_vol_sum = 0.0;
2126 longest = 0.0;
2127 xmax = xmin = ymax = ymin = zmax = zmin = 0.0;
2128
2129 insegments = 0l;
2130 hullsize = 0l;
2131 meshedges = meshhulledges = 0l;
2132 steinerleft = -1;
2133 dupverts = 0l;
2134 unuverts = 0l;
2135 nonregularcount = 0l;
2136 st_segref_count = st_facref_count = st_volref_count = 0l;
2137 fillregioncount = cavitycount = cavityexpcount = 0l;
2138 flip14count = flip26count = flipn2ncount = 0l;
2139 flip23count = flip32count = flip44count = flip41count = 0l;
2140 flip22count = flip31count = 0l;
2141 totalworkmemory = 0l;
2142
2143
2144 } // tetgenmesh()
2145
2146 void freememory()
2147 {
2148 if (bgm != NULL) {
2149 delete bgm;
2150 }
2151
2152 if (points != (memorypool *) NULL) {
2153 delete points;
2154 delete [] dummypoint;
2155 }
2156
2157 if (tetrahedrons != (memorypool *) NULL) {
2158 delete tetrahedrons;
2159 }
2160
2161 if (subfaces != (memorypool *) NULL) {
2162 delete subfaces;
2163 delete subsegs;
2164 }
2165
2166 if (tet2segpool != NULL) {
2167 delete tet2segpool;
2168 delete tet2subpool;
2169 }
2170
2171 if (flippool != NULL) {
2172 delete flippool;
2173 delete unflipqueue;
2174 }
2175
2176 if (cavetetlist != NULL) {
2177 delete cavetetlist;
2178 delete cavebdrylist;
2179 delete caveoldtetlist;
2180 delete cavetetvertlist;
2181 }
2182
2183 if (caveshlist != NULL) {
2184 delete caveshlist;
2185 delete caveshbdlist;
2186 delete cavesegshlist;
2187 delete cavetetshlist;
2188 delete cavetetseglist;
2189 delete caveencshlist;
2190 delete caveencseglist;
2191 }
2192
2193 if (subsegstack != NULL) {
2194 delete subsegstack;
2195 delete subfacstack;
2196 delete subvertstack;
2197 }
2198
2199 if (idx2facetlist != NULL) {
2200 delete [] idx2facetlist;
2201 delete [] facetverticeslist;
2202 }
2203
2204 if (segmentendpointslist != NULL) {
2205 delete [] segmentendpointslist;
2206 }
2207
2208 if (highordertable != NULL) {
2209 delete [] highordertable;
2210 }
2211 }
2212
2213 ~tetgenmesh()
2214 {
2215 freememory();
2216 } // ~tetgenmesh()
2217
2218}; // End of class tetgenmesh.
2219
2220///////////////////////////////////////////////////////////////////////////////
2221// //
2222// tetrahedralize() Interface for using TetGen's library to generate //
2223// Delaunay tetrahedralizations, constrained Delaunay //
2224// tetrahedralizations, quality tetrahedral meshes. //
2225// //
2226// 'in' is an object of 'tetgenio' which contains a PLC you want to tetrahed-//
2227// ralize or a previously generated tetrahedral mesh you want to refine. It //
2228// must not be a NULL. 'out' is another object of 'tetgenio' for storing the //
2229// generated tetrahedral mesh. It can be a NULL. If so, the output will be //
2230// saved to file(s). If 'bgmin' != NULL, it contains a background mesh which //
2231// defines a mesh size function. //
2232// //
2233///////////////////////////////////////////////////////////////////////////////
2234
2235void tetrahedralize(tetgenbehavior *b, tetgenio *in, tetgenio *out,
2236 tetgenio *addin = NULL, tetgenio *bgmin = NULL);
2237
2238#ifdef TETLIBRARY
2239void tetrahedralize(char *switches, tetgenio *in, tetgenio *out,
2240 tetgenio *addin = NULL, tetgenio *bgmin = NULL);
2241#endif // #ifdef TETLIBRARY
2242
2243///////////////////////////////////////////////////////////////////////////////
2244// //
2245// terminatetetgen() Terminate TetGen with a given exit code. //
2246// //
2247///////////////////////////////////////////////////////////////////////////////
2248
2249inline void terminatetetgen(tetgenmesh *m, int x)
2250{
2251 // Release the allocated memory.
2252 if (m) {
2253 m->freememory();
2254 }
2255#ifdef TETLIBRARY
2256 switch (x) {
2257 case 1: // Out of memory.
2258 printf("Error: Out of memory.\n");
2259 break;
2260 case 2: // Encounter an internal error.
2261 printf("Please report this bug to Hang.Si@wias-berlin.de. Include\n");
2262 printf(" the message above, your input data set, and the exact\n");
2263 printf(" command line you used to run this program, thank you.\n");
2264 break;
2265 case 3:
2266 printf("A self-intersection was detected. Program stopped.\n");
2267 printf("Hint: use -d option to detect all self-intersections.\n");
2268 break;
2269 case 4:
2270 printf("A very small input feature size was detected. Program stopped.\n");
2271 printf("Hint: use -T option to set a smaller tolerance.\n");
2272 break;
2273 case 5:
2274 printf("Two very close input facets were detected. Program stopped.\n");
2275 printf("Hint: use -Y option to avoid adding Steiner points in boundary.\n");
2276 break;
2277 case 10:
2278 printf("An input error was detected. Program stopped.\n");
2279 break;
2280 } // switch (x)
2281 exit(x);
2282#endif // #ifdef TETLIBRARY
2283}
2284
2285///////////////////////////////////////////////////////////////////////////////
2286// //
2287// Primitives for tetrahedra //
2288// //
2289///////////////////////////////////////////////////////////////////////////////
2290
2291// encode() compress a handle into a single pointer. It relies on the
2292// assumption that all addresses of tetrahedra are aligned to sixteen-
2293// byte boundaries, so that the last four significant bits are zero.
2294
2295inline tetgenmesh::tetrahedron tetgenmesh::encode(triface& t) {
2296 return (tetrahedron) ((uintptr_t) (t).tet | (uintptr_t) (t).ver);
2297}
2298
2299inline tetgenmesh::tetrahedron tetgenmesh::encode2(tetrahedron* ptr, int ver) {
2300 return (tetrahedron) ((uintptr_t) (ptr) | (uintptr_t) (ver));
2301}
2302
2303// decode() converts a pointer to a handle. The version is extracted from
2304// the four least significant bits of the pointer.
2305
2306inline void tetgenmesh::decode(tetrahedron ptr, triface& t) {
2307 (t).ver = (int) ((uintptr_t) (ptr) & (uintptr_t) 15);
2308 (t).tet = (tetrahedron *) ((uintptr_t) (ptr) ^ (uintptr_t) (t).ver);
2309}
2310
2311// bond() connects two tetrahedra together. (t1,v1) and (t2,v2) must
2312// refer to the same face and the same edge.
2313
2314inline void tetgenmesh::bond(triface& t1, triface& t2) {
2315 t1.tet[t1.ver & 3] = encode2(t2.tet, bondtbl[t1.ver][t2.ver]);
2316 t2.tet[t2.ver & 3] = encode2(t1.tet, bondtbl[t2.ver][t1.ver]);
2317}
2318
2319
2320// dissolve() a bond (from one side).
2321
2322inline void tetgenmesh::dissolve(triface& t) {
2323 t.tet[t.ver & 3] = NULL;
2324}
2325
2326// enext() finds the next edge (counterclockwise) in the same face.
2327
2328inline void tetgenmesh::enext(triface& t1, triface& t2) {
2329 t2.tet = t1.tet;
2330 t2.ver = enexttbl[t1.ver];
2331}
2332
2333inline void tetgenmesh::enextself(triface& t) {
2334 t.ver = enexttbl[t.ver];
2335}
2336
2337// eprev() finds the next edge (clockwise) in the same face.
2338
2339inline void tetgenmesh::eprev(triface& t1, triface& t2) {
2340 t2.tet = t1.tet;
2341 t2.ver = eprevtbl[t1.ver];
2342}
2343
2344inline void tetgenmesh::eprevself(triface& t) {
2345 t.ver = eprevtbl[t.ver];
2346}
2347
2348// esym() finds the reversed edge. It is in the other face of the
2349// same tetrahedron.
2350
2351inline void tetgenmesh::esym(triface& t1, triface& t2) {
2352 (t2).tet = (t1).tet;
2353 (t2).ver = esymtbl[(t1).ver];
2354}
2355
2356inline void tetgenmesh::esymself(triface& t) {
2357 (t).ver = esymtbl[(t).ver];
2358}
2359
2360// enextesym() finds the reversed edge of the next edge. It is in the other
2361// face of the same tetrahedron. It is the combination esym() * enext().
2362
2363inline void tetgenmesh::enextesym(triface& t1, triface& t2) {
2364 t2.tet = t1.tet;
2365 t2.ver = enextesymtbl[t1.ver];
2366}
2367
2368inline void tetgenmesh::enextesymself(triface& t) {
2369 t.ver = enextesymtbl[t.ver];
2370}
2371
2372// eprevesym() finds the reversed edge of the previous edge.
2373
2374inline void tetgenmesh::eprevesym(triface& t1, triface& t2) {
2375 t2.tet = t1.tet;
2376 t2.ver = eprevesymtbl[t1.ver];
2377}
2378
2379inline void tetgenmesh::eprevesymself(triface& t) {
2380 t.ver = eprevesymtbl[t.ver];
2381}
2382
2383// eorgoppo() Finds the opposite face of the origin of the current edge.
2384// Return the opposite edge of the current edge.
2385
2386inline void tetgenmesh::eorgoppo(triface& t1, triface& t2) {
2387 t2.tet = t1.tet;
2388 t2.ver = eorgoppotbl[t1.ver];
2389}
2390
2391inline void tetgenmesh::eorgoppoself(triface& t) {
2392 t.ver = eorgoppotbl[t.ver];
2393}
2394
2395// edestoppo() Finds the opposite face of the destination of the current
2396// edge. Return the opposite edge of the current edge.
2397
2398inline void tetgenmesh::edestoppo(triface& t1, triface& t2) {
2399 t2.tet = t1.tet;
2400 t2.ver = edestoppotbl[t1.ver];
2401}
2402
2403inline void tetgenmesh::edestoppoself(triface& t) {
2404 t.ver = edestoppotbl[t.ver];
2405}
2406
2407// fsym() finds the adjacent tetrahedron at the same face and the same edge.
2408
2409inline void tetgenmesh::fsym(triface& t1, triface& t2) {
2410 decode((t1).tet[(t1).ver & 3], t2);
2411 t2.ver = fsymtbl[t1.ver][t2.ver];
2412}
2413
2414
2415#define fsymself(t) \
2416 t1ver = (t).ver; \
2417 decode((t).tet[(t).ver & 3], (t));\
2418 (t).ver = fsymtbl[t1ver][(t).ver]
2419
2420// fnext() finds the next face while rotating about an edge according to
2421// a right-hand rule. The face is in the adjacent tetrahedron. It is
2422// the combination: fsym() * esym().
2423
2424inline void tetgenmesh::fnext(triface& t1, triface& t2) {
2425 decode(t1.tet[facepivot1[t1.ver]], t2);
2426 t2.ver = facepivot2[t1.ver][t2.ver];
2427}
2428
2429
2430#define fnextself(t) \
2431 t1ver = (t).ver; \
2432 decode((t).tet[facepivot1[(t).ver]], (t)); \
2433 (t).ver = facepivot2[t1ver][(t).ver]
2434
2435
2436// The following primtives get or set the origin, destination, face apex,
2437// or face opposite of an ordered tetrahedron.
2438
2439inline tetgenmesh::point tetgenmesh::org(triface& t) {
2440 return (point) (t).tet[orgpivot[(t).ver]];
2441}
2442
2443inline tetgenmesh::point tetgenmesh:: dest(triface& t) {
2444 return (point) (t).tet[destpivot[(t).ver]];
2445}
2446
2447inline tetgenmesh::point tetgenmesh:: apex(triface& t) {
2448 return (point) (t).tet[apexpivot[(t).ver]];
2449}
2450
2451inline tetgenmesh::point tetgenmesh:: oppo(triface& t) {
2452 return (point) (t).tet[oppopivot[(t).ver]];
2453}
2454
2455inline void tetgenmesh:: setorg(triface& t, point p) {
2456 (t).tet[orgpivot[(t).ver]] = (tetrahedron) (p);
2457}
2458
2459inline void tetgenmesh:: setdest(triface& t, point p) {
2460 (t).tet[destpivot[(t).ver]] = (tetrahedron) (p);
2461}
2462
2463inline void tetgenmesh:: setapex(triface& t, point p) {
2464 (t).tet[apexpivot[(t).ver]] = (tetrahedron) (p);
2465}
2466
2467inline void tetgenmesh:: setoppo(triface& t, point p) {
2468 (t).tet[oppopivot[(t).ver]] = (tetrahedron) (p);
2469}
2470
2471#define setvertices(t, torg, tdest, tapex, toppo) \
2472 (t).tet[orgpivot[(t).ver]] = (tetrahedron) (torg);\
2473 (t).tet[destpivot[(t).ver]] = (tetrahedron) (tdest); \
2474 (t).tet[apexpivot[(t).ver]] = (tetrahedron) (tapex); \
2475 (t).tet[oppopivot[(t).ver]] = (tetrahedron) (toppo)
2476
2477// Check or set a tetrahedron's attributes.
2478
2479inline REAL tetgenmesh::elemattribute(tetrahedron* ptr, int attnum) {
2480 return ((REAL *) (ptr))[elemattribindex + attnum];
2481}
2482
2483inline void tetgenmesh::setelemattribute(tetrahedron* ptr, int attnum,
2484 REAL value) {
2485 ((REAL *) (ptr))[elemattribindex + attnum] = value;
2486}
2487
2488// Check or set a tetrahedron's maximum volume bound.
2489
2490inline REAL tetgenmesh::volumebound(tetrahedron* ptr) {
2491 return ((REAL *) (ptr))[volumeboundindex];
2492}
2493
2494inline void tetgenmesh::setvolumebound(tetrahedron* ptr, REAL value) {
2495 ((REAL *) (ptr))[volumeboundindex] = value;
2496}
2497
2498// Get or set a tetrahedron's index (only used for output).
2499// These two routines use the reserved slot ptr[10].
2500
2501inline int tetgenmesh::elemindex(tetrahedron* ptr) {
2502 int *iptr = (int *) &(ptr[10]);
2503 return iptr[0];
2504}
2505
2506inline void tetgenmesh::setelemindex(tetrahedron* ptr, int value) {
2507 int *iptr = (int *) &(ptr[10]);
2508 iptr[0] = value;
2509}
2510
2511// Get or set a tetrahedron's marker.
2512// Set 'value = 0' cleans all the face/edge flags.
2513
2514inline int tetgenmesh::elemmarker(tetrahedron* ptr) {
2515 return ((int *) (ptr))[elemmarkerindex];
2516}
2517
2518inline void tetgenmesh::setelemmarker(tetrahedron* ptr, int value) {
2519 ((int *) (ptr))[elemmarkerindex] = value;
2520}
2521
2522// infect(), infected(), uninfect() -- primitives to flag or unflag a
2523// tetrahedron. The last bit of the element marker is flagged (1)
2524// or unflagged (0).
2525
2526inline void tetgenmesh::infect(triface& t) {
2527 ((int *) (t.tet))[elemmarkerindex] |= 1;
2528}
2529
2530inline void tetgenmesh::uninfect(triface& t) {
2531 ((int *) (t.tet))[elemmarkerindex] &= ~1;
2532}
2533
2534inline bool tetgenmesh::infected(triface& t) {
2535 return (((int *) (t.tet))[elemmarkerindex] & 1) != 0;
2536}
2537
2538// marktest(), marktested(), unmarktest() -- primitives to flag or unflag a
2539// tetrahedron. Use the second lowerest bit of the element marker.
2540
2541inline void tetgenmesh::marktest(triface& t) {
2542 ((int *) (t.tet))[elemmarkerindex] |= 2;
2543}
2544
2545inline void tetgenmesh::unmarktest(triface& t) {
2546 ((int *) (t.tet))[elemmarkerindex] &= ~2;
2547}
2548
2549inline bool tetgenmesh::marktested(triface& t) {
2550 return (((int *) (t.tet))[elemmarkerindex] & 2) != 0;
2551}
2552
2553// markface(), unmarkface(), facemarked() -- primitives to flag or unflag a
2554// face of a tetrahedron. From the last 3rd to 6th bits are used for
2555// face markers, e.g., the last third bit corresponds to loc = 0.
2556
2557inline void tetgenmesh::markface(triface& t) {
2558 ((int *) (t.tet))[elemmarkerindex] |= (4 << (t.ver & 3));
2559}
2560
2561inline void tetgenmesh::unmarkface(triface& t) {
2562 ((int *) (t.tet))[elemmarkerindex] &= ~(4 << (t.ver & 3));
2563}
2564
2565inline bool tetgenmesh::facemarked(triface& t) {
2566 return (((int *) (t.tet))[elemmarkerindex] & (4 << (t.ver & 3))) != 0;
2567}
2568
2569// markedge(), unmarkedge(), edgemarked() -- primitives to flag or unflag an
2570// edge of a tetrahedron. From the last 7th to 12th bits are used for
2571// edge markers, e.g., the last 7th bit corresponds to the 0th edge, etc.
2572// Remark: The last 7th bit is marked by 2^6 = 64.
2573
2574inline void tetgenmesh::markedge(triface& t) {
2575 ((int *) (t.tet))[elemmarkerindex] |= (int) (64 << ver2edge[(t).ver]);
2576}
2577
2578inline void tetgenmesh::unmarkedge(triface& t) {
2579 ((int *) (t.tet))[elemmarkerindex] &= ~(int) (64 << ver2edge[(t).ver]);
2580}
2581
2582inline bool tetgenmesh::edgemarked(triface& t) {
2583 return (((int *) (t.tet))[elemmarkerindex] &
2584 (int) (64 << ver2edge[(t).ver])) != 0;
2585}
2586
2587// marktest2(), unmarktest2(), marktest2ed() -- primitives to flag and unflag
2588// a tetrahedron. The 13th bit (2^12 = 4096) is used for this flag.
2589
2590inline void tetgenmesh::marktest2(triface& t) {
2591 ((int *) (t.tet))[elemmarkerindex] |= (int) (4096);
2592}
2593
2594inline void tetgenmesh::unmarktest2(triface& t) {
2595 ((int *) (t.tet))[elemmarkerindex] &= ~(int) (4096);
2596}
2597
2598inline bool tetgenmesh::marktest2ed(triface& t) {
2599 return (((int *) (t.tet))[elemmarkerindex] & (int) (4096)) != 0;
2600}
2601
2602// elemcounter(), setelemcounter() -- primitives to read or ser a (small)
2603// integer counter in this tet. It is saved from the 16th bit. On 32 bit
2604// system, the range of the counter is [0, 2^15 = 32768].
2605
2606inline int tetgenmesh::elemcounter(triface& t) {
2607 return (((int *) (t.tet))[elemmarkerindex]) >> 16;
2608}
2609
2610inline void tetgenmesh::setelemcounter(triface& t, int value) {
2611 int c = ((int *) (t.tet))[elemmarkerindex];
2612 // Clear the old counter while keep the other flags.
2613 c &= 65535; // sum_{i=0^15} 2^i
2614 c |= (value << 16);
2615 ((int *) (t.tet))[elemmarkerindex] = c;
2616}
2617
2618inline void tetgenmesh::increaseelemcounter(triface& t) {
2619 int c = elemcounter(t);
2620 setelemcounter(t, c + 1);
2621}
2622
2623inline void tetgenmesh::decreaseelemcounter(triface& t) {
2624 int c = elemcounter(t);
2625 setelemcounter(t, c - 1);
2626}
2627
2628// ishulltet() tests if t is a hull tetrahedron.
2629
2630inline bool tetgenmesh::ishulltet(triface& t) {
2631 return (point) (t).tet[7] == dummypoint;
2632}
2633
2634// isdeadtet() tests if t is a tetrahedron is dead.
2635
2636inline bool tetgenmesh::isdeadtet(triface& t) {
2637 return ((t.tet == NULL) || (t.tet[4] == NULL));
2638}
2639
2640///////////////////////////////////////////////////////////////////////////////
2641// //
2642// Primitives for subfaces and subsegments //
2643// //
2644///////////////////////////////////////////////////////////////////////////////
2645
2646// Each subface contains three pointers to its neighboring subfaces, with
2647// edge versions. To save memory, both information are kept in a single
2648// pointer. To make this possible, all subfaces are aligned to eight-byte
2649// boundaries, so that the last three bits of each pointer are zeros. An
2650// edge version (in the range 0 to 5) is compressed into the last three
2651// bits of each pointer by 'sencode()'. 'sdecode()' decodes a pointer,
2652// extracting an edge version and a pointer to the beginning of a subface.
2653
2654inline void tetgenmesh::sdecode(shellface sptr, face& s) {
2655 s.shver = (int) ((uintptr_t) (sptr) & (uintptr_t) 7);
2656 s.sh = (shellface *) ((uintptr_t) (sptr) ^ (uintptr_t) (s.shver));
2657}
2658
2659inline tetgenmesh::shellface tetgenmesh::sencode(face& s) {
2660 return (shellface) ((uintptr_t) s.sh | (uintptr_t) s.shver);
2661}
2662
2663inline tetgenmesh::shellface tetgenmesh::sencode2(shellface *sh, int shver) {
2664 return (shellface) ((uintptr_t) sh | (uintptr_t) shver);
2665}
2666
2667// sbond() bonds two subfaces (s1) and (s2) together. s1 and s2 must refer
2668// to the same edge. No requirement is needed on their orientations.
2669
2670inline void tetgenmesh::sbond(face& s1, face& s2)
2671{
2672 s1.sh[s1.shver >> 1] = sencode(s2);
2673 s2.sh[s2.shver >> 1] = sencode(s1);
2674}
2675
2676// sbond1() bonds s1 <== s2, i.e., after bonding, s1 is pointing to s2,
2677// but s2 is not pointing to s1. s1 and s2 must refer to the same edge.
2678// No requirement is needed on their orientations.
2679
2680inline void tetgenmesh::sbond1(face& s1, face& s2)
2681{
2682 s1.sh[s1.shver >> 1] = sencode(s2);
2683}
2684
2685// Dissolve a subface bond (from one side). Note that the other subface
2686// will still think it's connected to this subface.
2687
2688inline void tetgenmesh::sdissolve(face& s)
2689{
2690 s.sh[s.shver >> 1] = NULL;
2691}
2692
2693// spivot() finds the adjacent subface (s2) for a given subface (s1).
2694// s1 and s2 share at the same edge.
2695
2696inline void tetgenmesh::spivot(face& s1, face& s2)
2697{
2698 shellface sptr = s1.sh[s1.shver >> 1];
2699 sdecode(sptr, s2);
2700}
2701
2702inline void tetgenmesh::spivotself(face& s)
2703{
2704 shellface sptr = s.sh[s.shver >> 1];
2705 sdecode(sptr, s);
2706}
2707
2708// These primitives determine or set the origin, destination, or apex
2709// of a subface with respect to the edge version.
2710
2711inline tetgenmesh::point tetgenmesh::sorg(face& s)
2712{
2713 return (point) s.sh[sorgpivot[s.shver]];
2714}
2715
2716inline tetgenmesh::point tetgenmesh::sdest(face& s)
2717{
2718 return (point) s.sh[sdestpivot[s.shver]];
2719}
2720
2721inline tetgenmesh::point tetgenmesh::sapex(face& s)
2722{
2723 return (point) s.sh[sapexpivot[s.shver]];
2724}
2725
2726inline void tetgenmesh::setsorg(face& s, point pointptr)
2727{
2728 s.sh[sorgpivot[s.shver]] = (shellface) pointptr;
2729}
2730
2731inline void tetgenmesh::setsdest(face& s, point pointptr)
2732{
2733 s.sh[sdestpivot[s.shver]] = (shellface) pointptr;
2734}
2735
2736inline void tetgenmesh::setsapex(face& s, point pointptr)
2737{
2738 s.sh[sapexpivot[s.shver]] = (shellface) pointptr;
2739}
2740
2741#define setshvertices(s, pa, pb, pc)\
2742 setsorg(s, pa);\
2743 setsdest(s, pb);\
2744 setsapex(s, pc)
2745
2746// sesym() reserves the direction of the lead edge.
2747
2748inline void tetgenmesh::sesym(face& s1, face& s2)
2749{
2750 s2.sh = s1.sh;
2751 s2.shver = (s1.shver ^ 1); // Inverse the last bit.
2752}
2753
2754inline void tetgenmesh::sesymself(face& s)
2755{
2756 s.shver ^= 1;
2757}
2758
2759// senext() finds the next edge (counterclockwise) in the same orientation
2760// of this face.
2761
2762inline void tetgenmesh::senext(face& s1, face& s2)
2763{
2764 s2.sh = s1.sh;
2765 s2.shver = snextpivot[s1.shver];
2766}
2767
2768inline void tetgenmesh::senextself(face& s)
2769{
2770 s.shver = snextpivot[s.shver];
2771}
2772
2773inline void tetgenmesh::senext2(face& s1, face& s2)
2774{
2775 s2.sh = s1.sh;
2776 s2.shver = snextpivot[snextpivot[s1.shver]];
2777}
2778
2779inline void tetgenmesh::senext2self(face& s)
2780{
2781 s.shver = snextpivot[snextpivot[s.shver]];
2782}
2783
2784
2785// Check or set a subface's maximum area bound.
2786
2787inline REAL tetgenmesh::areabound(face& s)
2788{
2789 return ((REAL *) (s.sh))[areaboundindex];
2790}
2791
2792inline void tetgenmesh::setareabound(face& s, REAL value)
2793{
2794 ((REAL *) (s.sh))[areaboundindex] = value;
2795}
2796
2797// These two primitives read or set a shell marker. Shell markers are used
2798// to hold user boundary information.
2799
2800inline int tetgenmesh::shellmark(face& s)
2801{
2802 return ((int *) (s.sh))[shmarkindex];
2803}
2804
2805inline void tetgenmesh::setshellmark(face& s, int value)
2806{
2807 ((int *) (s.sh))[shmarkindex] = value;
2808}
2809
2810
2811
2812// sinfect(), sinfected(), suninfect() -- primitives to flag or unflag a
2813// subface. The last bit of ((int *) ((s).sh))[shmarkindex+1] is flagged.
2814
2815inline void tetgenmesh::sinfect(face& s)
2816{
2817 ((int *) ((s).sh))[shmarkindex+1] =
2818 (((int *) ((s).sh))[shmarkindex+1] | (int) 1);
2819}
2820
2821inline void tetgenmesh::suninfect(face& s)
2822{
2823 ((int *) ((s).sh))[shmarkindex+1] =
2824 (((int *) ((s).sh))[shmarkindex+1] & ~(int) 1);
2825}
2826
2827// Test a subface for viral infection.
2828
2829inline bool tetgenmesh::sinfected(face& s)
2830{
2831 return (((int *) ((s).sh))[shmarkindex+1] & (int) 1) != 0;
2832}
2833
2834// smarktest(), smarktested(), sunmarktest() -- primitives to flag or unflag
2835// a subface. The last 2nd bit of the integer is flagged.
2836
2837inline void tetgenmesh::smarktest(face& s)
2838{
2839 ((int *) ((s).sh))[shmarkindex+1] =
2840 (((int *)((s).sh))[shmarkindex+1] | (int) 2);
2841}
2842
2843inline void tetgenmesh::sunmarktest(face& s)
2844{
2845 ((int *) ((s).sh))[shmarkindex+1] =
2846 (((int *)((s).sh))[shmarkindex+1] & ~(int)2);
2847}
2848
2849inline bool tetgenmesh::smarktested(face& s)
2850{
2851 return ((((int *) ((s).sh))[shmarkindex+1] & (int) 2) != 0);
2852}
2853
2854// smarktest2(), smarktest2ed(), sunmarktest2() -- primitives to flag or
2855// unflag a subface. The last 3rd bit of the integer is flagged.
2856
2857inline void tetgenmesh::smarktest2(face& s)
2858{
2859 ((int *) ((s).sh))[shmarkindex+1] =
2860 (((int *)((s).sh))[shmarkindex+1] | (int) 4);
2861}
2862
2863inline void tetgenmesh::sunmarktest2(face& s)
2864{
2865 ((int *) ((s).sh))[shmarkindex+1] =
2866 (((int *)((s).sh))[shmarkindex+1] & ~(int)4);
2867}
2868
2869inline bool tetgenmesh::smarktest2ed(face& s)
2870{
2871 return ((((int *) ((s).sh))[shmarkindex+1] & (int) 4) != 0);
2872}
2873
2874// The last 4th bit of ((int *) ((s).sh))[shmarkindex+1] is flagged.
2875
2876inline void tetgenmesh::smarktest3(face& s)
2877{
2878 ((int *) ((s).sh))[shmarkindex+1] =
2879 (((int *)((s).sh))[shmarkindex+1] | (int) 8);
2880}
2881
2882inline void tetgenmesh::sunmarktest3(face& s)
2883{
2884 ((int *) ((s).sh))[shmarkindex+1] =
2885 (((int *)((s).sh))[shmarkindex+1] & ~(int)8);
2886}
2887
2888inline bool tetgenmesh::smarktest3ed(face& s)
2889{
2890 return ((((int *) ((s).sh))[shmarkindex+1] & (int) 8) != 0);
2891}
2892
2893
2894// Each facet has a unique index (automatically indexed). Starting from '0'.
2895// We save this index in the same field of the shell type.
2896
2897inline void tetgenmesh::setfacetindex(face& s, int value)
2898{
2899 ((int *) (s.sh))[shmarkindex + 2] = value;
2900}
2901
2902inline int tetgenmesh::getfacetindex(face& s)
2903{
2904 return ((int *) (s.sh))[shmarkindex + 2];
2905}
2906
2907///////////////////////////////////////////////////////////////////////////////
2908// //
2909// Primitives for interacting between tetrahedra and subfaces //
2910// //
2911///////////////////////////////////////////////////////////////////////////////
2912
2913// tsbond() bond a tetrahedron (t) and a subface (s) together.
2914// Note that t and s must be the same face and the same edge. Moreover,
2915// t and s have the same orientation.
2916// Since the edge number in t and in s can be any number in {0,1,2}. We bond
2917// the edge in s which corresponds to t's 0th edge, and vice versa.
2918
2919inline void tetgenmesh::tsbond(triface& t, face& s)
2920{
2921 if ((t).tet[9] == NULL) {
2922 // Allocate space for this tet.
2923 (t).tet[9] = (tetrahedron) tet2subpool->alloc();
2924 // Initialize.
2925 for (int i = 0; i < 4; i++) {
2926 ((shellface *) (t).tet[9])[i] = NULL;
2927 }
2928 }
2929 // Bond t <== s.
2930 ((shellface *) (t).tet[9])[(t).ver & 3] =
2931 sencode2((s).sh, tsbondtbl[t.ver][s.shver]);
2932 // Bond s <== t.
2933 s.sh[9 + ((s).shver & 1)] =
2934 (shellface) encode2((t).tet, stbondtbl[t.ver][s.shver]);
2935}
2936
2937// tspivot() finds a subface (s) abutting on the given tetrahdera (t).
2938// Return s.sh = NULL if there is no subface at t. Otherwise, return
2939// the subface s, and s and t must be at the same edge wth the same
2940// orientation.
2941
2942inline void tetgenmesh::tspivot(triface& t, face& s)
2943{
2944 if ((t).tet[9] == NULL) {
2945 (s).sh = NULL;
2946 return;
2947 }
2948 // Get the attached subface s.
2949 sdecode(((shellface *) (t).tet[9])[(t).ver & 3], (s));
2950 (s).shver = tspivottbl[t.ver][s.shver];
2951}
2952
2953// Quickly check if the handle (t, v) is a subface.
2954#define issubface(t) \
2955 ((t).tet[9] && ((t).tet[9])[(t).ver & 3])
2956
2957// stpivot() finds a tetrahedron (t) abutting a given subface (s).
2958// Return the t (if it exists) with the same edge and the same
2959// orientation of s.
2960
2961inline void tetgenmesh::stpivot(face& s, triface& t)
2962{
2963 decode((tetrahedron) s.sh[9 + (s.shver & 1)], t);
2964 if ((t).tet == NULL) {
2965 return;
2966 }
2967 (t).ver = stpivottbl[t.ver][s.shver];
2968}
2969
2970// Quickly check if this subface is attached to a tetrahedron.
2971
2972#define isshtet(s) \
2973 ((s).sh[9 + ((s).shver & 1)])
2974
2975// tsdissolve() dissolve a bond (from the tetrahedron side).
2976
2977inline void tetgenmesh::tsdissolve(triface& t)
2978{
2979 if ((t).tet[9] != NULL) {
2980 ((shellface *) (t).tet[9])[(t).ver & 3] = NULL;
2981 }
2982}
2983
2984// stdissolve() dissolve a bond (from the subface side).
2985
2986inline void tetgenmesh::stdissolve(face& s)
2987{
2988 (s).sh[9] = NULL;
2989 (s).sh[10] = NULL;
2990}
2991
2992///////////////////////////////////////////////////////////////////////////////
2993// //
2994// Primitives for interacting between subfaces and segments //
2995// //
2996///////////////////////////////////////////////////////////////////////////////
2997
2998// ssbond() bond a subface to a subsegment.
2999
3000inline void tetgenmesh::ssbond(face& s, face& edge)
3001{
3002 s.sh[6 + (s.shver >> 1)] = sencode(edge);
3003 edge.sh[0] = sencode(s);
3004}
3005
3006inline void tetgenmesh::ssbond1(face& s, face& edge)
3007{
3008 s.sh[6 + (s.shver >> 1)] = sencode(edge);
3009 //edge.sh[0] = sencode(s);
3010}
3011
3012// ssdisolve() dissolve a bond (from the subface side)
3013
3014inline void tetgenmesh::ssdissolve(face& s)
3015{
3016 s.sh[6 + (s.shver >> 1)] = NULL;
3017}
3018
3019// sspivot() finds a subsegment abutting a subface.
3020
3021inline void tetgenmesh::sspivot(face& s, face& edge)
3022{
3023 sdecode((shellface) s.sh[6 + (s.shver >> 1)], edge);
3024}
3025
3026// Quickly check if the edge is a subsegment.
3027
3028#define isshsubseg(s) \
3029 ((s).sh[6 + ((s).shver >> 1)])
3030
3031///////////////////////////////////////////////////////////////////////////////
3032// //
3033// Primitives for interacting between tetrahedra and segments //
3034// //
3035///////////////////////////////////////////////////////////////////////////////
3036
3037inline void tetgenmesh::tssbond1(triface& t, face& s)
3038{
3039 if ((t).tet[8] == NULL) {
3040 // Allocate space for this tet.
3041 (t).tet[8] = (tetrahedron) tet2segpool->alloc();
3042 // Initialization.
3043 for (int i = 0; i < 6; i++) {
3044 ((shellface *) (t).tet[8])[i] = NULL;
3045 }
3046 }
3047 ((shellface *) (t).tet[8])[ver2edge[(t).ver]] = sencode((s));
3048}
3049
3050inline void tetgenmesh::sstbond1(face& s, triface& t)
3051{
3052 ((tetrahedron *) (s).sh)[9] = encode(t);
3053}
3054
3055inline void tetgenmesh::tssdissolve1(triface& t)
3056{
3057 if ((t).tet[8] != NULL) {
3058 ((shellface *) (t).tet[8])[ver2edge[(t).ver]] = NULL;
3059 }
3060}
3061
3062inline void tetgenmesh::sstdissolve1(face& s)
3063{
3064 ((tetrahedron *) (s).sh)[9] = NULL;
3065}
3066
3067inline void tetgenmesh::tsspivot1(triface& t, face& s)
3068{
3069 if ((t).tet[8] != NULL) {
3070 sdecode(((shellface *) (t).tet[8])[ver2edge[(t).ver]], s);
3071 } else {
3072 (s).sh = NULL;
3073 }
3074}
3075
3076// Quickly check whether 't' is a segment or not.
3077
3078#define issubseg(t) \
3079 ((t).tet[8] && ((t).tet[8])[ver2edge[(t).ver]])
3080
3081inline void tetgenmesh::sstpivot1(face& s, triface& t)
3082{
3083 decode((tetrahedron) s.sh[9], t);
3084}
3085
3086///////////////////////////////////////////////////////////////////////////////
3087// //
3088// Primitives for points //
3089// //
3090///////////////////////////////////////////////////////////////////////////////
3091
3092inline int tetgenmesh::pointmark(point pt) {
3093 return ((int *) (pt))[pointmarkindex];
3094}
3095
3096inline void tetgenmesh::setpointmark(point pt, int value) {
3097 ((int *) (pt))[pointmarkindex] = value;
3098}
3099
3100
3101// These two primitives set and read the type of the point.
3102
3103inline enum tetgenmesh::verttype tetgenmesh::pointtype(point pt) {
3104 return (enum verttype) (((int *) (pt))[pointmarkindex + 1] >> (int) 8);
3105}
3106
3107inline void tetgenmesh::setpointtype(point pt, enum verttype value) {
3108 ((int *) (pt))[pointmarkindex + 1] =
3109 ((int) value << 8) + (((int *) (pt))[pointmarkindex + 1] & (int) 255);
3110}
3111
3112// Read and set the geometry tag of the point (used by -s option).
3113
3114inline int tetgenmesh::pointgeomtag(point pt) {
3115 return ((int *) (pt))[pointmarkindex + 2];
3116}
3117
3118inline void tetgenmesh::setpointgeomtag(point pt, int value) {
3119 ((int *) (pt))[pointmarkindex + 2] = value;
3120}
3121
3122// Read and set the u,v coordinates of the point (used by -s option).
3123
3124inline REAL tetgenmesh::pointgeomuv(point pt, int i) {
3125 return pt[pointparamindex + i];
3126}
3127
3128inline void tetgenmesh::setpointgeomuv(point pt, int i, REAL value) {
3129 pt[pointparamindex + i] = value;
3130}
3131
3132// pinfect(), puninfect(), pinfected() -- primitives to flag or unflag
3133// a point. The last bit of the integer '[pointindex+1]' is flagged.
3134
3135inline void tetgenmesh::pinfect(point pt) {
3136 ((int *) (pt))[pointmarkindex + 1] |= (int) 1;
3137}
3138
3139inline void tetgenmesh::puninfect(point pt) {
3140 ((int *) (pt))[pointmarkindex + 1] &= ~(int) 1;
3141}
3142
3143inline bool tetgenmesh::pinfected(point pt) {
3144 return (((int *) (pt))[pointmarkindex + 1] & (int) 1) != 0;
3145}
3146
3147// pmarktest(), punmarktest(), pmarktested() -- more primitives to
3148// flag or unflag a point.
3149
3150inline void tetgenmesh::pmarktest(point pt) {
3151 ((int *) (pt))[pointmarkindex + 1] |= (int) 2;
3152}
3153
3154inline void tetgenmesh::punmarktest(point pt) {
3155 ((int *) (pt))[pointmarkindex + 1] &= ~(int) 2;
3156}
3157
3158inline bool tetgenmesh::pmarktested(point pt) {
3159 return (((int *) (pt))[pointmarkindex + 1] & (int) 2) != 0;
3160}
3161
3162inline void tetgenmesh::pmarktest2(point pt) {
3163 ((int *) (pt))[pointmarkindex + 1] |= (int) 4;
3164}
3165
3166inline void tetgenmesh::punmarktest2(point pt) {
3167 ((int *) (pt))[pointmarkindex + 1] &= ~(int) 4;
3168}
3169
3170inline bool tetgenmesh::pmarktest2ed(point pt) {
3171 return (((int *) (pt))[pointmarkindex + 1] & (int) 4) != 0;
3172}
3173
3174inline void tetgenmesh::pmarktest3(point pt) {
3175 ((int *) (pt))[pointmarkindex + 1] |= (int) 8;
3176}
3177
3178inline void tetgenmesh::punmarktest3(point pt) {
3179 ((int *) (pt))[pointmarkindex + 1] &= ~(int) 8;
3180}
3181
3182inline bool tetgenmesh::pmarktest3ed(point pt) {
3183 return (((int *) (pt))[pointmarkindex + 1] & (int) 8) != 0;
3184}
3185
3186// These following primitives set and read a pointer to a tetrahedron
3187// a subface/subsegment, a point, or a tet of background mesh.
3188
3189inline tetgenmesh::tetrahedron tetgenmesh::point2tet(point pt) {
3190 return ((tetrahedron *) (pt))[point2simindex];
3191}
3192
3193inline void tetgenmesh::setpoint2tet(point pt, tetrahedron value) {
3194 ((tetrahedron *) (pt))[point2simindex] = value;
3195}
3196
3197inline tetgenmesh::point tetgenmesh::point2ppt(point pt) {
3198 return (point) ((tetrahedron *) (pt))[point2simindex + 1];
3199}
3200
3201inline void tetgenmesh::setpoint2ppt(point pt, point value) {
3202 ((tetrahedron *) (pt))[point2simindex + 1] = (tetrahedron) value;
3203}
3204
3205inline tetgenmesh::shellface tetgenmesh::point2sh(point pt) {
3206 return (shellface) ((tetrahedron *) (pt))[point2simindex + 2];
3207}
3208
3209inline void tetgenmesh::setpoint2sh(point pt, shellface value) {
3210 ((tetrahedron *) (pt))[point2simindex + 2] = (tetrahedron) value;
3211}
3212
3213
3214inline tetgenmesh::tetrahedron tetgenmesh::point2bgmtet(point pt) {
3215 return ((tetrahedron *) (pt))[point2simindex + 3];
3216}
3217
3218inline void tetgenmesh::setpoint2bgmtet(point pt, tetrahedron value) {
3219 ((tetrahedron *) (pt))[point2simindex + 3] = value;
3220}
3221
3222
3223// The primitives for saving and getting the insertion radius.
3224inline void tetgenmesh::setpointinsradius(point pt, REAL value)
3225{
3226 pt[pointmtrindex + sizeoftensor - 1] = value;
3227}
3228
3229inline REAL tetgenmesh::getpointinsradius(point pt)
3230{
3231 return pt[pointmtrindex + sizeoftensor - 1];
3232}
3233
3234// point2tetorg() Get the tetrahedron whose origin is the point.
3235
3236inline void tetgenmesh::point2tetorg(point pa, triface& searchtet)
3237{
3238 decode(point2tet(pa), searchtet);
3239 if ((point) searchtet.tet[4] == pa) {
3240 searchtet.ver = 11;
3241 } else if ((point) searchtet.tet[5] == pa) {
3242 searchtet.ver = 3;
3243 } else if ((point) searchtet.tet[6] == pa) {
3244 searchtet.ver = 7;
3245 } else {
3246 assert((point) searchtet.tet[7] == pa); // SELF_CHECK
3247 searchtet.ver = 0;
3248 }
3249}
3250
3251// point2shorg() Get the subface/segment whose origin is the point.
3252
3253inline void tetgenmesh::point2shorg(point pa, face& searchsh)
3254{
3255 sdecode(point2sh(pa), searchsh);
3256 if ((point) searchsh.sh[3] == pa) {
3257 searchsh.shver = 0;
3258 } else if ((point) searchsh.sh[4] == pa) {
3259 searchsh.shver = (searchsh.sh[5] != NULL ? 2 : 1);
3260 } else {
3261 assert((point) searchsh.sh[5] == pa); // SELF_CHECK
3262 searchsh.shver = 4;
3263 }
3264}
3265
3266// farsorg() Return the origin of the subsegment.
3267// farsdest() Return the destination of the subsegment.
3268
3269inline tetgenmesh::point tetgenmesh::farsorg(face& s)
3270{
3271 face travesh, neighsh;
3272
3273 travesh = s;
3274 while (1) {
3275 senext2(travesh, neighsh);
3276 spivotself(neighsh);
3277 if (neighsh.sh == NULL) break;
3278 if (sorg(neighsh) != sorg(travesh)) sesymself(neighsh);
3279 senext2(neighsh, travesh);
3280 }
3281 return sorg(travesh);
3282}
3283
3284inline tetgenmesh::point tetgenmesh::farsdest(face& s)
3285{
3286 face travesh, neighsh;
3287
3288 travesh = s;
3289 while (1) {
3290 senext(travesh, neighsh);
3291 spivotself(neighsh);
3292 if (neighsh.sh == NULL) break;
3293 if (sdest(neighsh) != sdest(travesh)) sesymself(neighsh);
3294 senext(neighsh, travesh);
3295 }
3296 return sdest(travesh);
3297}
3298
3299///////////////////////////////////////////////////////////////////////////////
3300// //
3301// Linear algebra operators. //
3302// //
3303///////////////////////////////////////////////////////////////////////////////
3304
3305// dot() returns the dot product: v1 dot v2.
3306inline REAL tetgenmesh::dot(REAL* v1, REAL* v2)
3307{
3308 return v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2];
3309}
3310
3311// cross() computes the cross product: n = v1 cross v2.
3312inline void tetgenmesh::cross(REAL* v1, REAL* v2, REAL* n)
3313{
3314 n[0] = v1[1] * v2[2] - v2[1] * v1[2];
3315 n[1] = -(v1[0] * v2[2] - v2[0] * v1[2]);
3316 n[2] = v1[0] * v2[1] - v2[0] * v1[1];
3317}
3318
3319// distance() computes the Euclidean distance between two points.
3320inline REAL tetgenmesh::distance(REAL* p1, REAL* p2)
3321{
3322 return sqrt((p2[0] - p1[0]) * (p2[0] - p1[0]) +
3323 (p2[1] - p1[1]) * (p2[1] - p1[1]) +
3324 (p2[2] - p1[2]) * (p2[2] - p1[2]));
3325}
3326
3327inline REAL tetgenmesh::norm2(REAL x, REAL y, REAL z)
3328{
3329 return (x) * (x) + (y) * (y) + (z) * (z);
3330}
3331
3332
3333#endif // #ifndef tetgenH
3334
3335