Chaste  Release::2018.1
tetgen.cpp
Go to the documentation of this file.
1 
7 #ifndef UNUSED_OPT
8 #ifdef NDEBUG
9 #define UNUSED_OPT(var) var=var
10 #else //NDEBUG
11 #define UNUSED_OPT(var)
12 #endif //NDEBUG
13 #endif //UNUSED_OPT
14 
15 // This include is here so that the chaste_libs=0 build correctly links predicates.o.
16 #include "predicates.hpp"
17 
19 // //
20 // TetGen //
21 // //
22 // A Quality Tetrahedral Mesh Generator and 3D Delaunay Triangulator //
23 // //
24 // Version 1.4 //
25 // April 16, 2007 //
26 // //
27 // Copyright (C) 2002--2007 //
28 // Hang Si //
29 // Research Group Numerical Mathematics and Scientific Computing //
30 // Weierstrass Institute for Applied Analysis and Stochastics //
31 // Mohrenstr. 39, 10117 Berlin, Germany //
32 // si@wias-berlin.de //
33 // //
34 // TetGen is freely available through the website: http://tetgen.berlios.de. //
35 // It may be copied, modified, and redistributed for non-commercial use. //
36 // Please consult the file LICENSE for the detailed copyright notices. //
37 // //
39 
41 // //
42 // tetgen.cxx //
43 // //
44 // The TetGen library and program. //
45 // //
47 
48 #include "tetgen.h"
49 namespace tetgen
50 { //Added namespace to avoid clash with triangle
51 
53 // //
54 // terminatetetgen() Terminate TetGen with a given exit code. //
55 // //
57 
58 void terminatetetgen(int x)
59 {
60 #ifdef TETLIBRARY
61  throw x;
62 #else
63  exit(x);
64 #endif // #ifdef TETLIBRARY
65 }
66 
67 //
68 // Begin of class 'tetgenio' implementation
69 //
70 
72 // //
73 // initialize() Initialize all variables of 'tetgenio'. //
74 // //
75 // It is called by the only class constructor 'tetgenio()' implicitly. Thus, //
76 // all variables are guaranteed to be initialized. Each array is initialized //
77 // to be a 'NULL' pointer, and its length is equal zero. Some variables have //
78 // their default value, 'firstnumber' equals zero, 'mesh_dim' equals 3, and //
79 // 'numberofcorners' equals 4. Another possible use of this routine is to //
80 // call it before to re-use an object. //
81 // //
83 
84 void tetgenio::initialize()
85 {
86  firstnumber = 0; // Default item index is numbered from Zero.
87  mesh_dim = 3; // Default mesh dimension is 3.
88  useindex = true;
89 
90  pointlist = (REAL *) NULL;
91  pointattributelist = (REAL *) NULL;
92  pointmtrlist = (REAL *) NULL;
93  pointmarkerlist = (int *) NULL;
94  numberofpoints = 0;
95  numberofpointattributes = 0;
96  numberofpointmtrs = 0;
97 
98  tetrahedronlist = (int *) NULL;
99  tetrahedronattributelist = (REAL *) NULL;
100  tetrahedronvolumelist = (REAL *) NULL;
101  neighborlist = (int *) NULL;
102  numberoftetrahedra = 0;
103  numberofcorners = 4; // Default is 4 nodes per element.
104  numberoftetrahedronattributes = 0;
105 
106  trifacelist = (int *) NULL;
107  adjtetlist = (int *) NULL;
108  trifacemarkerlist = (int *) NULL;
109  numberoftrifaces = 0;
110 
111  facetlist = (facet *) NULL;
112  facetmarkerlist = (int *) NULL;
113  numberoffacets = 0;
114 
115  edgelist = (int *) NULL;
116  edgemarkerlist = (int *) NULL;
117  numberofedges = 0;
118 
119  holelist = (REAL *) NULL;
120  numberofholes = 0;
121 
122  regionlist = (REAL *) NULL;
123  numberofregions = 0;
124 
125  facetconstraintlist = (REAL *) NULL;
126  numberoffacetconstraints = 0;
127  segmentconstraintlist = (REAL *) NULL;
128  numberofsegmentconstraints = 0;
129 
130  pbcgrouplist = (pbcgroup *) NULL;
131  numberofpbcgroups = 0;
132 
133  vpointlist = (REAL *) NULL;
134  vedgelist = (voroedge *) NULL;
135  vfacetlist = (vorofacet *) NULL;
136  vcelllist = (int **) NULL;
137  numberofvpoints = 0;
138  numberofvedges = 0;
139  numberofvfacets = 0;
140  numberofvcells = 0;
141 }
142 
144 // //
145 // deinitialize() Free the memory allocated in 'tetgenio'. //
146 // //
147 // It is called by the class destructor '~tetgenio()' implicitly. Hence, the //
148 // occupied memory by arrays of an object will be automatically released on //
149 // the deletion of the object. However, this routine assumes that the memory //
150 // is allocated by C++ memory allocation operator 'new', thus it is freed by //
151 // the C++ array deletor 'delete []'. If one uses the C/C++ library function //
152 // 'malloc()' to allocate memory for arrays, one has to free them with the //
153 // 'free()' function, and call routine 'initialize()' once to disable this //
154 // routine on deletion of the object. //
155 // //
157 
158 void tetgenio::deinitialize()
159 {
160  facet *f;
161  polygon *p;
162  pbcgroup *pg;
163  int i, j;
164 
165  if (pointlist != (REAL *) NULL) {
166  delete [] pointlist;
167  }
168  if (pointattributelist != (REAL *) NULL) {
169  delete [] pointattributelist;
170  }
171  if (pointmtrlist != (REAL *) NULL) {
172  delete [] pointmtrlist;
173  }
174  if (pointmarkerlist != (int *) NULL) {
175  delete [] pointmarkerlist;
176  }
177 
178  if (tetrahedronlist != (int *) NULL) {
179  delete [] tetrahedronlist;
180  }
181  if (tetrahedronattributelist != (REAL *) NULL) {
182  delete [] tetrahedronattributelist;
183  }
184  if (tetrahedronvolumelist != (REAL *) NULL) {
185  delete [] tetrahedronvolumelist;
186  }
187  if (neighborlist != (int *) NULL) {
188  delete [] neighborlist;
189  }
190 
191  if (trifacelist != (int *) NULL) {
192  delete [] trifacelist;
193  }
194  if (adjtetlist != (int *) NULL) {
195  delete [] adjtetlist;
196  }
197  if (trifacemarkerlist != (int *) NULL) {
198  delete [] trifacemarkerlist;
199  }
200 
201  if (edgelist != (int *) NULL) {
202  delete [] edgelist;
203  }
204  if (edgemarkerlist != (int *) NULL) {
205  delete [] edgemarkerlist;
206  }
207 
208  if (facetlist != (facet *) NULL) {
209  for (i = 0; i < numberoffacets; i++) {
210  f = &facetlist[i];
211  for (j = 0; j < f->numberofpolygons; j++) {
212  p = &f->polygonlist[j];
213  delete [] p->vertexlist;
214  }
215  delete [] f->polygonlist;
216  if (f->holelist != (REAL *) NULL) {
217  delete [] f->holelist;
218  }
219  }
220  delete [] facetlist;
221  }
222  if (facetmarkerlist != (int *) NULL) {
223  delete [] facetmarkerlist;
224  }
225 
226  if (holelist != (REAL *) NULL) {
227  delete [] holelist;
228  }
229  if (regionlist != (REAL *) NULL) {
230  delete [] regionlist;
231  }
232  if (facetconstraintlist != (REAL *) NULL) {
233  delete [] facetconstraintlist;
234  }
235  if (segmentconstraintlist != (REAL *) NULL) {
236  delete [] segmentconstraintlist;
237  }
238  if (pbcgrouplist != (pbcgroup *) NULL) {
239  for (i = 0; i < numberofpbcgroups; i++) {
240  pg = &(pbcgrouplist[i]);
241  if (pg->pointpairlist != (int *) NULL) {
242  delete [] pg->pointpairlist;
243  }
244  }
245  delete [] pbcgrouplist;
246  }
247  if (vpointlist != (REAL *) NULL) {
248  delete [] vpointlist;
249  }
250  if (vedgelist != (voroedge *) NULL) {
251  delete [] vedgelist;
252  }
253  if (vfacetlist != (vorofacet *) NULL) {
254  for (i = 0; i < numberofvfacets; i++) {
255  delete [] vfacetlist[i].elist;
256  }
257  delete [] vfacetlist;
258  }
259  if (vcelllist != (int **) NULL) {
260  for (i = 0; i < numberofvcells; i++) {
261  delete [] vcelllist[i];
262  }
263  delete [] vcelllist;
264  }
265 }
266 
268 // //
269 // load_node_call() Load a list of nodes. //
270 // //
271 // It is a support routine for routines: 'load_nodes()', 'load_poly()', and //
272 // 'load_tetmesh()'. 'infile' is the file handle contains the node list. It //
273 // may point to a .node, or .poly or .smesh file. 'markers' indicates each //
274 // node contains an additional marker (integer) or not. 'infilename' is the //
275 // name of the file being read, it is only appeared in error message. //
276 // //
277 // The 'firstnumber' (0 or 1) is automatically determined by the number of //
278 // the first index of the first point. //
279 // //
281 
282 bool tetgenio::load_node_call(FILE* infile, int markers, char* infilename)
283 {
284  char inputline[INPUTLINESIZE];
285  char *stringptr;
286  REAL x, y, z, attrib;
287  int firstnode, currentmarker;
288  int index, attribindex;
289  int i, j;
290 
291  // Initialize 'pointlist', 'pointattributelist', and 'pointmarkerlist'.
292  pointlist = new REAL[numberofpoints * 3];
293  if (pointlist == (REAL *) NULL) {
294  printf("Error: Out of memory.\n");
295  terminatetetgen(1);
296  }
297  if (numberofpointattributes > 0) {
298  pointattributelist = new REAL[numberofpoints * numberofpointattributes];
299  if (pointattributelist == (REAL *) NULL) {
300  printf("Error: Out of memory.\n");
301  terminatetetgen(1);
302  }
303  }
304  if (markers) {
305  pointmarkerlist = new int[numberofpoints];
306  if (pointmarkerlist == (int *) NULL) {
307  printf("Error: Out of memory.\n");
308  terminatetetgen(1);
309  }
310  }
311 
312  // Read the point section.
313  index = 0;
314  attribindex = 0;
315  for (i = 0; i < numberofpoints; i++) {
316  stringptr = readnumberline(inputline, infile, infilename);
317  if (useindex) {
318  if (i == 0) {
319  firstnode = (int) strtol (stringptr, &stringptr, 0);
320  if ((firstnode == 0) || (firstnode == 1)) {
321  firstnumber = firstnode;
322  }
323  }
324  stringptr = findnextnumber(stringptr);
325  } // if (useindex)
326  if (*stringptr == '\0') {
327  printf("Error: Point %d has no x coordinate.\n", firstnumber + i);
328  break;
329  }
330  x = (REAL) strtod(stringptr, &stringptr);
331  stringptr = findnextnumber(stringptr);
332  if (*stringptr == '\0') {
333  printf("Error: Point %d has no y coordinate.\n", firstnumber + i);
334  break;
335  }
336  y = (REAL) strtod(stringptr, &stringptr);
337  if (mesh_dim == 3) {
338  stringptr = findnextnumber(stringptr);
339  if (*stringptr == '\0') {
340  printf("Error: Point %d has no z coordinate.\n", firstnumber + i);
341  break;
342  }
343  z = (REAL) strtod(stringptr, &stringptr);
344  } else {
345  z = 0.0; // mesh_dim == 2;
346  }
347  pointlist[index++] = x;
348  pointlist[index++] = y;
349  pointlist[index++] = z;
350  // Read the point attributes.
351  for (j = 0; j < numberofpointattributes; j++) {
352  stringptr = findnextnumber(stringptr);
353  if (*stringptr == '\0') {
354  attrib = 0.0;
355  } else {
356  attrib = (REAL) strtod(stringptr, &stringptr);
357  }
358  pointattributelist[attribindex++] = attrib;
359  }
360  if (markers) {
361  // Read a point marker.
362  stringptr = findnextnumber(stringptr);
363  if (*stringptr == '\0') {
364  currentmarker = 0;
365  } else {
366  currentmarker = (int) strtol (stringptr, &stringptr, 0);
367  }
368  pointmarkerlist[i] = currentmarker;
369  }
370  }
371  if (i < numberofpoints) {
372  // Failed to read points due to some error.
373  delete [] pointlist;
374  pointlist = (REAL *) NULL;
375  if (markers) {
376  delete [] pointmarkerlist;
377  pointmarkerlist = (int *) NULL;
378  }
379  if (numberofpointattributes > 0) {
380  delete [] pointattributelist;
381  pointattributelist = (REAL *) NULL;
382  }
383  numberofpoints = 0;
384  return false;
385  }
386  return true;
387 }
388 
390 // //
391 // load_node() Load a list of nodes from a .node file. //
392 // //
393 // 'filename' is the inputfile without suffix. The node list is in 'filename.//
394 // node'. On completion, the node list is returned in 'pointlist'. //
395 // //
397 
398 bool tetgenio::load_node(char* filename)
399 {
400  FILE *infile;
401  char innodefilename[FILENAMESIZE];
402  char inputline[INPUTLINESIZE];
403  char *stringptr;
404  int markers;
405 
406  markers = 0;
407  // Assembling the actual file names we want to open.
408  strcpy(innodefilename, filename);
409  strcat(innodefilename, ".node");
410 
411  // Try to open a .node file.
412  infile = fopen(innodefilename, "r");
413  if (infile == (FILE *) NULL) {
414  printf("File I/O Error: Cannot access file %s.\n", innodefilename);
415  return false;
416  }
417  printf("Opening %s.\n", innodefilename);
418  // Read the first line of the file.
419  stringptr = readnumberline(inputline, infile, innodefilename);
420  // Is this list of points generated from rbox?
421  stringptr = strstr(inputline, "rbox");
422  if (stringptr == NULL) {
423  // Read number of points, number of dimensions, number of point
424  // attributes, and number of boundary markers.
425  stringptr = inputline;
426  numberofpoints = (int) strtol (stringptr, &stringptr, 0);
427  stringptr = findnextnumber(stringptr);
428  if (*stringptr == '\0') {
429  mesh_dim = 3;
430  } else {
431  mesh_dim = (int) strtol (stringptr, &stringptr, 0);
432  }
433  stringptr = findnextnumber(stringptr);
434  if (*stringptr == '\0') {
435  numberofpointattributes = 0;
436  } else {
437  numberofpointattributes = (int) strtol (stringptr, &stringptr, 0);
438  }
439  stringptr = findnextnumber(stringptr);
440  if (*stringptr == '\0') {
441  markers = 0;
442  } else {
443  markers = (int) strtol (stringptr, &stringptr, 0);
444  }
445  } else {
446  // It is a rbox (qhull) input file.
447  stringptr = inputline;
448  // Get the dimension.
449  mesh_dim = (int) strtol (stringptr, &stringptr, 0);
450  // Get the number of points.
451  stringptr = readnumberline(inputline, infile, innodefilename);
452  numberofpoints = (int) strtol (stringptr, &stringptr, 0);
453  // There is no index column.
454  useindex = 0;
455  }
456 
457  // if ((mesh_dim != 3) && (mesh_dim != 2)) {
458  // printf("Input error: TetGen only works for 2D & 3D point sets.\n");
459  // fclose(infile);
460  // return false;
461  // }
462  if (numberofpoints < (mesh_dim + 1)) {
463  printf("Input error: TetGen needs at least %d points.\n", mesh_dim + 1);
464  fclose(infile);
465  return false;
466  }
467 
468  // Load the list of nodes.
469  if (!load_node_call(infile, markers, innodefilename)) {
470  fclose(infile);
471  return false;
472  }
473  fclose(infile);
474  return true;
475 }
476 
478 // //
479 // load_pbc() Load a list of pbc groups into 'pbcgrouplist'. //
480 // //
481 // 'filename' is the filename of the original inputfile without suffix. The //
482 // pbc groups are found in file 'filename.pbc'. //
483 // //
484 // This routine will be called both in load_poly() and load_tetmesh(). //
485 // //
487 
488 bool tetgenio::load_pbc(char* filename)
489 {
490  FILE *infile;
491  char pbcfilename[FILENAMESIZE];
492  char inputline[INPUTLINESIZE];
493  char *stringptr;
494  pbcgroup *pg;
495  int index, p1, p2;
496  int i, j, k;
497 
498  // Pbc groups are saved in file "filename.pbc".
499  strcpy(pbcfilename, filename);
500  strcat(pbcfilename, ".pbc");
501  infile = fopen(pbcfilename, "r");
502  if (infile != (FILE *) NULL) {
503  printf("Opening %s.\n", pbcfilename);
504  } else {
505  // No such file. Return.
506  return false;
507  }
508 
509  // Read the number of pbc groups.
510  stringptr = readnumberline(inputline, infile, pbcfilename);
511  numberofpbcgroups = (int) strtol (stringptr, &stringptr, 0);
512  if (numberofpbcgroups == 0) {
513  // It looks this file contains no point.
514  fclose(infile);
515  return false;
516  }
517  // Initialize 'pbcgrouplist';
518  pbcgrouplist = new pbcgroup[numberofpbcgroups];
519 
520  // Read the list of pbc groups.
521  for (i = 0; i < numberofpbcgroups; i++) {
522  pg = &(pbcgrouplist[i]);
523  // Initialize pbcgroup i;
524  pg->numberofpointpairs = 0;
525  pg->pointpairlist = (int *) NULL;
526  // Read 'fmark1', 'fmark2'.
527  stringptr = readnumberline(inputline, infile, pbcfilename);
528  if (*stringptr == '\0') break;
529  pg->fmark1 = (int) strtol(stringptr, &stringptr, 0);
530  stringptr = findnextnumber(stringptr);
531  if (*stringptr == '\0') break;
532  pg->fmark2 = (int) strtol(stringptr, &stringptr, 0);
533  // Read 'transmat'.
534  do {
535  stringptr = readline(inputline, infile, NULL);
536  } while ((*stringptr != '[') && (*stringptr != '\0'));
537  if (*stringptr == '\0') break;
538  for (j = 0; j < 4; j++) {
539  for (k = 0; k < 4; k++) {
540  // Read the entry of [j, k].
541  stringptr = findnextnumber(stringptr);
542  if (*stringptr == '\0') {
543  // Try to read another line.
544  stringptr = readnumberline(inputline, infile, pbcfilename);
545  if (*stringptr == '\0') break;
546  }
547  pg->transmat[j][k] = (REAL) strtod(stringptr, &stringptr);
548  }
549  if (k < 4) break; // Not complete!
550  }
551  if (j < 4) break; // Not complete!
552  // Read 'numberofpointpairs'.
553  stringptr = readnumberline(inputline, infile, pbcfilename);
554  if (*stringptr == '\0') break;
555  pg->numberofpointpairs = (int) strtol(stringptr, &stringptr, 0);
556  if (pg->numberofpointpairs > 0) {
557  pg->pointpairlist = new int[pg->numberofpointpairs * 2];
558  // Read the point pairs.
559  index = 0;
560  for (j = 0; j < pg->numberofpointpairs; j++) {
561  stringptr = readnumberline(inputline, infile, pbcfilename);
562  p1 = (int) strtol(stringptr, &stringptr, 0);
563  stringptr = findnextnumber(stringptr);
564  p2 = (int) strtol(stringptr, &stringptr, 0);
565  pg->pointpairlist[index++] = p1;
566  pg->pointpairlist[index++] = p2;
567  }
568  }
569  }
570  fclose(infile);
571 
572  if (i < numberofpbcgroups) {
573  // Failed to read to additional points due to some error.
574  delete [] pbcgrouplist;
575  pbcgrouplist = (pbcgroup *) NULL;
576  numberofpbcgroups = 0;
577  return false;
578  }
579  return true;
580 }
581 
583 // //
584 // load_var() Load variant constraints applied on facets, segments, nodes.//
585 // //
586 // 'filename' is the filename of the original inputfile without suffix. The //
587 // constraints are found in file 'filename.var'. //
588 // //
590 
591 bool tetgenio::load_var(char* filename)
592 {
593  FILE *infile;
594  char varfilename[FILENAMESIZE];
595  char inputline[INPUTLINESIZE];
596  char *stringptr;
597  int index;
598  int i;
599 
600  // Variant constraints are saved in file "filename.var".
601  strcpy(varfilename, filename);
602  strcat(varfilename, ".var");
603  infile = fopen(varfilename, "r");
604  if (infile != (FILE *) NULL) {
605  printf("Opening %s.\n", varfilename);
606  } else {
607  // No such file. Return.
608  return false;
609  }
610 
611  // Read the facet constraint section.
612  stringptr = readnumberline(inputline, infile, varfilename);
613  if (*stringptr != '\0') {
614  numberoffacetconstraints = (int) strtol (stringptr, &stringptr, 0);
615  } else {
616  numberoffacetconstraints = 0;
617  }
618  if (numberoffacetconstraints > 0) {
619  // Initialize 'facetconstraintlist'.
620  facetconstraintlist = new REAL[numberoffacetconstraints * 2];
621  index = 0;
622  for (i = 0; i < numberoffacetconstraints; i++) {
623  stringptr = readnumberline(inputline, infile, varfilename);
624  stringptr = findnextnumber(stringptr);
625  if (*stringptr == '\0') {
626  printf("Error: facet constraint %d has no facet marker.\n",
627  firstnumber + i);
628  break;
629  } else {
630  facetconstraintlist[index++] = (REAL) strtod(stringptr, &stringptr);
631  }
632  stringptr = findnextnumber(stringptr);
633  if (*stringptr == '\0') {
634  printf("Error: facet constraint %d has no maximum area bound.\n",
635  firstnumber + i);
636  break;
637  } else {
638  facetconstraintlist[index++] = (REAL) strtod(stringptr, &stringptr);
639  }
640  }
641  if (i < numberoffacetconstraints) {
642  // This must be caused by an error.
643  fclose(infile);
644  return false;
645  }
646  }
647 
648  // Read the segment constraint section.
649  stringptr = readnumberline(inputline, infile, varfilename);
650  if (*stringptr != '\0') {
651  numberofsegmentconstraints = (int) strtol (stringptr, &stringptr, 0);
652  } else {
653  numberofsegmentconstraints = 0;
654  }
655  if (numberofsegmentconstraints > 0) {
656  // Initialize 'segmentconstraintlist'.
657  segmentconstraintlist = new REAL[numberofsegmentconstraints * 3];
658  index = 0;
659  for (i = 0; i < numberofsegmentconstraints; i++) {
660  stringptr = readnumberline(inputline, infile, varfilename);
661  stringptr = findnextnumber(stringptr);
662  if (*stringptr == '\0') {
663  printf("Error: segment constraint %d has no frist endpoint.\n",
664  firstnumber + i);
665  break;
666  } else {
667  segmentconstraintlist[index++] = (REAL) strtod(stringptr, &stringptr);
668  }
669  stringptr = findnextnumber(stringptr);
670  if (*stringptr == '\0') {
671  printf("Error: segment constraint %d has no second endpoint.\n",
672  firstnumber + i);
673  break;
674  } else {
675  segmentconstraintlist[index++] = (REAL) strtod(stringptr, &stringptr);
676  }
677  stringptr = findnextnumber(stringptr);
678  if (*stringptr == '\0') {
679  printf("Error: segment constraint %d has no maximum length bound.\n",
680  firstnumber + i);
681  break;
682  } else {
683  segmentconstraintlist[index++] = (REAL) strtod(stringptr, &stringptr);
684  }
685  }
686  if (i < numberofsegmentconstraints) {
687  // This must be caused by an error.
688  fclose(infile);
689  return false;
690  }
691  }
692 
693  fclose(infile);
694  return true;
695 }
696 
698 // //
699 // load_mtr() Load a size specification map from file. //
700 // //
701 // 'filename' is the filename of the original inputfile without suffix. The //
702 // size map is found in file 'filename.mtr'. //
703 // //
705 
706 bool tetgenio::load_mtr(char* filename)
707 {
708  FILE *infile;
709  char mtrfilename[FILENAMESIZE];
710  char inputline[INPUTLINESIZE];
711  char *stringptr;
712  REAL mtr;
713  int mtrindex;
714  int i, j;
715 
716  strcpy(mtrfilename, filename);
717  strcat(mtrfilename, ".mtr");
718  infile = fopen(mtrfilename, "r");
719  if (infile != (FILE *) NULL) {
720  printf("Opening %s.\n", mtrfilename);
721  } else {
722  // No such file. Return.
723  return false;
724  }
725 
726  // Read number of points, number of columns (1, 3, or 6).
727  stringptr = readnumberline(inputline, infile, mtrfilename);
728  stringptr = findnextnumber(stringptr); // Skip number of points.
729  if (*stringptr != '\0') {
730  numberofpointmtrs = (int) strtol (stringptr, &stringptr, 0);
731  }
732  if (numberofpointmtrs == 0) {
733  // Column number doesn't match. Set a default number (1).
734  numberofpointmtrs = 1;
735  }
736 
737  // Allocate space for pointmtrlist.
738  pointmtrlist = new REAL[numberofpoints * numberofpointmtrs];
739  if (pointmtrlist == (REAL *) NULL) {
740  printf("Error: Out of memory.\n");
741  terminatetetgen(1);
742  }
743  mtrindex = 0;
744  for (i = 0; i < numberofpoints; i++) {
745  // Read metrics.
746  stringptr = readnumberline(inputline, infile, mtrfilename);
747  for (j = 0; j < numberofpointmtrs; j++) {
748  if (*stringptr == '\0') {
749  printf("Error: Metric %d is missing value #%d in %s.\n",
750  i + firstnumber, j + 1, mtrfilename);
751  terminatetetgen(1);
752  }
753  mtr = (REAL) strtod(stringptr, &stringptr);
754  pointmtrlist[mtrindex++] = mtr;
755  stringptr = findnextnumber(stringptr);
756  }
757  }
758 
759  fclose(infile);
760  return true;
761 }
762 
764 // //
765 // load_poly() Load a piecewise linear complex from a .poly or .smesh. //
766 // //
767 // 'filename' is the inputfile without suffix. The PLC is in 'filename.poly' //
768 // or 'filename.smesh', and possibly plus 'filename.node' (when the first //
769 // line of the file starts with a zero). //
770 // //
772 
773 bool tetgenio::load_poly(char* filename)
774 {
775  FILE *infile, *polyfile;
776  char innodefilename[FILENAMESIZE];
777  char inpolyfilename[FILENAMESIZE];
778  char insmeshfilename[FILENAMESIZE];
779  char inputline[INPUTLINESIZE];
780  char *stringptr, *infilename;
781  int smesh, markers, currentmarker;
782  int readnodefile, index;
783  int i, j, k;
784 
785  // Assembling the actual file names we want to open.
786  strcpy(innodefilename, filename);
787  strcpy(inpolyfilename, filename);
788  strcpy(insmeshfilename, filename);
789  strcat(innodefilename, ".node");
790  strcat(inpolyfilename, ".poly");
791  strcat(insmeshfilename, ".smesh");
792 
793  // First assume it is a .poly file.
794  smesh = 0;
795  // Try to open a .poly file.
796  polyfile = fopen(inpolyfilename, "r");
797  if (polyfile == (FILE *) NULL) {
798  // .poly doesn't exist! Try to open a .smesh file.
799  polyfile = fopen(insmeshfilename, "r");
800  if (polyfile == (FILE *) NULL) {
801  printf("File I/O Error: Cannot access file %s and %s.\n",
802  inpolyfilename, insmeshfilename);
803  return false;
804  } else {
805  printf("Opening %s.\n", insmeshfilename);
806  }
807  smesh = 1;
808  } else {
809  printf("Opening %s.\n", inpolyfilename);
810  }
811  // Initialize the default values.
812  mesh_dim = 3; // Three-dimemsional accoordinates.
813  numberofpointattributes = 0; // no point attribute.
814  markers = 0; // no boundary marker.
815  // Read number of points, number of dimensions, number of point
816  // attributes, and number of boundary markers.
817  stringptr = readnumberline(inputline, polyfile, inpolyfilename);
818  numberofpoints = (int) strtol (stringptr, &stringptr, 0);
819  stringptr = findnextnumber(stringptr);
820  if (*stringptr != '\0') {
821  mesh_dim = (int) strtol (stringptr, &stringptr, 0);
822  }
823  stringptr = findnextnumber(stringptr);
824  if (*stringptr != '\0') {
825  numberofpointattributes = (int) strtol (stringptr, &stringptr, 0);
826  }
827  stringptr = findnextnumber(stringptr);
828  if (*stringptr != '\0') {
829  markers = (int) strtol (stringptr, &stringptr, 0);
830  }
831  if (numberofpoints > 0) {
832  readnodefile = 0;
833  if (smesh) {
834  infilename = insmeshfilename;
835  } else {
836  infilename = inpolyfilename;
837  }
838  infile = polyfile;
839  } else {
840  // If the .poly or .smesh file claims there are zero points, that
841  // means the points should be read from a separate .node file.
842  readnodefile = 1;
843  infilename = innodefilename;
844  }
845 
846  if (readnodefile) {
847  // Read the points from the .node file.
848  printf("Opening %s.\n", innodefilename);
849  infile = fopen(innodefilename, "r");
850  if (infile == (FILE *) NULL) {
851  printf("File I/O Error: Cannot access file %s.\n", innodefilename);
852  return false;
853  }
854  // Initialize the default values.
855  mesh_dim = 3; // Three-dimemsional accoordinates.
856  numberofpointattributes = 0; // no point attribute.
857  markers = 0; // no boundary marker.
858  // Read number of points, number of dimensions, number of point
859  // attributes, and number of boundary markers.
860  stringptr = readnumberline(inputline, infile, innodefilename);
861  numberofpoints = (int) strtol (stringptr, &stringptr, 0);
862  stringptr = findnextnumber(stringptr);
863  if (*stringptr != '\0') {
864  mesh_dim = (int) strtol (stringptr, &stringptr, 0);
865  }
866  stringptr = findnextnumber(stringptr);
867  if (*stringptr != '\0') {
868  numberofpointattributes = (int) strtol (stringptr, &stringptr, 0);
869  }
870  stringptr = findnextnumber(stringptr);
871  if (*stringptr != '\0') {
872  markers = (int) strtol (stringptr, &stringptr, 0);
873  }
874  }
875 
876  if ((mesh_dim != 3) && (mesh_dim != 2)) {
877  printf("Input error: TetGen only works for 2D & 3D point sets.\n");
878  fclose(infile);
879  return false;
880  }
881  if (numberofpoints < (mesh_dim + 1)) {
882  printf("Input error: TetGen needs at least %d points.\n", mesh_dim + 1);
883  fclose(infile);
884  return false;
885  }
886 
887  // Load the list of nodes.
888  if (!load_node_call(infile, markers, infilename)) {
889  fclose(infile);
890  return false;
891  }
892 
893  if (readnodefile) {
894  fclose(infile);
895  }
896 
897  facet *f;
898  polygon *p;
899 
900  if (mesh_dim == 3) {
901 
902  // Read number of facets and number of boundary markers.
903  stringptr = readnumberline(inputline, polyfile, inpolyfilename);
904  numberoffacets = (int) strtol (stringptr, &stringptr, 0);
905  if (numberoffacets <= 0) {
906  // No facet list, return.
907  fclose(polyfile);
908  return true;
909  }
910  stringptr = findnextnumber(stringptr);
911  if (*stringptr == '\0') {
912  markers = 0; // no boundary marker.
913  } else {
914  markers = (int) strtol (stringptr, &stringptr, 0);
915  }
916 
917  // Initialize the 'facetlist', 'facetmarkerlist'.
918  facetlist = new facet[numberoffacets];
919  if (markers == 1) {
920  facetmarkerlist = new int[numberoffacets];
921  }
922 
923  // Read data into 'facetlist', 'facetmarkerlist'.
924  if (smesh == 0) {
925  // Facets are in .poly file format.
926  for (i = 1; i <= numberoffacets; i++) {
927  f = &(facetlist[i - 1]);
928  init(f);
929  f->numberofholes = 0;
930  currentmarker = 0;
931  // Read number of polygons, number of holes, and a boundary marker.
932  stringptr = readnumberline(inputline, polyfile, inpolyfilename);
933  f->numberofpolygons = (int) strtol (stringptr, &stringptr, 0);
934  stringptr = findnextnumber(stringptr);
935  if (*stringptr != '\0') {
936  f->numberofholes = (int) strtol (stringptr, &stringptr, 0);
937  if (markers == 1) {
938  stringptr = findnextnumber(stringptr);
939  if (*stringptr != '\0') {
940  currentmarker = (int) strtol(stringptr, &stringptr, 0);
941  }
942  }
943  }
944  // Initialize facetmarker if it needs.
945  if (markers == 1) {
946  facetmarkerlist[i - 1] = currentmarker;
947  }
948  // Each facet should has at least one polygon.
949  if (f->numberofpolygons <= 0) {
950  printf("Error: Wrong number of polygon in %d facet.\n", i);
951  break;
952  }
953  // Initialize the 'f->polygonlist'.
954  f->polygonlist = new polygon[f->numberofpolygons];
955  // Go through all polygons, read in their vertices.
956  for (j = 1; j <= f->numberofpolygons; j++) {
957  p = &(f->polygonlist[j - 1]);
958  init(p);
959  // Read number of vertices of this polygon.
960  stringptr = readnumberline(inputline, polyfile, inpolyfilename);
961  p->numberofvertices = (int) strtol(stringptr, &stringptr, 0);
962  if (p->numberofvertices < 1) {
963  printf("Error: Wrong polygon %d in facet %d\n", j, i);
964  break;
965  }
966  // Initialize 'p->vertexlist'.
967  p->vertexlist = new int[p->numberofvertices];
968  // Read all vertices of this polygon.
969  for (k = 1; k <= p->numberofvertices; k++) {
970  stringptr = findnextnumber(stringptr);
971  if (*stringptr == '\0') {
972  // Try to load another non-empty line and continue to read the
973  // rest of vertices.
974  stringptr = readnumberline(inputline, polyfile, inpolyfilename);
975  if (*stringptr == '\0') {
976  printf("Error: Missing %d endpoints of polygon %d in facet %d",
977  p->numberofvertices - k, j, i);
978  break;
979  }
980  }
981  p->vertexlist[k - 1] = (int) strtol (stringptr, &stringptr, 0);
982  }
983  }
984  if (j <= f->numberofpolygons) {
985  // This must be caused by an error. However, there're j - 1
986  // polygons have been read. Reset the 'f->numberofpolygon'.
987  if (j == 1) {
988  // This is the first polygon.
989  delete [] f->polygonlist;
990  }
991  f->numberofpolygons = j - 1;
992  // No hole will be read even it exists.
993  f->numberofholes = 0;
994  break;
995  }
996  // If this facet has hole pints defined, read them.
997  if (f->numberofholes > 0) {
998  // Initialize 'f->holelist'.
999  f->holelist = new REAL[f->numberofholes * 3];
1000  // Read the holes' coordinates.
1001  index = 0;
1002  for (j = 1; j <= f->numberofholes; j++) {
1003  stringptr = readnumberline(inputline, polyfile, inpolyfilename);
1004  for (k = 1; k <= 3; k++) {
1005  stringptr = findnextnumber(stringptr);
1006  if (*stringptr == '\0') {
1007  printf("Error: Hole %d in facet %d has no coordinates", j, i);
1008  break;
1009  }
1010  f->holelist[index++] = (REAL) strtod (stringptr, &stringptr);
1011  }
1012  if (k <= 3) {
1013  // This must be caused by an error.
1014  break;
1015  }
1016  }
1017  if (j <= f->numberofholes) {
1018  // This must be caused by an error.
1019  break;
1020  }
1021  }
1022  }
1023  if (i <= numberoffacets) {
1024  // This must be caused by an error.
1025  numberoffacets = i - 1;
1026  fclose(polyfile);
1027  return false;
1028  }
1029  } else { // poly == 0
1030  // Read the facets from a .smesh file.
1031  for (i = 1; i <= numberoffacets; i++) {
1032  f = &(facetlist[i - 1]);
1033  init(f);
1034  // Initialize 'f->facetlist'. In a .smesh file, each facetlist only
1035  // contains exactly one polygon, no hole.
1036  f->numberofpolygons = 1;
1037  f->polygonlist = new polygon[f->numberofpolygons];
1038  p = &(f->polygonlist[0]);
1039  init(p);
1040  // Read number of vertices of this polygon.
1041  stringptr = readnumberline(inputline, polyfile, insmeshfilename);
1042  p->numberofvertices = (int) strtol (stringptr, &stringptr, 0);
1043  if (p->numberofvertices < 1) {
1044  printf("Error: Wrong number of vertex in facet %d\n", i);
1045  break;
1046  }
1047  // Initialize 'p->vertexlist'.
1048  p->vertexlist = new int[p->numberofvertices];
1049  for (k = 1; k <= p->numberofvertices; k++) {
1050  stringptr = findnextnumber(stringptr);
1051  if (*stringptr == '\0') {
1052  // Try to load another non-empty line and continue to read the
1053  // rest of vertices.
1054  stringptr = readnumberline(inputline, polyfile, inpolyfilename);
1055  if (*stringptr == '\0') {
1056  printf("Error: Missing %d endpoints in facet %d",
1057  p->numberofvertices - k, i);
1058  break;
1059  }
1060  }
1061  p->vertexlist[k - 1] = (int) strtol (stringptr, &stringptr, 0);
1062  }
1063  if (k <= p->numberofvertices) {
1064  // This must be caused by an error.
1065  break;
1066  }
1067  // Read facet's boundary marker at last.
1068  if (markers == 1) {
1069  stringptr = findnextnumber(stringptr);
1070  if (*stringptr == '\0') {
1071  currentmarker = 0;
1072  } else {
1073  currentmarker = (int) strtol(stringptr, &stringptr, 0);
1074  }
1075  facetmarkerlist[i - 1] = currentmarker;
1076  }
1077  }
1078  if (i <= numberoffacets) {
1079  // This must be caused by an error.
1080  numberoffacets = i - 1;
1081  fclose(polyfile);
1082  return false;
1083  }
1084  }
1085 
1086  // Read the hole section.
1087  stringptr = readnumberline(inputline, polyfile, inpolyfilename);
1088  if (*stringptr != '\0') {
1089  numberofholes = (int) strtol (stringptr, &stringptr, 0);
1090  } else {
1091  numberofholes = 0;
1092  }
1093  if (numberofholes > 0) {
1094  // Initialize 'holelist'.
1095  holelist = new REAL[numberofholes * 3];
1096  for (i = 0; i < 3 * numberofholes; i += 3) {
1097  stringptr = readnumberline(inputline, polyfile, inpolyfilename);
1098  stringptr = findnextnumber(stringptr);
1099  if (*stringptr == '\0') {
1100  printf("Error: Hole %d has no x coord.\n", firstnumber + (i / 3));
1101  break;
1102  } else {
1103  holelist[i] = (REAL) strtod(stringptr, &stringptr);
1104  }
1105  stringptr = findnextnumber(stringptr);
1106  if (*stringptr == '\0') {
1107  printf("Error: Hole %d has no y coord.\n", firstnumber + (i / 3));
1108  break;
1109  } else {
1110  holelist[i + 1] = (REAL) strtod(stringptr, &stringptr);
1111  }
1112  stringptr = findnextnumber(stringptr);
1113  if (*stringptr == '\0') {
1114  printf("Error: Hole %d has no z coord.\n", firstnumber + (i / 3));
1115  break;
1116  } else {
1117  holelist[i + 2] = (REAL) strtod(stringptr, &stringptr);
1118  }
1119  }
1120  if (i < 3 * numberofholes) {
1121  // This must be caused by an error.
1122  fclose(polyfile);
1123  return false;
1124  }
1125  }
1126 
1127  // Read the region section. The 'region' section is optional, if we
1128  // don't reach the end-of-file, try read it in.
1129  stringptr = readnumberline(inputline, polyfile, NULL);
1130  if (stringptr != (char *) NULL && *stringptr != '\0') {
1131  numberofregions = (int) strtol (stringptr, &stringptr, 0);
1132  } else {
1133  numberofregions = 0;
1134  }
1135  if (numberofregions > 0) {
1136  // Initialize 'regionlist'.
1137  regionlist = new REAL[numberofregions * 5];
1138  index = 0;
1139  for (i = 0; i < numberofregions; i++) {
1140  stringptr = readnumberline(inputline, polyfile, inpolyfilename);
1141  stringptr = findnextnumber(stringptr);
1142  if (*stringptr == '\0') {
1143  printf("Error: Region %d has no x coordinate.\n", firstnumber + i);
1144  break;
1145  } else {
1146  regionlist[index++] = (REAL) strtod(stringptr, &stringptr);
1147  }
1148  stringptr = findnextnumber(stringptr);
1149  if (*stringptr == '\0') {
1150  printf("Error: Region %d has no y coordinate.\n", firstnumber + i);
1151  break;
1152  } else {
1153  regionlist[index++] = (REAL) strtod(stringptr, &stringptr);
1154  }
1155  stringptr = findnextnumber(stringptr);
1156  if (*stringptr == '\0') {
1157  printf("Error: Region %d has no z coordinate.\n", firstnumber + i);
1158  break;
1159  } else {
1160  regionlist[index++] = (REAL) strtod(stringptr, &stringptr);
1161  }
1162  stringptr = findnextnumber(stringptr);
1163  if (*stringptr == '\0') {
1164  printf("Error: Region %d has no region attrib.\n", firstnumber + i);
1165  break;
1166  } else {
1167  regionlist[index++] = (REAL) strtod(stringptr, &stringptr);
1168  }
1169  stringptr = findnextnumber(stringptr);
1170  if (*stringptr == '\0') {
1171  regionlist[index] = regionlist[index - 1];
1172  } else {
1173  regionlist[index] = (REAL) strtod(stringptr, &stringptr);
1174  }
1175  index++;
1176  }
1177  if (i < numberofregions) {
1178  // This must be caused by an error.
1179  fclose(polyfile);
1180  return false;
1181  }
1182  }
1183 
1184  } else {
1185 
1186  // Read a PSLG from Triangle's poly file.
1187  assert(mesh_dim == 2);
1188  // A PSLG is a facet of a PLC.
1189  numberoffacets = 1;
1190  // Initialize the 'facetlist'.
1191  facetlist = new facet[numberoffacets];
1192  facetmarkerlist = (int *) NULL; // No facet markers.
1193  f = &(facetlist[0]);
1194  init(f);
1195  // Read number of segments.
1196  stringptr = readnumberline(inputline, polyfile, inpolyfilename);
1197  // Segments are degenerate polygons.
1198  f->numberofpolygons = (int) strtol (stringptr, &stringptr, 0);
1199  if (f->numberofpolygons > 0) {
1200  f->polygonlist = new polygon[f->numberofpolygons];
1201  }
1202  // Go through all segments, read in their vertices.
1203  for (j = 0; j < f->numberofpolygons; j++) {
1204  p = &(f->polygonlist[j]);
1205  init(p);
1206  // Read in a segment.
1207  stringptr = readnumberline(inputline, polyfile, inpolyfilename);
1208  stringptr = findnextnumber(stringptr); // Skip its index.
1209  p->numberofvertices = 2; // A segment always has two vertices.
1210  p->vertexlist = new int[p->numberofvertices];
1211  p->vertexlist[0] = (int) strtol (stringptr, &stringptr, 0);
1212  stringptr = findnextnumber(stringptr);
1213  p->vertexlist[1] = (int) strtol (stringptr, &stringptr, 0);
1214  }
1215  // Read number of holes.
1216  stringptr = readnumberline(inputline, polyfile, inpolyfilename);
1217  f->numberofholes = (int) strtol (stringptr, &stringptr, 0);
1218  if (f->numberofholes > 0) {
1219  // Initialize 'f->holelist'.
1220  f->holelist = new REAL[f->numberofholes * 3];
1221  // Read the holes' coordinates.
1222  for (j = 0; j < f->numberofholes; j++) {
1223  // Read a 2D hole point.
1224  stringptr = readnumberline(inputline, polyfile, inpolyfilename);
1225  stringptr = findnextnumber(stringptr); // Skip its index.
1226  f->holelist[j * 3] = (REAL) strtod (stringptr, &stringptr);
1227  stringptr = findnextnumber(stringptr);
1228  f->holelist[j * 3 + 1] = (REAL) strtod (stringptr, &stringptr);
1229  f->holelist[j * 3 + 2] = 0.0; // The z-coord.
1230  }
1231  }
1232  // The regions are skipped.
1233 
1234  }
1235 
1236  // End of reading poly/smesh file.
1237  fclose(polyfile);
1238 
1239  // Try to load a .var file if it exists.
1240  load_var(filename);
1241  // Try to load a .mtr file if it exists.
1242  load_mtr(filename);
1243  // Try to read a .pbc file if it exists.
1244  load_pbc(filename);
1245 
1246  return true;
1247 }
1248 
1250 // //
1251 // load_off() Load a polyhedron described in a .off file. //
1252 // //
1253 // The .off format is one of file formats of the Geomview, an interactive //
1254 // program for viewing and manipulating geometric objects. More information //
1255 // is available form: http://www.geomview.org. //
1256 // //
1257 // 'filename' is a input filename with extension .off or without extension ( //
1258 // the .off will be added in this case). On completion, the polyhedron is //
1259 // returned in 'pointlist' and 'facetlist'. //
1260 // //
1262 
1263 bool tetgenio::load_off(char* filename)
1264 {
1265  FILE *fp;
1266  tetgenio::facet *f;
1267  tetgenio::polygon *p;
1268  char infilename[FILENAMESIZE];
1269  char buffer[INPUTLINESIZE];
1270  char *bufferp;
1271  double *coord;
1272  int nverts = 0, iverts = 0;
1273  int nfaces = 0, ifaces = 0;
1274  int nedges = 0;
1275  int line_count = 0, i;
1276 
1277  strncpy(infilename, filename, 1024 - 1);
1278  infilename[FILENAMESIZE - 1] = '\0';
1279  if (infilename[0] == '\0') {
1280  printf("Error: No filename.\n");
1281  return false;
1282  }
1283  if (strcmp(&infilename[strlen(infilename) - 4], ".off") != 0) {
1284  strcat(infilename, ".off");
1285  }
1286 
1287  if (!(fp = fopen(infilename, "r"))) {
1288  printf("File I/O Error: Unable to open file %s\n", infilename);
1289  return false;
1290  }
1291  printf("Opening %s.\n", infilename);
1292 
1293  // OFF requires the index starts from '0'.
1294  firstnumber = 0;
1295 
1296  while ((bufferp = readline(buffer, fp, &line_count)) != NULL) {
1297  // Check section
1298  if (nverts == 0) {
1299  // Read header
1300  bufferp = strstr(bufferp, "OFF");
1301  if (bufferp != NULL) {
1302  // Read mesh counts
1303  bufferp = findnextnumber(bufferp); // Skip field "OFF".
1304  if (*bufferp == '\0') {
1305  // Read a non-empty line.
1306  bufferp = readline(buffer, fp, &line_count);
1307  }
1308  if ((sscanf(bufferp, "%d%d%d", &nverts, &nfaces, &nedges) != 3)
1309  || (nverts == 0)) {
1310  printf("Syntax error reading header on line %d in file %s\n",
1311  line_count, infilename);
1312  fclose(fp);
1313  return false;
1314  }
1315  // Allocate memory for 'tetgenio'
1316  if (nverts > 0) {
1317  numberofpoints = nverts;
1318  pointlist = new REAL[nverts * 3];
1319  }
1320  if (nfaces > 0) {
1321  numberoffacets = nfaces;
1322  facetlist = new tetgenio::facet[nfaces];
1323  }
1324  }
1325  } else if (iverts < nverts) {
1326  // Read vertex coordinates
1327  coord = &pointlist[iverts * 3];
1328  for (i = 0; i < 3; i++) {
1329  if (*bufferp == '\0') {
1330  printf("Syntax error reading vertex coords on line %d in file %s\n",
1331  line_count, infilename);
1332  fclose(fp);
1333  return false;
1334  }
1335  coord[i] = (REAL) strtod(bufferp, &bufferp);
1336  bufferp = findnextnumber(bufferp);
1337  }
1338  iverts++;
1339  } else if (ifaces < nfaces) {
1340  // Get next face
1341  f = &facetlist[ifaces];
1342  init(f);
1343  // In .off format, each facet has one polygon, no hole.
1344  f->numberofpolygons = 1;
1345  f->polygonlist = new tetgenio::polygon[1];
1346  p = &f->polygonlist[0];
1347  init(p);
1348  // Read the number of vertices, it should be greater than 0.
1349  p->numberofvertices = (int) strtol(bufferp, &bufferp, 0);
1350  if (p->numberofvertices == 0) {
1351  printf("Syntax error reading polygon on line %d in file %s\n",
1352  line_count, infilename);
1353  fclose(fp);
1354  return false;
1355  }
1356  // Allocate memory for face vertices
1357  p->vertexlist = new int[p->numberofvertices];
1358  for (i = 0; i < p->numberofvertices; i++) {
1359  bufferp = findnextnumber(bufferp);
1360  if (*bufferp == '\0') {
1361  printf("Syntax error reading polygon on line %d in file %s\n",
1362  line_count, infilename);
1363  fclose(fp);
1364  return false;
1365  }
1366  p->vertexlist[i] = (int) strtol(bufferp, &bufferp, 0);
1367  }
1368  ifaces++;
1369  } else {
1370  // Should never get here
1371  printf("Found extra text starting at line %d in file %s\n", line_count,
1372  infilename);
1373  break;
1374  }
1375  }
1376 
1377  // Close file
1378  fclose(fp);
1379 
1380  // Check whether read all points
1381  if (iverts != nverts) {
1382  printf("Expected %d vertices, but read only %d vertices in file %s\n",
1383  nverts, iverts, infilename);
1384  return false;
1385  }
1386 
1387  // Check whether read all faces
1388  if (ifaces != nfaces) {
1389  printf("Expected %d faces, but read only %d faces in file %s\n",
1390  nfaces, ifaces, infilename);
1391  return false;
1392  }
1393 
1394  return true;
1395 }
1396 
1398 // //
1399 // load_ply() Load a polyhedron described in a .ply file. //
1400 // //
1401 // 'filename' is the file name with extension .ply or without extension (the //
1402 // .ply will be added in this case). //
1403 // //
1404 // This is a simplified version of reading .ply files, which only reads the //
1405 // set of vertices and the set of faces. Other informations (such as color, //
1406 // material, texture, etc) in .ply file are ignored. Complete routines for //
1407 // reading and writing ,ply files are available from: http://www.cc.gatech. //
1408 // edu/projects/large_models/ply.html. Except the header section, ply file //
1409 // format has exactly the same format for listing vertices and polygons as //
1410 // off file format. //
1411 // //
1412 // On completion, 'pointlist' and 'facetlist' together return the polyhedron.//
1413 // //
1415 
1416 bool tetgenio::load_ply(char* filename)
1417 {
1418  FILE *fp;
1419  tetgenio::facet *f;
1420  tetgenio::polygon *p;
1421  char infilename[FILENAMESIZE];
1422  char buffer[INPUTLINESIZE];
1423  char *bufferp, *str;
1424  double *coord;
1425  int endheader = 0, format = 0;
1426  int nverts = 0, iverts = 0;
1427  int nfaces = 0, ifaces = 0;
1428  int line_count = 0, i;
1429 
1430  strncpy(infilename, filename, FILENAMESIZE - 1);
1431  infilename[FILENAMESIZE - 1] = '\0';
1432  if (infilename[0] == '\0') {
1433  printf("Error: No filename.\n");
1434  return false;
1435  }
1436  if (strcmp(&infilename[strlen(infilename) - 4], ".ply") != 0) {
1437  strcat(infilename, ".ply");
1438  }
1439 
1440  if (!(fp = fopen(infilename, "r"))) {
1441  printf("Error: Unable to open file %s\n", infilename);
1442  return false;
1443  }
1444  printf("Opening %s.\n", infilename);
1445 
1446  // PLY requires the index starts from '0'.
1447  firstnumber = 0;
1448 
1449  while ((bufferp = readline(buffer, fp, &line_count)) != NULL) {
1450  if (!endheader) {
1451  // Find if it is the keyword "end_header".
1452  str = strstr(bufferp, "end_header");
1453  // strstr() is case sensitive.
1454  if (!str) str = strstr(bufferp, "End_header");
1455  if (!str) str = strstr(bufferp, "End_Header");
1456  if (str) {
1457  // This is the end of the header section.
1458  endheader = 1;
1459  continue;
1460  }
1461  // Parse the number of vertices and the number of faces.
1462  if (nverts == 0 || nfaces == 0) {
1463  // Find if it si the keyword "element".
1464  str = strstr(bufferp, "element");
1465  if (!str) str = strstr(bufferp, "Element");
1466  if (str) {
1467  bufferp = findnextfield(str);
1468  if (*bufferp == '\0') {
1469  printf("Syntax error reading element type on line%d in file %s\n",
1470  line_count, infilename);
1471  fclose(fp);
1472  return false;
1473  }
1474  if (nverts == 0) {
1475  // Find if it is the keyword "vertex".
1476  str = strstr(bufferp, "vertex");
1477  if (!str) str = strstr(bufferp, "Vertex");
1478  if (str) {
1479  bufferp = findnextnumber(str);
1480  if (*bufferp == '\0') {
1481  printf("Syntax error reading vertex number on line");
1482  printf(" %d in file %s\n", line_count, infilename);
1483  fclose(fp);
1484  return false;
1485  }
1486  nverts = (int) strtol(bufferp, &bufferp, 0);
1487  // Allocate memory for 'tetgenio'
1488  if (nverts > 0) {
1489  numberofpoints = nverts;
1490  pointlist = new REAL[nverts * 3];
1491  }
1492  }
1493  }
1494  if (nfaces == 0) {
1495  // Find if it is the keyword "face".
1496  str = strstr(bufferp, "face");
1497  if (!str) str = strstr(bufferp, "Face");
1498  if (str) {
1499  bufferp = findnextnumber(str);
1500  if (*bufferp == '\0') {
1501  printf("Syntax error reading face number on line");
1502  printf(" %d in file %s\n", line_count, infilename);
1503  fclose(fp);
1504  return false;
1505  }
1506  nfaces = (int) strtol(bufferp, &bufferp, 0);
1507  // Allocate memory for 'tetgenio'
1508  if (nfaces > 0) {
1509  numberoffacets = nfaces;
1510  facetlist = new tetgenio::facet[nfaces];
1511  }
1512  }
1513  }
1514  } // It is not the string "element".
1515  }
1516  if (format == 0) {
1517  // Find the keyword "format".
1518  str = strstr(bufferp, "format");
1519  if (!str) str = strstr(bufferp, "Format");
1520  if (str) {
1521  format = 1;
1522  bufferp = findnextfield(str);
1523  // Find if it is the string "ascii".
1524  str = strstr(bufferp, "ascii");
1525  if (!str) str = strstr(bufferp, "ASCII");
1526  if (!str) {
1527  printf("This routine only reads ascii format of ply files.\n");
1528  printf("Hint: You can convert the binary to ascii format by\n");
1529  printf(" using the provided ply tools:\n");
1530  printf(" ply2ascii < %s > ascii_%s\n", infilename, infilename);
1531  fclose(fp);
1532  return false;
1533  }
1534  }
1535  }
1536  } else if (iverts < nverts) {
1537  // Read vertex coordinates
1538  coord = &pointlist[iverts * 3];
1539  for (i = 0; i < 3; i++) {
1540  if (*bufferp == '\0') {
1541  printf("Syntax error reading vertex coords on line %d in file %s\n",
1542  line_count, infilename);
1543  fclose(fp);
1544  return false;
1545  }
1546  coord[i] = (REAL) strtod(bufferp, &bufferp);
1547  bufferp = findnextnumber(bufferp);
1548  }
1549  iverts++;
1550  } else if (ifaces < nfaces) {
1551  // Get next face
1552  f = &facetlist[ifaces];
1553  init(f);
1554  // In .off format, each facet has one polygon, no hole.
1555  f->numberofpolygons = 1;
1556  f->polygonlist = new tetgenio::polygon[1];
1557  p = &f->polygonlist[0];
1558  init(p);
1559  // Read the number of vertices, it should be greater than 0.
1560  p->numberofvertices = (int) strtol(bufferp, &bufferp, 0);
1561  if (p->numberofvertices == 0) {
1562  printf("Syntax error reading polygon on line %d in file %s\n",
1563  line_count, infilename);
1564  fclose(fp);
1565  return false;
1566  }
1567  // Allocate memory for face vertices
1568  p->vertexlist = new int[p->numberofvertices];
1569  for (i = 0; i < p->numberofvertices; i++) {
1570  bufferp = findnextnumber(bufferp);
1571  if (*bufferp == '\0') {
1572  printf("Syntax error reading polygon on line %d in file %s\n",
1573  line_count, infilename);
1574  fclose(fp);
1575  return false;
1576  }
1577  p->vertexlist[i] = (int) strtol(bufferp, &bufferp, 0);
1578  }
1579  ifaces++;
1580  } else {
1581  // Should never get here
1582  printf("Found extra text starting at line %d in file %s\n", line_count,
1583  infilename);
1584  break;
1585  }
1586  }
1587 
1588  // Close file
1589  fclose(fp);
1590 
1591  // Check whether read all points
1592  if (iverts != nverts) {
1593  printf("Expected %d vertices, but read only %d vertices in file %s\n",
1594  nverts, iverts, infilename);
1595  return false;
1596  }
1597 
1598  // Check whether read all faces
1599  if (ifaces != nfaces) {
1600  printf("Expected %d faces, but read only %d faces in file %s\n",
1601  nfaces, ifaces, infilename);
1602  return false;
1603  }
1604 
1605  return true;
1606 }
1607 
1609 // //
1610 // load_stl() Load a surface mesh described in a .stl file. //
1611 // //
1612 // 'filename' is the file name with extension .stl or without extension (the //
1613 // .stl will be added in this case). //
1614 // //
1615 // The .stl or stereolithography format is an ASCII or binary file used in //
1616 // manufacturing. It is a list of the triangular surfaces that describe a //
1617 // computer generated solid model. This is the standard input for most rapid //
1618 // prototyping machines. //
1619 // //
1620 // On completion, 'pointlist' and 'facetlist' together return the polyhedron.//
1621 // Note: After load_stl(), there exist many duplicated points in 'pointlist'.//
1622 // They will be unified during the Delaunay tetrahedralization process. //
1623 // //
1625 
1626 bool tetgenio::load_stl(char* filename)
1627 {
1628  FILE *fp;
1629  tetgenmesh::list *plist;
1630  tetgenio::facet *f;
1631  tetgenio::polygon *p;
1632  char infilename[FILENAMESIZE];
1633  char buffer[INPUTLINESIZE];
1634  char *bufferp, *str;
1635  double *coord;
1636  int solid = 0;
1637  int nverts = 0, iverts = 0;
1638  int nfaces = 0;
1639  int line_count = 0, i;
1640 
1641  strncpy(infilename, filename, FILENAMESIZE - 1);
1642  infilename[FILENAMESIZE - 1] = '\0';
1643  if (infilename[0] == '\0') {
1644  printf("Error: No filename.\n");
1645  return false;
1646  }
1647  if (strcmp(&infilename[strlen(infilename) - 4], ".stl") != 0) {
1648  strcat(infilename, ".stl");
1649  }
1650 
1651  if (!(fp = fopen(infilename, "r"))) {
1652  printf("Error: Unable to open file %s\n", infilename);
1653  return false;
1654  }
1655  printf("Opening %s.\n", infilename);
1656 
1657  // STL file has no number of points available. Use a list to read points.
1658  plist = new tetgenmesh::list(sizeof(double) * 3, NULL, 1024);
1659 
1660  while ((bufferp = readline(buffer, fp, &line_count)) != NULL) {
1661  // The ASCII .stl file must start with the lower case keyword solid and
1662  // end with endsolid.
1663  if (solid == 0) {
1664  // Read header
1665  bufferp = strstr(bufferp, "solid");
1666  if (bufferp != NULL) {
1667  solid = 1;
1668  }
1669  } else {
1670  // We're inside the block of the solid.
1671  str = bufferp;
1672  // Is this the end of the solid.
1673  bufferp = strstr(bufferp, "endsolid");
1674  if (bufferp != NULL) {
1675  solid = 0;
1676  } else {
1677  // Read the XYZ coordinates if it is a vertex.
1678  bufferp = str;
1679  bufferp = strstr(bufferp, "vertex");
1680  if (bufferp != NULL) {
1681  coord = (double *) plist->append(NULL);
1682  for (i = 0; i < 3; i++) {
1683  bufferp = findnextnumber(bufferp);
1684  if (*bufferp == '\0') {
1685  printf("Syntax error reading vertex coords on line %d\n",
1686  line_count);
1687  delete plist;
1688  fclose(fp);
1689  return false;
1690  }
1691  coord[i] = (REAL) strtod(bufferp, &bufferp);
1692  }
1693  }
1694  }
1695  }
1696  }
1697  fclose(fp);
1698 
1699  nverts = plist->len();
1700  // nverts should be an integer times 3 (every 3 vertices denote a face).
1701  if (nverts == 0 || (nverts % 3 != 0)) {
1702  printf("Error: Wrong number of vertices in file %s.\n", infilename);
1703  delete plist;
1704  return false;
1705  }
1706  numberofpoints = nverts;
1707  pointlist = new REAL[nverts * 3];
1708  for (i = 0; i < nverts; i++) {
1709  coord = (double *) (* plist)[i];
1710  iverts = i * 3;
1711  pointlist[iverts] = (REAL) coord[0];
1712  pointlist[iverts + 1] = (REAL) coord[1];
1713  pointlist[iverts + 2] = (REAL) coord[2];
1714  }
1715 
1716  nfaces = (int) (nverts / 3);
1717  numberoffacets = nfaces;
1718  facetlist = new tetgenio::facet[nfaces];
1719 
1720  // Default use '1' as the array starting index.
1721  firstnumber = 1;
1722  iverts = firstnumber;
1723  for (i = 0; i < nfaces; i++) {
1724  f = &facetlist[i];
1725  init(f);
1726  // In .stl format, each facet has one polygon, no hole.
1727  f->numberofpolygons = 1;
1728  f->polygonlist = new tetgenio::polygon[1];
1729  p = &f->polygonlist[0];
1730  init(p);
1731  // Each polygon has three vertices.
1732  p->numberofvertices = 3;
1733  p->vertexlist = new int[p->numberofvertices];
1734  p->vertexlist[0] = iverts;
1735  p->vertexlist[1] = iverts + 1;
1736  p->vertexlist[2] = iverts + 2;
1737  iverts += 3;
1738  }
1739 
1740  delete plist;
1741  return true;
1742 }
1743 
1745 // //
1746 // load_medit() Load a surface mesh described in .mesh file. //
1747 // //
1748 // 'filename' is the file name with extension .mesh or without entension ( //
1749 // the .mesh will be added in this case). .mesh is the file format of Medit, //
1750 // a user-friendly interactive mesh viewing program. //
1751 // //
1752 // This routine ONLY reads the sections containing vertices, triangles, and //
1753 // quadrilaters. Other sections (such as tetrahedra, edges, ...) are ignored.//
1754 // //
1756 
1757 bool tetgenio::load_medit(char* filename)
1758 {
1759  FILE *fp;
1760  tetgenio::facet *tmpflist, *f;
1761  tetgenio::polygon *p;
1762  char infilename[FILENAMESIZE];
1763  char buffer[INPUTLINESIZE];
1764  char *bufferp, *str;
1765  double *coord;
1766  int *tmpfmlist;
1767  int dimension = 0;
1768  int nverts = 0;
1769  int nfaces = 0;
1770  int line_count = 0;
1771  int corners = 0; // 3 (triangle) or 4 (quad).
1772  int i, j;
1773 
1774  strncpy(infilename, filename, FILENAMESIZE - 1);
1775  infilename[FILENAMESIZE - 1] = '\0';
1776  if (infilename[0] == '\0') {
1777  printf("Error: No filename.\n");
1778  return false;
1779  }
1780  if (strcmp(&infilename[strlen(infilename) - 5], ".mesh") != 0) {
1781  strcat(infilename, ".mesh");
1782  }
1783 
1784  if (!(fp = fopen(infilename, "r"))) {
1785  printf("Error: Unable to open file %s\n", infilename);
1786  return false;
1787  }
1788  printf("Opening %s.\n", infilename);
1789 
1790  // Default uses the index starts from '1'.
1791  firstnumber = 1;
1792 
1793  while ((bufferp = readline(buffer, fp, &line_count)) != NULL) {
1794  if (*bufferp == '#') continue; // A comment line is skipped.
1795  if (dimension == 0) {
1796  // Find if it is the keyword "Dimension".
1797  str = strstr(bufferp, "Dimension");
1798  if (!str) str = strstr(bufferp, "dimension");
1799  if (!str) str = strstr(bufferp, "DIMENSION");
1800  if (str) {
1801  // Read the dimensions
1802  bufferp = findnextnumber(str); // Skip field "Dimension".
1803  if (*bufferp == '\0') {
1804  // Read a non-empty line.
1805  bufferp = readline(buffer, fp, &line_count);
1806  }
1807  dimension = (int) strtol(bufferp, &bufferp, 0);
1808  if (dimension != 2 && dimension != 3) {
1809  printf("Unknown dimension in file on line %d in file %s\n",
1810  line_count, infilename);
1811  fclose(fp);
1812  return false;
1813  }
1814  mesh_dim = dimension;
1815  }
1816  }
1817  if (nverts == 0) {
1818  // Find if it is the keyword "Vertices".
1819  str = strstr(bufferp, "Vertices");
1820  if (!str) str = strstr(bufferp, "vertices");
1821  if (!str) str = strstr(bufferp, "VERTICES");
1822  if (str) {
1823  // Read the number of vertices.
1824  bufferp = findnextnumber(str); // Skip field "Vertices".
1825  if (*bufferp == '\0') {
1826  // Read a non-empty line.
1827  bufferp = readline(buffer, fp, &line_count);
1828  }
1829  nverts = (int) strtol(bufferp, &bufferp, 0);
1830  // Allocate memory for 'tetgenio'
1831  if (nverts > 0) {
1832  numberofpoints = nverts;
1833  pointlist = new REAL[nverts * 3];
1834  }
1835  // Read the follwoing node list.
1836  for (i = 0; i < nverts; i++) {
1837  bufferp = readline(buffer, fp, &line_count);
1838  if (bufferp == NULL) {
1839  printf("Unexpected end of file on line %d in file %s\n",
1840  line_count, infilename);
1841  fclose(fp);
1842  return false;
1843  }
1844  // Read vertex coordinates
1845  coord = &pointlist[i * 3];
1846  for (j = 0; j < 3; j++) {
1847  if (*bufferp == '\0') {
1848  printf("Syntax error reading vertex coords on line");
1849  printf(" %d in file %s\n", line_count, infilename);
1850  fclose(fp);
1851  return false;
1852  }
1853  if ((j < 2) || (dimension == 3)) {
1854  coord[j] = (REAL) strtod(bufferp, &bufferp);
1855  } else {
1856  assert((j == 2) && (dimension == 2));
1857  coord[j] = 0.0;
1858  }
1859  bufferp = findnextnumber(bufferp);
1860  }
1861  }
1862  continue;
1863  }
1864  }
1865  if (nfaces == 0) {
1866  // Find if it is the keyword "Triangles" or "Quadrilaterals".
1867  corners = 0;
1868  str = strstr(bufferp, "Triangles");
1869  if (!str) str = strstr(bufferp, "triangles");
1870  if (!str) str = strstr(bufferp, "TRIANGLES");
1871  if (str) {
1872  corners = 3;
1873  } else {
1874  str = strstr(bufferp, "Quadrilaterals");
1875  if (!str) str = strstr(bufferp, "quadrilaterals");
1876  if (!str) str = strstr(bufferp, "QUADRILATERALS");
1877  if (str) {
1878  corners = 4;
1879  }
1880  }
1881  if (corners == 3 || corners == 4) {
1882  // Read the number of triangles (or quadrilaterals).
1883  bufferp = findnextnumber(str); // Skip field "Triangles".
1884  if (*bufferp == '\0') {
1885  // Read a non-empty line.
1886  bufferp = readline(buffer, fp, &line_count);
1887  }
1888  nfaces = strtol(bufferp, &bufferp, 0);
1889  // Allocate memory for 'tetgenio'
1890  if (nfaces > 0) {
1891  if (numberoffacets > 0) {
1892  // facetlist has already been allocated. Enlarge arrays.
1893  tmpflist = new tetgenio::facet[numberoffacets + nfaces];
1894  tmpfmlist = new int[numberoffacets + nfaces];
1895  // Copy the data of old arrays into new arrays.
1896  for (i = 0; i < numberoffacets; i++) {
1897  f = &(tmpflist[i]);
1898  tetgenio::init(f);
1899  *f = facetlist[i];
1900  tmpfmlist[i] = facetmarkerlist[i];
1901  }
1902  // Release old arrays.
1903  delete [] facetlist;
1904  delete [] facetmarkerlist;
1905  // Remember the new arrays.
1906  facetlist = tmpflist;
1907  facetmarkerlist = tmpfmlist;
1908  } else {
1909  // This is the first time to allocate facetlist.
1910  facetlist = new tetgenio::facet[nfaces];
1911  facetmarkerlist = new int[nfaces];
1912  }
1913  }
1914  // Read the following list of faces.
1915  for (i = numberoffacets; i < numberoffacets + nfaces; i++) {
1916  bufferp = readline(buffer, fp, &line_count);
1917  if (bufferp == NULL) {
1918  printf("Unexpected end of file on line %d in file %s\n",
1919  line_count, infilename);
1920  fclose(fp);
1921  return false;
1922  }
1923  f = &facetlist[i];
1924  tetgenio::init(f);
1925  // In .mesh format, each facet has one polygon, no hole.
1926  f->numberofpolygons = 1;
1927  f->polygonlist = new tetgenio::polygon[1];
1928  p = &f->polygonlist[0];
1929  tetgenio::init(p);
1930  p->numberofvertices = corners;
1931  // Allocate memory for face vertices
1932  p->vertexlist = new int[p->numberofvertices];
1933  // Read the vertices of the face.
1934  for (j = 0; j < corners; j++) {
1935  if (*bufferp == '\0') {
1936  printf("Syntax error reading face on line %d in file %s\n",
1937  line_count, infilename);
1938  fclose(fp);
1939  return false;
1940  }
1941  p->vertexlist[j] = (int) strtol(bufferp, &bufferp, 0);
1942  if (firstnumber == 1) {
1943  // Check if a '0' index appears.
1944  if (p->vertexlist[j] == 0) {
1945  // The first index is set to be 0.
1946  firstnumber = 0;
1947  }
1948  }
1949  bufferp = findnextnumber(bufferp);
1950  }
1951  // Read the marker of the face if it exists.
1952  facetmarkerlist[i] = 0;
1953  if (*bufferp != '\0') {
1954  facetmarkerlist[i] = (int) strtol(bufferp, &bufferp, 0);
1955  }
1956  }
1957  // Have read in a list of triangles/quads.
1958  numberoffacets += nfaces;
1959  nfaces = 0;
1960  }
1961  }
1962  // if (nverts > 0 && nfaces > 0) break; // Ignore other data.
1963  }
1964 
1965  // Close file
1966  fclose(fp);
1967 
1968  return true;
1969 }
1970 
1972 // //
1973 // load_plc() Load a piecewise linear complex from file. //
1974 // //
1975 // This is main entrance for loading plcs from different file formats into //
1976 // tetgenio. 'filename' is the input file name without extention. 'object' //
1977 // indicates which file format is used to describ the plc. //
1978 // //
1980 
1981 bool tetgenio::load_plc(char* filename, int object)
1982 {
1983  enum tetgenbehavior::objecttype type;
1984 
1985  type = (enum tetgenbehavior::objecttype) object;
1986  switch (type) {
1987  case tetgenbehavior::NODES:
1988  return load_node(filename);
1989  case tetgenbehavior::POLY:
1990  return load_poly(filename);
1991  case tetgenbehavior::OFF:
1992  return load_off(filename);
1993  case tetgenbehavior::PLY:
1994  return load_ply(filename);
1995  case tetgenbehavior::STL:
1996  return load_stl(filename);
1997  case tetgenbehavior::MEDIT:
1998  return load_medit(filename);
1999  default:
2000  return load_poly(filename);
2001  }
2002 }
2003 
2005 // //
2006 // load_tetmesh() Load a tetrahedral mesh from files. //
2007 // //
2008 // 'filename' is the inputfile without suffix. The nodes of the tetrahedral //
2009 // mesh is in "filename.node", the elements is in "filename.ele", if the //
2010 // "filename.face" and "filename.vol" exists, they will also be read. //
2011 // //
2013 
2014 bool tetgenio::load_tetmesh(char* filename)
2015 {
2016  FILE *infile;
2017  char innodefilename[FILENAMESIZE];
2018  char inelefilename[FILENAMESIZE];
2019  char infacefilename[FILENAMESIZE];
2020  char inedgefilename[FILENAMESIZE];
2021  char involfilename[FILENAMESIZE];
2022  char inputline[INPUTLINESIZE];
2023  char *stringptr, *infilename;
2024  REAL attrib, volume;
2025  int volelements;
2026  int markers, corner;
2027  int index, attribindex;
2028  int i, j;
2029 
2030  markers = 0;
2031 
2032  // Assembling the actual file names we want to open.
2033  strcpy(innodefilename, filename);
2034  strcpy(inelefilename, filename);
2035  strcpy(infacefilename, filename);
2036  strcpy(inedgefilename, filename);
2037  strcpy(involfilename, filename);
2038  strcat(innodefilename, ".node");
2039  strcat(inelefilename, ".ele");
2040  strcat(infacefilename, ".face");
2041  strcat(inedgefilename, ".edge");
2042  strcat(involfilename, ".vol");
2043 
2044  // Read the points from a .node file.
2045  infilename = innodefilename;
2046  printf("Opening %s.\n", infilename);
2047  infile = fopen(infilename, "r");
2048  if (infile == (FILE *) NULL) {
2049  printf("File I/O Error: Cannot access file %s.\n", infilename);
2050  return false;
2051  }
2052  // Read the first line of the file.
2053  stringptr = readnumberline(inputline, infile, infilename);
2054  // Is this list of points generated from rbox?
2055  stringptr = strstr(inputline, "rbox");
2056  if (stringptr == NULL) {
2057  // Read number of points, number of dimensions, number of point
2058  // attributes, and number of boundary markers.
2059  stringptr = inputline;
2060  numberofpoints = (int) strtol (stringptr, &stringptr, 0);
2061  stringptr = findnextnumber(stringptr);
2062  if (*stringptr == '\0') {
2063  mesh_dim = 3;
2064  } else {
2065  mesh_dim = (int) strtol (stringptr, &stringptr, 0);
2066  }
2067  stringptr = findnextnumber(stringptr);
2068  if (*stringptr == '\0') {
2069  numberofpointattributes = 0;
2070  } else {
2071  numberofpointattributes = (int) strtol (stringptr, &stringptr, 0);
2072  }
2073  stringptr = findnextnumber(stringptr);
2074  if (*stringptr == '\0') {
2075  markers = 0; // Default value.
2076  } else {
2077  markers = (int) strtol (stringptr, &stringptr, 0);
2078  }
2079  } else {
2080  // It is a rbox (qhull) input file.
2081  stringptr = inputline;
2082  // Get the dimension.
2083  mesh_dim = (int) strtol (stringptr, &stringptr, 0);
2084  // Get the number of points.
2085  stringptr = readnumberline(inputline, infile, infilename);
2086  numberofpoints = (int) strtol (stringptr, &stringptr, 0);
2087  // There is no index column.
2088  useindex = 0;
2089  }
2090 
2091  // Load the list of nodes.
2092  if (!load_node_call(infile, markers, infilename)) {
2093  fclose(infile);
2094  return false;
2095  }
2096  fclose(infile);
2097 
2098  // Read the elements from an .ele file.
2099  if (mesh_dim == 3) {
2100  infilename = inelefilename;
2101  infile = fopen(infilename, "r");
2102  if (infile != (FILE *) NULL) {
2103  printf("Opening %s.\n", infilename);
2104  // Read number of elements, number of corners (4 or 10), number of
2105  // element attributes.
2106  stringptr = readnumberline(inputline, infile, infilename);
2107  numberoftetrahedra = (int) strtol (stringptr, &stringptr, 0);
2108  stringptr = findnextnumber(stringptr);
2109  if (*stringptr == '\0') {
2110  numberofcorners = 4; // Default read 4 nodes per element.
2111  } else {
2112  numberofcorners = (int) strtol(stringptr, &stringptr, 0);
2113  }
2114  stringptr = findnextnumber(stringptr);
2115  if (*stringptr == '\0') {
2116  numberoftetrahedronattributes = 0; // Default no attribute.
2117  } else {
2118  numberoftetrahedronattributes = (int) strtol(stringptr, &stringptr, 0);
2119  }
2120  if (numberofcorners != 4 && numberofcorners != 10) {
2121  printf("Error: Wrong number of corners %d (should be 4 or 10).\n",
2122  numberofcorners);
2123  fclose(infile);
2124  return false;
2125  }
2126  // Allocate memory for tetrahedra.
2127  if (numberoftetrahedra > 0) {
2128  tetrahedronlist = new int[numberoftetrahedra * numberofcorners];
2129  if (tetrahedronlist == (int *) NULL) {
2130  printf("Error: Out of memory.\n");
2131  terminatetetgen(1);
2132  }
2133  // Allocate memory for output tetrahedron attributes if necessary.
2134  if (numberoftetrahedronattributes > 0) {
2135  tetrahedronattributelist = new REAL[numberoftetrahedra *
2136  numberoftetrahedronattributes];
2137  if (tetrahedronattributelist == (REAL *) NULL) {
2138  printf("Error: Out of memory.\n");
2139  terminatetetgen(1);
2140  }
2141  }
2142  }
2143  // Read the list of tetrahedra.
2144  index = 0;
2145  attribindex = 0;
2146  for (i = 0; i < numberoftetrahedra; i++) {
2147  // Read tetrahedron index and the tetrahedron's corners.
2148  stringptr = readnumberline(inputline, infile, infilename);
2149  for (j = 0; j < numberofcorners; j++) {
2150  stringptr = findnextnumber(stringptr);
2151  if (*stringptr == '\0') {
2152  printf("Error: Tetrahedron %d is missing vertex %d in %s.\n",
2153  i + firstnumber, j + 1, infilename);
2154  terminatetetgen(1);
2155  }
2156  corner = (int) strtol(stringptr, &stringptr, 0);
2157  if (corner < firstnumber || corner >= numberofpoints + firstnumber) {
2158  printf("Error: Tetrahedron %d has an invalid vertex index.\n",
2159  i + firstnumber);
2160  terminatetetgen(1);
2161  }
2162  tetrahedronlist[index++] = corner;
2163  }
2164  // Read the tetrahedron's attributes.
2165  for (j = 0; j < numberoftetrahedronattributes; j++) {
2166  stringptr = findnextnumber(stringptr);
2167  if (*stringptr == '\0') {
2168  attrib = 0.0;
2169  } else {
2170  attrib = (REAL) strtod(stringptr, &stringptr);
2171  }
2172  tetrahedronattributelist[attribindex++] = attrib;
2173  }
2174  }
2175  fclose(infile);
2176  }
2177  } // if (meshdim == 3)
2178 
2179  // Read the hullfaces or subfaces from a .face file if it exists.
2180  if (mesh_dim == 3) {
2181  infilename = infacefilename;
2182  } else {
2183  infilename = inelefilename;
2184  }
2185  infile = fopen(infilename, "r");
2186  if (infile != (FILE *) NULL) {
2187  printf("Opening %s.\n", infilename);
2188  // Read number of faces, boundary markers.
2189  stringptr = readnumberline(inputline, infile, infilename);
2190  numberoftrifaces = (int) strtol (stringptr, &stringptr, 0);
2191  stringptr = findnextnumber(stringptr);
2192  if (mesh_dim == 2) {
2193  // Skip a number.
2194  stringptr = findnextnumber(stringptr);
2195  }
2196  if (*stringptr == '\0') {
2197  markers = 0; // Default there is no marker per face.
2198  } else {
2199  markers = (int) strtol (stringptr, &stringptr, 0);
2200  }
2201  if (numberoftrifaces > 0) {
2202  trifacelist = new int[numberoftrifaces * 3];
2203  if (trifacelist == (int *) NULL) {
2204  printf("Error: Out of memory.\n");
2205  terminatetetgen(1);
2206  }
2207  if (markers) {
2208  trifacemarkerlist = new int[numberoftrifaces * 3];
2209  if (trifacemarkerlist == (int *) NULL) {
2210  printf("Error: Out of memory.\n");
2211  terminatetetgen(1);
2212  }
2213  }
2214  }
2215  // Read the list of faces.
2216  index = 0;
2217  for (i = 0; i < numberoftrifaces; i++) {
2218  // Read face index and the face's three corners.
2219  stringptr = readnumberline(inputline, infile, infilename);
2220  for (j = 0; j < 3; j++) {
2221  stringptr = findnextnumber(stringptr);
2222  if (*stringptr == '\0') {
2223  printf("Error: Face %d is missing vertex %d in %s.\n",
2224  i + firstnumber, j + 1, infilename);
2225  terminatetetgen(1);
2226  }
2227  corner = (int) strtol(stringptr, &stringptr, 0);
2228  if (corner < firstnumber || corner >= numberofpoints + firstnumber) {
2229  printf("Error: Face %d has an invalid vertex index.\n",
2230  i + firstnumber);
2231  terminatetetgen(1);
2232  }
2233  trifacelist[index++] = corner;
2234  }
2235  // Read the boundary marker if it exists.
2236  if (markers) {
2237  stringptr = findnextnumber(stringptr);
2238  if (*stringptr == '\0') {
2239  attrib = 0.0;
2240  } else {
2241  attrib = (REAL) strtod(stringptr, &stringptr);
2242  }
2243  trifacemarkerlist[i] = (int) attrib;
2244  }
2245  }
2246  fclose(infile);
2247  }
2248 
2249  // Read the boundary edges from a .edge file if it exists.
2250  infilename = inedgefilename;
2251  infile = fopen(infilename, "r");
2252  if (infile != (FILE *) NULL) {
2253  printf("Opening %s.\n", infilename);
2254  // Read number of boundary edges.
2255  stringptr = readnumberline(inputline, infile, infilename);
2256  numberofedges = (int) strtol (stringptr, &stringptr, 0);
2257  if (numberofedges > 0) {
2258  edgelist = new int[numberofedges * 2];
2259  if (edgelist == (int *) NULL) {
2260  printf("Error: Out of memory.\n");
2261  terminatetetgen(1);
2262  }
2263  }
2264  // Read the list of faces.
2265  index = 0;
2266  for (i = 0; i < numberofedges; i++) {
2267  // Read face index and the edge's two endpoints.
2268  stringptr = readnumberline(inputline, infile, infilename);
2269  for (j = 0; j < 2; j++) {
2270  stringptr = findnextnumber(stringptr);
2271  if (*stringptr == '\0') {
2272  printf("Error: Edge %d is missing vertex %d in %s.\n",
2273  i + firstnumber, j + 1, infilename);
2274  terminatetetgen(1);
2275  }
2276  corner = (int) strtol(stringptr, &stringptr, 0);
2277  if (corner < firstnumber || corner >= numberofpoints + firstnumber) {
2278  printf("Error: Edge %d has an invalid vertex index.\n",
2279  i + firstnumber);
2280  terminatetetgen(1);
2281  }
2282  edgelist[index++] = corner;
2283  }
2284  }
2285  fclose(infile);
2286  }
2287 
2288  // Read the volume constraints from a .vol file if it exists.
2289  infilename = involfilename;
2290  infile = fopen(infilename, "r");
2291  if (infile != (FILE *) NULL) {
2292  printf("Opening %s.\n", infilename);
2293  // Read number of tetrahedra.
2294  stringptr = readnumberline(inputline, infile, infilename);
2295  volelements = (int) strtol (stringptr, &stringptr, 0);
2296  if (volelements != numberoftetrahedra) {
2297  printf("Warning: %s and %s disagree on number of tetrahedra.\n",
2298  inelefilename, involfilename);
2299  volelements = 0;
2300  }
2301  if (volelements > 0) {
2302  tetrahedronvolumelist = new REAL[volelements];
2303  if (tetrahedronvolumelist == (REAL *) NULL) {
2304  printf("Error: Out of memory.\n");
2305  terminatetetgen(1);
2306  }
2307  }
2308  // Read the list of volume constraints.
2309  for (i = 0; i < volelements; i++) {
2310  stringptr = readnumberline(inputline, infile, infilename);
2311  stringptr = findnextnumber(stringptr);
2312  if (*stringptr == '\0') {
2313  volume = -1.0; // No constraint on this tetrahedron.
2314  } else {
2315  volume = (REAL) strtod(stringptr, &stringptr);
2316  }
2317  tetrahedronvolumelist[i] = volume;
2318  }
2319  fclose(infile);
2320  }
2321 
2322  // Try to load a .mtr file if it exists.
2323  load_mtr(filename);
2324  // Try to read a .pbc file if it exists.
2325  load_pbc(filename);
2326 
2327  return true;
2328 }
2329 
2331 // //
2332 // load_voronoi() Load a Voronoi diagram from files. //
2333 // //
2334 // 'filename' is the inputfile without suffix. The Voronoi diagram is read //
2335 // from files: filename.v.node, filename.v.edge, and filename.v.face. //
2336 // //
2338 
2339 bool tetgenio::load_voronoi(char* filename)
2340 {
2341  FILE *infile;
2342  char innodefilename[FILENAMESIZE];
2343  char inedgefilename[FILENAMESIZE];
2344  char inputline[INPUTLINESIZE];
2345  char *stringptr, *infilename;
2346  voroedge *vedge;
2347  REAL x, y, z;
2348  int firstnode, corner;
2349  int index;
2350  int i, j;
2351 
2352  // Assembling the actual file names we want to open.
2353  strcpy(innodefilename, filename);
2354  strcpy(inedgefilename, filename);
2355  strcat(innodefilename, ".v.node");
2356  strcat(inedgefilename, ".v.edge");
2357 
2358  // Read the points from a .v.node file.
2359  infilename = innodefilename;
2360  printf("Opening %s.\n", infilename);
2361  infile = fopen(infilename, "r");
2362  if (infile == (FILE *) NULL) {
2363  printf("File I/O Error: Cannot access file %s.\n", infilename);
2364  return false;
2365  }
2366  // Read the first line of the file.
2367  stringptr = readnumberline(inputline, infile, infilename);
2368  // Is this list of points generated from rbox?
2369  stringptr = strstr(inputline, "rbox");
2370  if (stringptr == NULL) {
2371  // Read number of points, number of dimensions, number of point
2372  // attributes, and number of boundary markers.
2373  stringptr = inputline;
2374  numberofvpoints = (int) strtol (stringptr, &stringptr, 0);
2375  stringptr = findnextnumber(stringptr);
2376  if (*stringptr == '\0') {
2377  mesh_dim = 3; // Default.
2378  } else {
2379  mesh_dim = (int) strtol (stringptr, &stringptr, 0);
2380  }
2381  useindex = 1; // There is an index column.
2382  } else {
2383  // It is a rbox (qhull) input file.
2384  stringptr = inputline;
2385  // Get the dimension.
2386  mesh_dim = (int) strtol (stringptr, &stringptr, 0);
2387  // Get the number of points.
2388  stringptr = readnumberline(inputline, infile, infilename);
2389  numberofvpoints = (int) strtol (stringptr, &stringptr, 0);
2390  useindex = 0; // No index column.
2391  }
2392  // Initialize 'vpointlist'.
2393  vpointlist = new REAL[numberofvpoints * 3];
2394  if (vpointlist == (REAL *) NULL) {
2395  printf("Error: Out of memory.\n");
2396  terminatetetgen(1);
2397  }
2398  // Read the point section.
2399  index = 0;
2400  for (i = 0; i < numberofvpoints; i++) {
2401  stringptr = readnumberline(inputline, infile, infilename);
2402  if (useindex) {
2403  if (i == 0) {
2404  firstnode = (int) strtol (stringptr, &stringptr, 0);
2405  if ((firstnode == 0) || (firstnode == 1)) {
2406  firstnumber = firstnode;
2407  }
2408  }
2409  stringptr = findnextnumber(stringptr);
2410  } // if (useindex)
2411  if (*stringptr == '\0') {
2412  printf("Error: Point %d has no x coordinate.\n", firstnumber + i);
2413  terminatetetgen(1);
2414  }
2415  x = (REAL) strtod(stringptr, &stringptr);
2416  stringptr = findnextnumber(stringptr);
2417  if (*stringptr == '\0') {
2418  printf("Error: Point %d has no y coordinate.\n", firstnumber + i);
2419  terminatetetgen(1);
2420  }
2421  y = (REAL) strtod(stringptr, &stringptr);
2422  if (mesh_dim == 3) {
2423  stringptr = findnextnumber(stringptr);
2424  if (*stringptr == '\0') {
2425  printf("Error: Point %d has no z coordinate.\n", firstnumber + i);
2426  terminatetetgen(1);
2427  }
2428  z = (REAL) strtod(stringptr, &stringptr);
2429  } else {
2430  z = 0.0; // mesh_dim == 2;
2431  }
2432  vpointlist[index++] = x;
2433  vpointlist[index++] = y;
2434  vpointlist[index++] = z;
2435  }
2436  fclose(infile);
2437 
2438  // Read the Voronoi edges from a .v.edge file if it exists.
2439  infilename = inedgefilename;
2440  infile = fopen(infilename, "r");
2441  if (infile != (FILE *) NULL) {
2442  printf("Opening %s.\n", infilename);
2443  // Read number of boundary edges.
2444  stringptr = readnumberline(inputline, infile, infilename);
2445  numberofvedges = (int) strtol (stringptr, &stringptr, 0);
2446  if (numberofvedges > 0) {
2447  vedgelist = new voroedge[numberofvedges];
2448  }
2449  // Read the list of faces.
2450  index = 0;
2451  for (i = 0; i < numberofvedges; i++) {
2452  // Read edge index and the edge's two endpoints.
2453  stringptr = readnumberline(inputline, infile, infilename);
2454  vedge = &(vedgelist[i]);
2455  for (j = 0; j < 2; j++) {
2456  stringptr = findnextnumber(stringptr);
2457  if (*stringptr == '\0') {
2458  printf("Error: Edge %d is missing vertex %d in %s.\n",
2459  i + firstnumber, j + 1, infilename);
2460  terminatetetgen(1);
2461  }
2462  corner = (int) strtol(stringptr, &stringptr, 0);
2463  j == 0 ? vedge->v1 = corner : vedge->v2 = corner;
2464  }
2465  if (vedge->v2 < 0) {
2466  for (j = 0; j < mesh_dim; j++) {
2467  stringptr = findnextnumber(stringptr);
2468  if (*stringptr == '\0') {
2469  printf("Error: Edge %d is missing normal in %s.\n",
2470  i + firstnumber, infilename);
2471  terminatetetgen(1);
2472  }
2473  vedge->vnormal[j] = (REAL) strtod(stringptr, &stringptr);
2474  }
2475  if (mesh_dim == 2) {
2476  vedge->vnormal[2] = 0.0;
2477  }
2478  } else {
2479  vedge->vnormal[0] = 0.0;
2480  vedge->vnormal[1] = 0.0;
2481  vedge->vnormal[2] = 0.0;
2482  }
2483  }
2484  fclose(infile);
2485  }
2486 
2487  return true;
2488 }
2489 
2491 // //
2492 // save_nodes() Save points to a .node file. //
2493 // //
2494 // 'filename' is a string containing the file name without suffix. //
2495 // //
2497 
2498 void tetgenio::save_nodes(char* filename)
2499 {
2500  FILE *fout;
2501  char outnodefilename[FILENAMESIZE];
2502  char outmtrfilename[FILENAMESIZE];
2503  int i, j;
2504 
2505  sprintf(outnodefilename, "%s.node", filename);
2506  printf("Saving nodes to %s\n", outnodefilename);
2507  fout = fopen(outnodefilename, "w");
2508  fprintf(fout, "%d %d %d %d\n", numberofpoints, mesh_dim,
2509  numberofpointattributes, pointmarkerlist != NULL ? 1 : 0);
2510  for (i = 0; i < numberofpoints; i++) {
2511  if (mesh_dim == 2) {
2512  fprintf(fout, "%d %.16g %.16g", i + firstnumber, pointlist[i * 2],
2513  pointlist[i * 2 + 1]);
2514  } else {
2515  fprintf(fout, "%d %.16g %.16g %.16g", i + firstnumber,
2516  pointlist[i * 3], pointlist[i * 3 + 1], pointlist[i * 3 + 2]);
2517  }
2518  for (j = 0; j < numberofpointattributes; j++) {
2519  fprintf(fout, " %.16g",
2520  pointattributelist[i * numberofpointattributes + j]);
2521  }
2522  if (pointmarkerlist != NULL) {
2523  fprintf(fout, " %d", pointmarkerlist[i]);
2524  }
2525  fprintf(fout, "\n");
2526  }
2527  fclose(fout);
2528 
2529  // If the point metrics exist, output them to a .mtr file.
2530  if ((numberofpointmtrs > 0) && (pointmtrlist != (REAL *) NULL)) {
2531  sprintf(outmtrfilename, "%s.mtr", filename);
2532  printf("Saving metrics to %s\n", outmtrfilename);
2533  fout = fopen(outmtrfilename, "w");
2534  fprintf(fout, "%d %d\n", numberofpoints, numberofpointmtrs);
2535  for (i = 0; i < numberofpoints; i++) {
2536  for (j = 0; j < numberofpointmtrs; j++) {
2537  fprintf(fout, "%.16g ", pointmtrlist[i * numberofpointmtrs + j]);
2538  }
2539  fprintf(fout, "\n");
2540  }
2541  fclose(fout);
2542  }
2543 }
2544 
2546 // //
2547 // save_elements() Save elements to a .ele file. //
2548 // //
2549 // 'filename' is a string containing the file name without suffix. //
2550 // //
2552 
2553 void tetgenio::save_elements(char* filename)
2554 {
2555  FILE *fout;
2556  char outelefilename[FILENAMESIZE];
2557  int i, j;
2558 
2559  sprintf(outelefilename, "%s.ele", filename);
2560  printf("Saving elements to %s\n", outelefilename);
2561  fout = fopen(outelefilename, "w");
2562  fprintf(fout, "%d %d %d\n", numberoftetrahedra, numberofcorners,
2563  numberoftetrahedronattributes);
2564  for (i = 0; i < numberoftetrahedra; i++) {
2565  fprintf(fout, "%d", i + firstnumber);
2566  for (j = 0; j < numberofcorners; j++) {
2567  fprintf(fout, " %5d", tetrahedronlist[i * numberofcorners + j]);
2568  }
2569  for (j = 0; j < numberoftetrahedronattributes; j++) {
2570  fprintf(fout, " %g",
2571  tetrahedronattributelist[i * numberoftetrahedronattributes + j]);
2572  }
2573  fprintf(fout, "\n");
2574  }
2575 
2576  fclose(fout);
2577 }
2578 
2580 // //
2581 // save_faces() Save faces to a .face file. //
2582 // //
2583 // 'filename' is a string containing the file name without suffix. //
2584 // //
2586 
2587 void tetgenio::save_faces(char* filename)
2588 {
2589  FILE *fout;
2590  char outfacefilename[FILENAMESIZE];
2591  int i;
2592 
2593  sprintf(outfacefilename, "%s.face", filename);
2594  printf("Saving faces to %s\n", outfacefilename);
2595  fout = fopen(outfacefilename, "w");
2596  fprintf(fout, "%d %d\n", numberoftrifaces,
2597  trifacemarkerlist != NULL ? 1 : 0);
2598  for (i = 0; i < numberoftrifaces; i++) {
2599  fprintf(fout, "%d %5d %5d %5d", i + firstnumber, trifacelist[i * 3],
2600  trifacelist[i * 3 + 1], trifacelist[i * 3 + 2]);
2601  if (trifacemarkerlist != NULL) {
2602  fprintf(fout, " %d", trifacemarkerlist[i]);
2603  }
2604  fprintf(fout, "\n");
2605  }
2606 
2607  fclose(fout);
2608 }
2609 
2611 // //
2612 // save_edges() Save egdes to a .edge file. //
2613 // //
2614 // 'filename' is a string containing the file name without suffix. //
2615 // //
2617 
2618 void tetgenio::save_edges(char* filename)
2619 {
2620  FILE *fout;
2621  char outedgefilename[FILENAMESIZE];
2622  int i;
2623 
2624  sprintf(outedgefilename, "%s.edge", filename);
2625  printf("Saving edges to %s\n", outedgefilename);
2626  fout = fopen(outedgefilename, "w");
2627  fprintf(fout, "%d %d\n", numberofedges, edgemarkerlist != NULL ? 1 : 0);
2628  for (i = 0; i < numberofedges; i++) {
2629  fprintf(fout, "%d %4d %4d", i + firstnumber, edgelist[i * 2],
2630  edgelist[i * 2 + 1]);
2631  if (edgemarkerlist != NULL) {
2632  fprintf(fout, " %d", edgemarkerlist[i]);
2633  }
2634  fprintf(fout, "\n");
2635  }
2636 
2637  fclose(fout);
2638 }
2639 
2641 // //
2642 // save_neighbors() Save egdes to a .neigh file. //
2643 // //
2644 // 'filename' is a string containing the file name without suffix. //
2645 // //
2647 
2648 void tetgenio::save_neighbors(char* filename)
2649 {
2650  FILE *fout;
2651  char outneighborfilename[FILENAMESIZE];
2652  int i;
2653 
2654  sprintf(outneighborfilename, "%s.neigh", filename);
2655  printf("Saving neighbors to %s\n", outneighborfilename);
2656  fout = fopen(outneighborfilename, "w");
2657  fprintf(fout, "%d %d\n", numberoftetrahedra, mesh_dim + 1);
2658  for (i = 0; i < numberoftetrahedra; i++) {
2659  if (mesh_dim == 2) {
2660  fprintf(fout, "%d %5d %5d %5d", i + firstnumber, neighborlist[i * 3],
2661  neighborlist[i * 3 + 1], neighborlist[i * 3 + 2]);
2662  } else {
2663  fprintf(fout, "%d %5d %5d %5d %5d", i + firstnumber,
2664  neighborlist[i * 4], neighborlist[i * 4 + 1],
2665  neighborlist[i * 4 + 2], neighborlist[i * 4 + 3]);
2666  }
2667  fprintf(fout, "\n");
2668  }
2669 
2670  fclose(fout);
2671 }
2672 
2674 // //
2675 // save_poly() Save segments or facets to a .poly file. //
2676 // //
2677 // 'filename' is a string containing the file name without suffix. It only //
2678 // save the facets, holes and regions. The nodes are saved in a .node file //
2679 // by routine save_nodes(). //
2680 // //
2682 
2683 void tetgenio::save_poly(char* filename)
2684 {
2685  FILE *fout;
2686  facet *f;
2687  polygon *p;
2688  char outpolyfilename[FILENAMESIZE];
2689  int i, j, k;
2690 
2691  sprintf(outpolyfilename, "%s.poly", filename);
2692  printf("Saving poly to %s\n", outpolyfilename);
2693  fout = fopen(outpolyfilename, "w");
2694 
2695  // The zero indicates that the vertices are in a separate .node file.
2696  // Followed by number of dimensions, number of vertex attributes,
2697  // and number of boundary markers (zero or one).
2698  fprintf(fout, "%d %d %d %d\n", 0, mesh_dim, numberofpointattributes,
2699  pointmarkerlist != NULL ? 1 : 0);
2700 
2701  // Save segments or facets.
2702  if (mesh_dim == 2) {
2703  // Number of segments, number of boundary markers (zero or one).
2704  fprintf(fout, "%d %d\n", numberofedges, edgemarkerlist != NULL ? 1 : 0);
2705  for (i = 0; i < numberofedges; i++) {
2706  fprintf(fout, "%d %4d %4d", i + firstnumber, edgelist[i * 2],
2707  edgelist[i * 2 + 1]);
2708  if (edgemarkerlist != NULL) {
2709  fprintf(fout, " %d", edgemarkerlist[i]);
2710  }
2711  fprintf(fout, "\n");
2712  }
2713  } else {
2714  // Number of facets, number of boundary markers (zero or one).
2715  fprintf(fout, "%d %d\n", numberoffacets, facetmarkerlist != NULL ? 1 : 0);
2716  for (i = 0; i < numberoffacets; i++) {
2717  f = &(facetlist[i]);
2718  fprintf(fout, "%d %d %d # %d\n", f->numberofpolygons,f->numberofholes,
2719  facetmarkerlist != NULL ? facetmarkerlist[i] : 0, i + firstnumber);
2720  // Output polygons of this facet.
2721  for (j = 0; j < f->numberofpolygons; j++) {
2722  p = &(f->polygonlist[j]);
2723  fprintf(fout, "%d ", p->numberofvertices);
2724  for (k = 0; k < p->numberofvertices; k++) {
2725  if (((k + 1) % 10) == 0) {
2726  fprintf(fout, "\n ");
2727  }
2728  fprintf(fout, " %d", p->vertexlist[k]);
2729  }
2730  fprintf(fout, "\n");
2731  }
2732  // Output holes of this facet.
2733  for (j = 0; j < f->numberofholes; j++) {
2734  fprintf(fout, "%d %.12g %.12g %.12g\n", j + firstnumber,
2735  f->holelist[j * 3], f->holelist[j * 3 + 1], f->holelist[j * 3 + 2]);
2736  }
2737  }
2738  }
2739 
2740  // Save holes.
2741  fprintf(fout, "%d\n", numberofholes);
2742  for (i = 0; i < numberofholes; i++) {
2743  // Output x, y coordinates.
2744  fprintf(fout, "%d %.12g %.12g", i + firstnumber, holelist[i * mesh_dim],
2745  holelist[i * mesh_dim + 1]);
2746  if (mesh_dim == 3) {
2747  // Output z coordinate.
2748  fprintf(fout, " %.12g", holelist[i * mesh_dim + 2]);
2749  }
2750  fprintf(fout, "\n");
2751  }
2752 
2753  // Save regions.
2754  fprintf(fout, "%d\n", numberofregions);
2755  for (i = 0; i < numberofregions; i++) {
2756  if (mesh_dim == 2) {
2757  // Output the index, x, y coordinates, attribute (region number)
2758  // and maximum area constraint (maybe -1).
2759  fprintf(fout, "%d %.12g %.12g %.12g %.12g\n", i + firstnumber,
2760  regionlist[i * 4], regionlist[i * 4 + 1],
2761  regionlist[i * 4 + 2], regionlist[i * 4 + 3]);
2762  } else {
2763  // Output the index, x, y, z coordinates, attribute (region number)
2764  // and maximum volume constraint (maybe -1).
2765  fprintf(fout, "%d %.12g %.12g %.12g %.12g %.12g\n", i + firstnumber,
2766  regionlist[i * 5], regionlist[i * 5 + 1],
2767  regionlist[i * 5 + 2], regionlist[i * 5 + 3],
2768  regionlist[i * 5 + 4]);
2769  }
2770  }
2771 
2772  fclose(fout);
2773 }
2774 
2776 // //
2777 // readline() Read a nonempty line from a file. //
2778 // //
2779 // A line is considered "nonempty" if it contains something more than white //
2780 // spaces. If a line is considered empty, it will be dropped and the next //
2781 // line will be read, this process ends until reaching the end-of-file or a //
2782 // non-empty line. Return NULL if it is the end-of-file, otherwise, return //
2783 // a pointer to the first non-whitespace character of the line. //
2784 // //
2786 
2787 char* tetgenio::readline(char *string, FILE *infile, int *linenumber)
2788 {
2789  char *result;
2790 
2791  // Search for a non-empty line.
2792  do {
2793  result = fgets(string, INPUTLINESIZE - 1, infile);
2794  if (linenumber) (*linenumber)++;
2795  if (result == (char *) NULL) {
2796  return (char *) NULL;
2797  }
2798  // Skip white spaces.
2799  while ((*result == ' ') || (*result == '\t')) result++;
2800  // If it's end of line, read another line and try again.
2801  } while (*result == '\0');
2802  return result;
2803 }
2804 
2806 // //
2807 // findnextfield() Find the next field of a string. //
2808 // //
2809 // Jumps past the current field by searching for whitespace or a comma, then //
2810 // jumps past the whitespace or the comma to find the next field. //
2811 // //
2813 
2814 char* tetgenio::findnextfield(char *string)
2815 {
2816  char *result;
2817 
2818  result = string;
2819  // Skip the current field. Stop upon reaching whitespace or a comma.
2820  while ((*result != '\0') && (*result != ' ') && (*result != '\t') &&
2821  (*result != ',') && (*result != ';')) {
2822  result++;
2823  }
2824  // Now skip the whitespace or the comma, stop at anything else that looks
2825  // like a character, or the end of a line.
2826  while ((*result == ' ') || (*result == '\t') || (*result == ',') ||
2827  (*result == ';')) {
2828  result++;
2829  }
2830  return result;
2831 }
2832 
2834 // //
2835 // readnumberline() Read a nonempty number line from a file. //
2836 // //
2837 // A line is considered "nonempty" if it contains something that looks like //
2838 // a number. Comments (prefaced by `#') are ignored. //
2839 // //
2841 
2842 char* tetgenio::readnumberline(char *string, FILE *infile, char *infilename)
2843 {
2844  char *result;
2845 
2846  // Search for something that looks like a number.
2847  do {
2848  result = fgets(string, INPUTLINESIZE, infile);
2849  if (result == (char *) NULL) {
2850  if (infilename != (char *) NULL) {
2851  printf(" Error: Unexpected end of file in %s.\n", infilename);
2852  terminatetetgen(1);
2853  }
2854  return result;
2855  }
2856  // Skip anything that doesn't look like a number, a comment,
2857  // or the end of a line.
2858  while ((*result != '\0') && (*result != '#')
2859  && (*result != '.') && (*result != '+') && (*result != '-')
2860  && ((*result < '0') || (*result > '9'))) {
2861  result++;
2862  }
2863  // If it's a comment or end of line, read another line and try again.
2864  } while ((*result == '#') || (*result == '\0'));
2865  return result;
2866 }
2867 
2869 // //
2870 // findnextnumber() Find the next field of a number string. //
2871 // //
2872 // Jumps past the current field by searching for whitespace or a comma, then //
2873 // jumps past the whitespace or the comma to find the next field that looks //
2874 // like a number. //
2875 // //
2877 
2878 char* tetgenio::findnextnumber(char *string)
2879 {
2880  char *result;
2881 
2882  result = string;
2883  // Skip the current field. Stop upon reaching whitespace or a comma.
2884  while ((*result != '\0') && (*result != '#') && (*result != ' ') &&
2885  (*result != '\t') && (*result != ',')) {
2886  result++;
2887  }
2888  // Now skip the whitespace and anything else that doesn't look like a
2889  // number, a comment, or the end of a line.
2890  while ((*result != '\0') && (*result != '#')
2891  && (*result != '.') && (*result != '+') && (*result != '-')
2892  && ((*result < '0') || (*result > '9'))) {
2893  result++;
2894  }
2895  // Check for a comment (prefixed with `#').
2896  if (*result == '#') {
2897  *result = '\0';
2898  }
2899  return result;
2900 }
2901 
2902 //
2903 // End of class 'tetgenio' implementation
2904 //
2905 
2906 static REAL PI = 3.14159265358979323846264338327950288419716939937510582;
2907 
2908 //
2909 // Begin of class 'tetgenbehavior' implementation
2910 //
2911 
2913 // //
2914 // tetgenbehavior() Initialize veriables of 'tetgenbehavior'. //
2915 // //
2917 
2918 tetgenbehavior::tetgenbehavior()
2919 {
2920  // Initialize command line switches.
2921  plc = 0;
2922  quality = 0;
2923  refine = 0;
2924  coarse = 0;
2925  metric = 0;
2926  minratio = 2.0;
2927  goodratio = 0.0;
2928  minangle = 20.0;
2929  goodangle = 0.0;
2930  maxdihedral = 165.0;
2931  mindihedral = 5.0;
2932  varvolume = 0;
2933  fixedvolume = 0;
2934  maxvolume = -1.0;
2935  regionattrib = 0;
2936  insertaddpoints = 0;
2937  diagnose = 0;
2938  offcenter = 0;
2939  conformdel = 0;
2940  alpha1 = sqrt(2.0);
2941  alpha2 = 1.0;
2942  alpha3 = 0.6;
2943  zeroindex = 0;
2944  facesout = 0;
2945  edgesout = 0;
2946  neighout = 0;
2947  voroout = 0;
2948  meditview = 0;
2949  gidview = 0;
2950  geomview = 0;
2951  optlevel = 3;
2952  optpasses = 3;
2953  order = 1;
2954  nojettison = 0;
2955  nobound = 0;
2956  nonodewritten = 0;
2957  noelewritten = 0;
2958  nofacewritten = 0;
2959  noiterationnum = 0;
2960  nobisect = 0;
2961  noflip = 0;
2962  steiner = -1;
2963  fliprepair = 1;
2964  nomerge = 0;
2965  docheck = 0;
2966  quiet = 0;
2967  verbose = 0;
2968  useshelles = 0;
2969  epsilon = 1.0e-8;
2970  epsilon2 = 1.0e-5;
2971  object = NONE;
2972  // Initialize strings
2973  commandline[0] = '\0';
2974  infilename[0] = '\0';
2975  outfilename[0] = '\0';
2976  addinfilename[0] = '\0';
2977  bgmeshfilename[0] = '\0';
2978 }
2979 
2981 // //
2982 // versioninfo() Print the version information of TetGen. //
2983 // //
2985 
2986 void tetgenbehavior::versioninfo()
2987 {
2988  printf("Version 1.4.2 (April 16, 2007).\n");
2989  printf("\n");
2990  printf("Copyright (C) 2002 - 2007\n");
2991  printf("Hang Si\n");
2992  printf("Mohrenstr. 39, 10117 Berlin, Germany\n");
2993  printf("si@wias-berlin.de\n");
2994 }
2995 
2997 // //
2998 // syntax() Print list of command line switches and exit the program. //
2999 // //
3001 
3002 void tetgenbehavior::syntax()
3003 {
3004  printf(" tetgen [-prq_Ra_AiMYS_T_dzo_fenvgGOJBNEFICQVh] input_file\n");
3005  printf(" -p Tetrahedralizes a piecewise linear complex (PLC).\n");
3006  printf(" -r Reconstructs a previously generated mesh.\n");
3007  printf(" -q Quality mesh generation (adding new mesh points to ");
3008  printf("improve mesh quality).\n");
3009  printf(" -R Mesh coarsening (deleting redundant mesh points).\n");
3010  printf(" -a Applies a maximum tetrahedron volume constraint.\n");
3011  printf(" -A Assigns attributes to identify tetrahedra in different ");
3012  printf("regions.\n");
3013  printf(" -i Inserts a list of additional points into mesh.\n");
3014  printf(" -M Does not merge coplanar facets.\n");
3015  printf(" -Y Suppresses boundary facets/segments splitting.\n");
3016  printf(" -S Specifies maximum number of added points.\n");
3017  printf(" -T Sets a tolerance for coplanar test (default 1e-8).\n");
3018  printf(" -d Detects self-intersections of facets of the PLC.\n");
3019  printf(" -z Numbers all output items starting from zero.\n");
3020  printf(" -o2 Generates second-order subparametric elements.\n");
3021  printf(" -f Outputs all faces to .face file.");
3022  printf("file.\n");
3023  printf(" -e Outputs all edges to .edge file.\n");
3024  printf(" -n Outputs tetrahedra neighbors to .neigh file.\n");
3025  printf(" -v Outputs Voronoi diagram to files.\n");
3026  printf(" -g Outputs mesh to .mesh file for viewing by Medit.\n");
3027  printf(" -G Outputs mesh to .msh file for viewing by Gid.\n");
3028  printf(" -O Outputs mesh to .off file for viewing by Geomview.\n");
3029  printf(" -J No jettison of unused vertices from output .node file.\n");
3030  printf(" -B Suppresses output of boundary information.\n");
3031  printf(" -N Suppresses output of .node file.\n");
3032  printf(" -E Suppresses output of .ele file.\n");
3033  printf(" -F Suppresses output of .face file.\n");
3034  printf(" -I Suppresses mesh iteration numbers.\n");
3035  printf(" -C Checks the consistency of the final mesh.\n");
3036  printf(" -Q Quiet: No terminal output except errors.\n");
3037  printf(" -V Verbose: Detailed information, more terminal output.\n");
3038  printf(" -h Help: A brief instruction for using TetGen.\n");
3039 }
3040 
3042 // //
3043 // usage() Print a brief instruction for using TetGen. //
3044 // //
3046 
3047 void tetgenbehavior::usage()
3048 {
3049  printf("TetGen\n");
3050  printf("A Quality Tetrahedral Mesh Generator and 3D Delaunay ");
3051  printf("Triangulator\n");
3052  versioninfo();
3053  printf("\n");
3054  printf("What Can TetGen Do?\n");
3055  printf("\n");
3056  printf(" TetGen generates exact Delaunay tetrahedralizations, exact\n");
3057  printf(" constrained Delaunay tetrahedralizations, and quality ");
3058  printf("tetrahedral\n meshes. The latter are nicely graded and whose ");
3059  printf("tetrahedra have\n radius-edge ratio bounded, thus are suitable ");
3060  printf("for finite element and\n finite volume analysis.\n");
3061  printf("\n");
3062  printf("Command Line Syntax:\n");
3063  printf("\n");
3064  printf(" Below is the command line syntax of TetGen with a list of ");
3065  printf("short\n");
3066  printf(" descriptions. Underscores indicate that numbers may optionally\n");
3067  printf(" follow certain switches. Do not leave any space between a ");
3068  printf("switch\n");
3069  printf(" and its numeric parameter. \'input_file\' contains input data\n");
3070  printf(" depending on the switches you supplied which may be a ");
3071  printf(" piecewise\n");
3072  printf(" linear complex or a list of nodes. File formats and detailed\n");
3073  printf(" description of command line switches are found in user's ");
3074  printf("manual.\n");
3075  printf("\n");
3076  syntax();
3077  printf("\n");
3078  printf("Examples of How to Use TetGen:\n");
3079  printf("\n");
3080  printf(" \'tetgen object\' reads vertices from object.node, and writes ");
3081  printf("their\n Delaunay tetrahedralization to object.1.node and ");
3082  printf("object.1.ele.\n");
3083  printf("\n");
3084  printf(" \'tetgen -p object\' reads a PLC from object.poly or object.");
3085  printf("smesh (and\n possibly object.node) and writes its constrained ");
3086  printf("Delaunay\n tetrahedralization to object.1.node, object.1.ele and ");
3087  printf("object.1.face.\n");
3088  printf("\n");
3089  printf(" \'tetgen -pq1.414a.1 object\' reads a PLC from object.poly or\n");
3090  printf(" object.smesh (and possibly object.node), generates a mesh ");
3091  printf("whose\n tetrahedra have radius-edge ratio smaller than 1.414 and ");
3092  printf("have volume\n of 0.1 or less, and writes the mesh to ");
3093  printf("object.1.node, object.1.ele\n and object.1.face.\n");
3094  printf("\n");
3095  printf("Please send bugs/comments to Hang Si <si@wias-berlin.de>\n");
3096 }
3097 
3099 // //
3100 // parse_commandline() Read the command line, identify switches, and set //
3101 // up options and file names. //
3102 // //
3103 // 'argc' and 'argv' are the same parameters passed to the function main() //
3104 // of a C/C++ program. They together represent the command line user invoked //
3105 // from an environment in which TetGen is running. //
3106 // //
3107 // When TetGen is invoked from an environment. 'argc' is nonzero, switches //
3108 // and input filename should be supplied as zero-terminated strings in //
3109 // argv[0] through argv[argc - 1] and argv[0] shall be the name used to //
3110 // invoke TetGen, i.e. "tetgen". Switches are previously started with a //
3111 // dash '-' to identify them from the input filename. //
3112 // //
3113 // When TetGen is called from within another program. 'argc' is set to zero. //
3114 // switches are given in one zero-terminated string (no previous dash is //
3115 // required.), and 'argv' is a pointer points to this string. No input //
3116 // filename is required (usually the input data has been directly created by //
3117 // user in the 'tetgenio' structure). A default filename 'tetgen-tmpfile' //
3118 // will be created for debugging output purpose. //
3119 // //
3121 
3122 bool tetgenbehavior::parse_commandline(int argc, char **argv)
3123 {
3124  int startindex;
3125  int increment;
3126  int meshnumber;
3127  int scount;
3128  int i, j, k;
3129  char workstring[1024];
3130 
3131  // First determine the input style of the switches.
3132  if (argc == 0) {
3133  startindex = 0; // Switches are given without a dash.
3134  argc = 1; // For running the following for-loop once.
3135  commandline[0] = '\0';
3136  } else {
3137  startindex = 1;
3138  strcpy(commandline, argv[0]);
3139  strcat(commandline, " ");
3140  }
3141 
3142  // Rcount used to count the number of '-R' be used.
3143  scount = 0;
3144 
3145  for (i = startindex; i < argc; i++) {
3146  // Remember the command line switches.
3147  strcat(commandline, argv[i]);
3148  strcat(commandline, " ");
3149  if (startindex == 1) {
3150  // Is this string a filename?
3151  if (argv[i][0] != '-') {
3152  strncpy(infilename, argv[i], 1024 - 1);
3153  infilename[1024 - 1] = '\0';
3154  // Go to the next string directly.
3155  continue;
3156  }
3157  }
3158  // Parse the individual switch from the string.
3159  for (j = startindex; argv[i][j] != '\0'; j++) {
3160  if (argv[i][j] == 'p') {
3161  plc = 1;
3162  } else if (argv[i][j] == 'r') {
3163  refine = 1;
3164  } else if (argv[i][j] == 'R') {
3165  coarse = 1;
3166  } else if (argv[i][j] == 'q') {
3167  quality++;
3168  if (((argv[i][j + 1] >= '0') && (argv[i][j + 1] <= '9')) ||
3169  (argv[i][j + 1] == '.')) {
3170  k = 0;
3171  while (((argv[i][j + 1] >= '0') && (argv[i][j + 1] <= '9')) ||
3172  (argv[i][j + 1] == '.')) {
3173  j++;
3174  workstring[k] = argv[i][j];
3175  k++;
3176  }
3177  workstring[k] = '\0';
3178  if (quality == 1) {
3179  minratio = (REAL) strtod(workstring, (char **) NULL);
3180  } else if (quality == 2) {
3181  mindihedral = (REAL) strtod(workstring, (char **) NULL);
3182  } else if (quality == 3) {
3183  maxdihedral = (REAL) strtod(workstring, (char **) NULL);
3184  } else if (quality == 4) {
3185  alpha2 = (REAL) strtod(workstring, (char **) NULL);
3186  } else if (quality == 5) {
3187  alpha1 = (REAL) strtod(workstring, (char **) NULL);
3188  }
3189  }
3190  } else if (argv[i][j] == 'm') {
3191  metric++;
3192  if (((argv[i][j + 1] >= '0') && (argv[i][j + 1] <= '9')) ||
3193  (argv[i][j + 1] == '.')) {
3194  k = 0;
3195  while (((argv[i][j + 1] >= '0') && (argv[i][j + 1] <= '9')) ||
3196  (argv[i][j + 1] == '.')) {
3197  j++;
3198  workstring[k] = argv[i][j];
3199  k++;
3200  }
3201  workstring[k] = '\0';
3202  if (metric == 1) {
3203  alpha1 = (REAL) strtod(workstring, (char **) NULL);
3204  } else if (metric == 2) {
3205  alpha2 = (REAL) strtod(workstring, (char **) NULL);
3206  }
3207  }
3208  } else if (argv[i][j] == 'a') {
3209  if (((argv[i][j + 1] >= '0') && (argv[i][j + 1] <= '9')) ||
3210  (argv[i][j + 1] == '.')) {
3211  fixedvolume = 1;
3212  k = 0;
3213  while (((argv[i][j + 1] >= '0') && (argv[i][j + 1] <= '9')) ||
3214  (argv[i][j + 1] == '.') || (argv[i][j + 1] == 'e') ||
3215  (argv[i][j + 1] == '-') || (argv[i][j + 1] == '+')) {
3216  j++;
3217  workstring[k] = argv[i][j];
3218  k++;
3219  }
3220  workstring[k] = '\0';
3221  maxvolume = (REAL) strtod(workstring, (char **) NULL);
3222  } else {
3223  varvolume = 1;
3224  }
3225  } else if (argv[i][j] == 'A') {
3226  regionattrib++;
3227  } else if (argv[i][j] == 'i') {
3228  insertaddpoints = 1;
3229  } else if (argv[i][j] == 'd') {
3230  diagnose = 1;
3231  } else if (argv[i][j] == 'z') {
3232  zeroindex = 1;
3233  } else if (argv[i][j] == 'f') {
3234  facesout = 1;
3235  } else if (argv[i][j] == 'e') {
3236  edgesout++;
3237  } else if (argv[i][j] == 'n') {
3238  neighout++;
3239  } else if (argv[i][j] == 'v') {
3240  voroout = 1;
3241  } else if (argv[i][j] == 'g') {
3242  meditview = 1;
3243  } else if (argv[i][j] == 'G') {
3244  gidview = 1;
3245  } else if (argv[i][j] == 'O') {
3246  geomview = 1;
3247  } else if (argv[i][j] == 'M') {
3248  nomerge = 1;
3249  } else if (argv[i][j] == 'Y') {
3250  nobisect++;
3251  } else if (argv[i][j] == 'J') {
3252  nojettison = 1;
3253  } else if (argv[i][j] == 'B') {
3254  nobound = 1;
3255  } else if (argv[i][j] == 'N') {
3256  nonodewritten = 1;
3257  } else if (argv[i][j] == 'E') {
3258  noelewritten = 1;
3259  } else if (argv[i][j] == 'F') {
3260  nofacewritten = 1;
3261  } else if (argv[i][j] == 'I') {
3262  noiterationnum = 1;
3263  } else if (argv[i][j] == 'o') {
3264  if (argv[i][j + 1] == '2') {
3265  j++;
3266  order = 2;
3267  }
3268  } else if (argv[i][j] == 'S') {
3269  if (((argv[i][j + 1] >= '0') && (argv[i][j + 1] <= '9')) ||
3270  (argv[i][j + 1] == '.')) {
3271  k = 0;
3272  while (((argv[i][j + 1] >= '0') && (argv[i][j + 1] <= '9')) ||
3273  (argv[i][j + 1] == '.') || (argv[i][j + 1] == 'e') ||
3274  (argv[i][j + 1] == '-') || (argv[i][j + 1] == '+')) {
3275  j++;
3276  workstring[k] = argv[i][j];
3277  k++;
3278  }
3279  workstring[k] = '\0';
3280  steiner = (int) strtol(workstring, (char **) NULL, 0);
3281  }
3282  } else if (argv[i][j] == 's') {
3283  scount++;
3284  if (((argv[i][j + 1] >= '0') && (argv[i][j + 1] <= '9')) ||
3285  (argv[i][j + 1] == '.')) {
3286  k = 0;
3287  while (((argv[i][j + 1] >= '0') && (argv[i][j + 1] <= '9')) ||
3288  (argv[i][j + 1] == '.') || (argv[i][j + 1] == 'e') ||
3289  (argv[i][j + 1] == '-') || (argv[i][j + 1] == '+')) {
3290  j++;
3291  workstring[k] = argv[i][j];
3292  k++;
3293  }
3294  workstring[k] = '\0';
3295  if (scount == 1) {
3296  optlevel = (int) strtol(workstring, (char **) NULL, 0);
3297  } else if (scount == 2) {
3298  optpasses = (int) strtol(workstring, (char **) NULL, 0);
3299  }
3300  }
3301  } else if (argv[i][j] == 'D') {
3302  conformdel++;
3303  } else if (argv[i][j] == 'T') {
3304  if (((argv[i][j + 1] >= '0') && (argv[i][j + 1] <= '9')) ||
3305  (argv[i][j + 1] == '.')) {
3306  k = 0;
3307  while (((argv[i][j + 1] >= '0') && (argv[i][j + 1] <= '9')) ||
3308  (argv[i][j + 1] == '.') || (argv[i][j + 1] == 'e') ||
3309  (argv[i][j + 1] == '-') || (argv[i][j + 1] == '+')) {
3310  j++;
3311  workstring[k] = argv[i][j];
3312  k++;
3313  }
3314  workstring[k] = '\0';
3315  epsilon = (REAL) strtod(workstring, (char **) NULL);
3316  }
3317  } else if (argv[i][j] == 'C') {
3318  docheck++;
3319  } else if (argv[i][j] == 'X') {
3320  fliprepair = 0;
3321  } else if (argv[i][j] == 'Q') {
3322  quiet = 1;
3323  } else if (argv[i][j] == 'V') {
3324  verbose++;
3325  // } else if (argv[i][j] == 'v') {
3326  // versioninfo();
3327  // terminatetetgen(0);
3328  } else if ((argv[i][j] == 'h') || (argv[i][j] == 'H') ||
3329  (argv[i][j] == '?')) {
3330  usage();
3331  terminatetetgen(0);
3332  } else {
3333  printf("Warning: Unknown switch -%c.\n", argv[i][j]);
3334  }
3335  }
3336  }
3337 
3338  if (startindex == 0) {
3339  // Set a temporary filename for debugging output.
3340  strcpy(infilename, "tetgen-tmpfile");
3341  } else {
3342  if (infilename[0] == '\0') {
3343  // No input file name. Print the syntax and exit.
3344  syntax();
3345  terminatetetgen(0);
3346  }
3347  // Recognize the object from file extension if it is available.
3348  if (!strcmp(&infilename[strlen(infilename) - 5], ".node")) {
3349  infilename[strlen(infilename) - 5] = '\0';
3350  object = NODES;
3351  } else if (!strcmp(&infilename[strlen(infilename) - 5], ".poly")) {
3352  infilename[strlen(infilename) - 5] = '\0';
3353  object = POLY;
3354  plc = 1;
3355  } else if (!strcmp(&infilename[strlen(infilename) - 6], ".smesh")) {
3356  infilename[strlen(infilename) - 6] = '\0';
3357  object = POLY;
3358  plc = 1;
3359  } else if (!strcmp(&infilename[strlen(infilename) - 4], ".off")) {
3360  infilename[strlen(infilename) - 4] = '\0';
3361  object = OFF;
3362  plc = 1;
3363  } else if (!strcmp(&infilename[strlen(infilename) - 4], ".ply")) {
3364  infilename[strlen(infilename) - 4] = '\0';
3365  object = PLY;
3366  plc = 1;
3367  } else if (!strcmp(&infilename[strlen(infilename) - 4], ".stl")) {
3368  infilename[strlen(infilename) - 4] = '\0';
3369  object = STL;
3370  plc = 1;
3371  } else if (!strcmp(&infilename[strlen(infilename) - 5], ".mesh")) {
3372  infilename[strlen(infilename) - 5] = '\0';
3373  object = MEDIT;
3374  plc = 1;
3375  } else if (!strcmp(&infilename[strlen(infilename) - 4], ".ele")) {
3376  infilename[strlen(infilename) - 4] = '\0';
3377  object = MESH;
3378  refine = 1;
3379  }
3380  }
3381  plc = plc || diagnose;
3382  useshelles = plc || refine || coarse || quality;
3383  goodratio = minratio;
3384  goodratio *= goodratio;
3385 
3386  // Detect improper combinations of switches.
3387  if (plc && refine) {
3388  printf("Error: Switch -r cannot use together with -p.\n");
3389  return false;
3390  }
3391  if (refine && (plc || noiterationnum)) {
3392  printf("Error: Switches %s cannot use together with -r.\n",
3393  "-p, -d, and -I");
3394  return false;
3395  }
3396  if (diagnose && (quality || insertaddpoints || (order == 2) || neighout
3397  || docheck)) {
3398  printf("Error: Switches %s cannot use together with -d.\n",
3399  "-q, -i, -o2, -n, and -C");
3400  return false;
3401  }
3402 
3403  // Be careful not to allocate space for element area constraints that
3404  // will never be assigned any value (other than the default -1.0).
3405  if (!refine && !plc) {
3406  varvolume = 0;
3407  }
3408  // Be careful not to add an extra attribute to each element unless the
3409  // input supports it (PLC in, but not refining a preexisting mesh).
3410  if (refine || !plc) {
3411  regionattrib = 0;
3412  }
3413  // If '-a' or '-aa' is in use, enable '-q' option too.
3414  if (fixedvolume || varvolume) {
3415  if (quality == 0) {
3416  quality = 1;
3417  }
3418  }
3419  // Calculate the goodangle for testing bad subfaces.
3420  goodangle = cos(minangle * PI / 180.0);
3421  goodangle *= goodangle;
3422 
3423  increment = 0;
3424  strcpy(workstring, infilename);
3425  j = 1;
3426  while (workstring[j] != '\0') {
3427  if ((workstring[j] == '.') && (workstring[j + 1] != '\0')) {
3428  increment = j + 1;
3429  }
3430  j++;
3431  }
3432  meshnumber = 0;
3433  if (increment > 0) {
3434  j = increment;
3435  do {
3436  if ((workstring[j] >= '0') && (workstring[j] <= '9')) {
3437  meshnumber = meshnumber * 10 + (int) (workstring[j] - '0');
3438  } else {
3439  increment = 0;
3440  }
3441  j++;
3442  } while (workstring[j] != '\0');
3443  }
3444  if (noiterationnum) {
3445  strcpy(outfilename, infilename);
3446  } else if (increment == 0) {
3447  strcpy(outfilename, infilename);
3448  strcat(outfilename, ".1");
3449  } else {
3450  workstring[increment] = '%';
3451  workstring[increment + 1] = 'd';
3452  workstring[increment + 2] = '\0';
3453  sprintf(outfilename, workstring, meshnumber + 1);
3454  }
3455  // Additional input file name has the end ".a".
3456  strcpy(addinfilename, infilename);
3457  strcat(addinfilename, ".a");
3458  // Background filename has the form "*.b.ele", "*.b.node", ...
3459  strcpy(bgmeshfilename, infilename);
3460  strcat(bgmeshfilename, ".b");
3461 
3462  return true;
3463 }
3464 
3465 //
3466 // End of class 'tetgenbehavior' implementation
3467 //
3468 
3469 //
3470 // Begin of class 'tetgenmesh' implementation
3471 //
3472 
3473 //
3474 // Begin of class 'list', 'memorypool' and 'link' implementation
3475 //
3476 
3477 // Following are predefined compare functions for primitive data types.
3478 // These functions take two pointers of the corresponding date type,
3479 // perform the comparation. Return -1, 0 or 1 indicating the default
3480 // linear order of two operators.
3481 
3482 // Compare two 'integers'.
3483 int tetgenmesh::compare_2_ints(const void* x, const void* y) {
3484  if (* (int *) x < * (int *) y) {
3485  return -1;
3486  } else if (* (int *) x > * (int *) y) {
3487  return 1;
3488  } else {
3489  return 0;
3490  }
3491 }
3492 
3493 // Compare two 'longs'. Note: in 64-bit machine the 'long' type is 64-bit
3494 // (8-byte) where the 'int' only 32-bit (4-byte).
3495 int tetgenmesh::compare_2_longs(const void* x, const void* y) {
3496  if (* (long *) x < * (long *) y) {
3497  return -1;
3498  } else if (* (long *) x > * (long *) y) {
3499  return 1;
3500  } else {
3501  return 0;
3502  }
3503 }
3504 
3505 // Compare two 'unsigned longs'.
3506 int tetgenmesh::compare_2_unsignedlongs(const void* x, const void* y) {
3507  if (* (unsigned long *) x < * (unsigned long *) y) {
3508  return -1;
3509  } else if (* (unsigned long *) x > * (unsigned long *) y) {
3510  return 1;
3511  } else {
3512  return 0;
3513  }
3514 }
3515 
3517 // //
3518 // set_compfunc() Determine the size of primitive data types and set the //
3519 // corresponding predefined linear order functions. //
3520 // //
3521 // 'str' is a zero-end string indicating a primitive data type, like 'int', //
3522 // 'long' or 'unsigned long'. Every string ending with a '*' is though as a //
3523 // type of pointer and the type 'unsign long' is used for it. //
3524 // //
3525 // When the type of 'str' is determined, the size of this type (in byte) is //
3526 // returned in 'itbytes', and the pointer of corresponding predefined linear //
3527 // order functions is returned in 'pcomp'. //
3528 // //
3530 
3531 void tetgenmesh::set_compfunc(char* str, int* itbytes, compfunc* pcomp)
3532 {
3533  // First figure out whether it is a pointer or not.
3534  if (str[strlen(str) - 1] == '*') {
3535  *itbytes = sizeof(unsigned long);
3536  *pcomp = &compare_2_unsignedlongs;
3537  return;
3538  }
3539  // Then determine other types.
3540  if (strcmp(str, "int") == 0) {
3541  *itbytes = sizeof(int);
3542  *pcomp = &compare_2_ints;
3543  } else if (strcmp(str, "long") == 0) {
3544  *itbytes = sizeof(long);
3545  *pcomp = &compare_2_longs;
3546  } else if (strcmp(str, "unsigned long") == 0) {
3547  *itbytes = sizeof(unsigned long);
3548  *pcomp = &compare_2_unsignedlongs;
3549  } else {
3550  // It is an unknown type.
3551  printf("Error in set_compfunc(): unknown type %s.\n", str);
3552  terminatetetgen(1);
3553  }
3554 }
3555 
3557 // //
3558 // listinit() Initialize a list for storing a data type. //
3559 // //
3560 // Determine the size of each item, set the maximum size allocated at onece, //
3561 // set the expand size in case the list is full, and set the linear order //
3562 // function if it is provided (default is NULL). //
3563 // //
3565 
3566 void tetgenmesh::list::
3567 listinit(int itbytes, compfunc pcomp, int mitems,int exsize)
3568 {
3569 #ifdef SELF_CHECK
3570  assert(itbytes > 0 && mitems > 0 && exsize > 0);
3571 #endif
3572  itembytes = itbytes;
3573  comp = pcomp;
3574  maxitems = mitems;
3575  expandsize = exsize;
3576  base = (char *) malloc(maxitems * itembytes);
3577  if (base == (char *) NULL) {
3578  printf("Error: Out of memory.\n");
3579  terminatetetgen(1);
3580  }
3581  items = 0;
3582 }
3583 
3585 // //
3586 // append() Add a new item at the end of the list. //
3587 // //
3588 // A new space at the end of this list will be allocated for storing the new //
3589 // item. If the memory is not sufficient, reallocation will be performed. If //
3590 // 'appitem' is not NULL, the contents of this pointer will be copied to the //
3591 // new allocated space. Returns the pointer to the new allocated space. //
3592 // //
3594 
3595 void* tetgenmesh::list::append(void *appitem)
3596 {
3597  // Do we have enough space?
3598  if (items == maxitems) {
3599  char* newbase = (char *) realloc(base, (maxitems + expandsize) *
3600  itembytes);
3601  if (newbase == (char *) NULL) {
3602  printf("Error: Out of memory.\n");
3603  terminatetetgen(1);
3604  }
3605  base = newbase;
3606  maxitems += expandsize;
3607  }
3608  if (appitem != (void *) NULL) {
3609  memcpy(base + items * itembytes, appitem, itembytes);
3610  }
3611  items++;
3612  return (void *) (base + (items - 1) * itembytes);
3613 }
3614 
3616 // //
3617 // insert() Insert an item before 'pos' (range from 0 to items - 1). //
3618 // //
3619 // A new space will be inserted at the position 'pos', that is, items lie //
3620 // after pos (including the item at pos) will be moved one space downwords. //
3621 // If 'insitem' is not NULL, its contents will be copied into the new //
3622 // inserted space. Return a pointer to the new inserted space. //
3623 // //
3625 
3626 void* tetgenmesh::list::insert(int pos, void* insitem)
3627 {
3628  if (pos >= items) {
3629  return append(insitem);
3630  }
3631  // Do we have enough space.
3632  if (items == maxitems) {
3633  char* newbase = (char *) realloc(base, (maxitems + expandsize) *
3634  itembytes);
3635  if (newbase == (char *) NULL) {
3636  printf("Error: Out of memory.\n");
3637  terminatetetgen(1);
3638  }
3639  base = newbase;
3640  maxitems += expandsize;
3641  }
3642  // Do block move.
3643  memmove(base + (pos + 1) * itembytes, // dest
3644  base + pos * itembytes, // src
3645  (items - pos) * itembytes); // size in bytes
3646  // Insert the item.
3647  if (insitem != (void *) NULL) {
3648  memcpy(base + pos * itembytes, insitem, itembytes);
3649  }
3650  items++;
3651  return (void *) (base + pos * itembytes);
3652 }
3653 
3655 // //
3656 // del() Delete an item at 'pos' (range from 0 to items - 1). //
3657 // //
3658 // The space at 'pos' will be overlapped by other item. If 'order' is 1, the //
3659 // remaining items of the list have the same order as usual, i.e., items lie //
3660 // after pos will be moved one space upwords. If 'order' is 0, the last item //
3661 // of the list will be moved up to pos. //
3662 // //
3664 
3665 void tetgenmesh::list::del(int pos, int order)
3666 {
3667  // If 'pos' is the last item of the list, nothing need to do.
3668  if (pos >= 0 && pos < items - 1) {
3669  if (order == 1) {
3670  // Do block move.
3671  memmove(base + pos * itembytes, // dest
3672  base + (pos + 1) * itembytes, // src
3673  (items - pos - 1) * itembytes);
3674  } else {
3675  // Use the last item to overlap the del item.
3676  memcpy(base + pos * itembytes, // item at pos
3677  base + (items - 1) * itembytes, // item at last
3678  itembytes);
3679  }
3680  }
3681  if (items > 0) {
3682  items--;
3683  }
3684 }
3685 
3687 // //
3688 // hasitem() Search in this list to find if 'checkitem' exists. //
3689 // //
3690 // This routine assumes that a linear order function has been set. It loops //
3691 // through the entire list, compares each item to 'checkitem'. If it exists, //
3692 // return its position (between 0 to items - 1), otherwise, return -1. //
3693 // //
3695 
3696 int tetgenmesh::list::hasitem(void* checkitem)
3697 {
3698  int i;
3699 
3700  for (i = 0; i < items; i++) {
3701  if (comp != (compfunc) NULL) {
3702  if ((* comp)((void *)(base + i * itembytes), checkitem) == 0) {
3703  return i;
3704  }
3705  }
3706  }
3707  return -1;
3708 }
3709 
3711 // //
3712 // sort() Sort the items with respect to a linear order function. //
3713 // //
3714 // Uses QuickSort routines (qsort) of the standard C/C++ library (stdlib.h). //
3715 // //
3717 
3718 void tetgenmesh::list::sort()
3719 {
3720  qsort((void *) base, (size_t) items, (size_t) itembytes, comp);
3721 }
3722 
3724 // //
3725 // memorypool() The constructors of memorypool. //
3726 // //
3728 
3729 tetgenmesh::memorypool::memorypool()
3730 {
3731  firstblock = nowblock = (void **) NULL;
3732  nextitem = (void *) NULL;
3733  deaditemstack = (void *) NULL;
3734  pathblock = (void **) NULL;
3735  pathitem = (void *) NULL;
3736  itemwordtype = POINTER;
3737  alignbytes = 0;
3738  itembytes = itemwords = 0;
3739  itemsperblock = 0;
3740  items = maxitems = 0l;
3741  unallocateditems = 0;
3742  pathitemsleft = 0;
3743 }
3744 
3745 tetgenmesh::memorypool::
3746 memorypool(int bytecount, int itemcount, enum wordtype wtype, int alignment)
3747 {
3748  poolinit(bytecount, itemcount, wtype, alignment);
3749 }
3750 
3752 // //
3753 // ~memorypool() Free to the operating system all memory taken by a pool. //
3754 // //
3756 
3757 tetgenmesh::memorypool::~memorypool()
3758 {
3759  while (firstblock != (void **) NULL) {
3760  nowblock = (void **) *(firstblock);
3761  free(firstblock);
3762  firstblock = nowblock;
3763  }
3764 }
3765 
3767 // //
3768 // poolinit() Initialize a pool of memory for allocation of items. //
3769 // //
3770 // A `pool' is created whose records have size at least `bytecount'. Items //
3771 // will be allocated in `itemcount'-item blocks. Each item is assumed to be //
3772 // a collection of words, and either pointers or floating-point values are //
3773 // assumed to be the "primary" word type. (The "primary" word type is used //
3774 // to determine alignment of items.) If `alignment' isn't zero, all items //
3775 // will be `alignment'-byte aligned in memory. `alignment' must be either a //
3776 // multiple or a factor of the primary word size; powers of two are safe. //
3777 // `alignment' is normally used to create a few unused bits at the bottom of //
3778 // each item's pointer, in which information may be stored. //
3779 // //
3781 
3782 void tetgenmesh::memorypool::
3783 poolinit(int bytecount, int itemcount, enum wordtype wtype, int alignment)
3784 {
3785  int wordsize;
3786 
3787  // Initialize values in the pool.
3788  itemwordtype = wtype;
3789  wordsize = (itemwordtype == POINTER) ? sizeof(void *) : sizeof(REAL);
3790  // Find the proper alignment, which must be at least as large as:
3791  // - The parameter `alignment'.
3792  // - The primary word type, to avoid unaligned accesses.
3793  // - sizeof(void *), so the stack of dead items can be maintained
3794  // without unaligned accesses.
3795  if (alignment > wordsize) {
3796  alignbytes = alignment;
3797  } else {
3798  alignbytes = wordsize;
3799  }
3800  if ((int) sizeof(void *) > alignbytes) {
3801  alignbytes = (int) sizeof(void *);
3802  }
3803  itemwords = ((bytecount + alignbytes - 1) / alignbytes)
3804  * (alignbytes / wordsize);
3805  itembytes = itemwords * wordsize;
3806  itemsperblock = itemcount;
3807 
3808  // Allocate a block of items. Space for `itemsperblock' items and one
3809  // pointer (to point to the next block) are allocated, as well as space
3810  // to ensure alignment of the items.
3811  firstblock = (void **) malloc(itemsperblock * itembytes + sizeof(void *)
3812  + alignbytes);
3813  if (firstblock == (void **) NULL) {
3814  printf("Error: Out of memory.\n");
3815  terminatetetgen(1);
3816  }
3817  // Set the next block pointer to NULL.
3818  *(firstblock) = (void *) NULL;
3819  restart();
3820 }
3821 
3823 // //
3824 // restart() Deallocate all items in this pool. //
3825 // //
3826 // The pool is returned to its starting state, except that no memory is //
3827 // freed to the operating system. Rather, the previously allocated blocks //
3828 // are ready to be reused. //
3829 // //
3831 
3832 void tetgenmesh::memorypool::restart()
3833 {
3834  unsigned long alignptr;
3835 
3836  items = 0;
3837  maxitems = 0;
3838 
3839  // Set the currently active block.
3840  nowblock = firstblock;
3841  // Find the first item in the pool. Increment by the size of (void *).
3842  alignptr = (unsigned long) (nowblock + 1);
3843  // Align the item on an `alignbytes'-byte boundary.
3844  nextitem = (void *)
3845  (alignptr + (unsigned long) alignbytes -
3846  (alignptr % (unsigned long) alignbytes));
3847  // There are lots of unallocated items left in this block.
3848  unallocateditems = itemsperblock;
3849  // The stack of deallocated items is empty.
3850  deaditemstack = (void *) NULL;
3851 }
3852 
3854 // //
3855 // alloc() Allocate space for an item. //
3856 // //
3858 
3859 void* tetgenmesh::memorypool::alloc()
3860 {
3861  void *newitem;
3862  void **newblock;
3863  unsigned long alignptr;
3864 
3865  // First check the linked list of dead items. If the list is not
3866  // empty, allocate an item from the list rather than a fresh one.
3867  if (deaditemstack != (void *) NULL) {
3868  newitem = deaditemstack; // Take first item in list.
3869  deaditemstack = * (void **) deaditemstack;
3870  } else {
3871  // Check if there are any free items left in the current block.
3872  if (unallocateditems == 0) {
3873  // Check if another block must be allocated.
3874  if (*nowblock == (void *) NULL) {
3875  // Allocate a new block of items, pointed to by the previous block.
3876  newblock = (void **) malloc(itemsperblock * itembytes + sizeof(void *)
3877  + alignbytes);
3878  if (newblock == (void **) NULL) {
3879  printf("Error: Out of memory.\n");
3880  terminatetetgen(1);
3881  }
3882  *nowblock = (void *) newblock;
3883  // The next block pointer is NULL.
3884  *newblock = (void *) NULL;
3885  }
3886  // Move to the new block.
3887  nowblock = (void **) *nowblock;
3888  // Find the first item in the block.
3889  // Increment by the size of (void *).
3890  alignptr = (unsigned long) (nowblock + 1);
3891  // Align the item on an `alignbytes'-byte boundary.
3892  nextitem = (void *)
3893  (alignptr + (unsigned long) alignbytes -
3894  (alignptr % (unsigned long) alignbytes));
3895  // There are lots of unallocated items left in this block.
3896  unallocateditems = itemsperblock;
3897  }
3898  // Allocate a new item.
3899  newitem = nextitem;
3900  // Advance `nextitem' pointer to next free item in block.
3901  if (itemwordtype == POINTER) {
3902  nextitem = (void *) ((void **) nextitem + itemwords);
3903  } else {
3904  nextitem = (void *) ((REAL *) nextitem + itemwords);
3905  }
3906  unallocateditems--;
3907  maxitems++;
3908  }
3909  items++;
3910  return newitem;
3911 }
3912 
3914 // //
3915 // dealloc() Deallocate space for an item. //
3916 // //
3917 // The deallocated space is stored in a queue for later reuse. //
3918 // //
3920 
3921 void tetgenmesh::memorypool::dealloc(void *dyingitem)
3922 {
3923  // Push freshly killed item onto stack.
3924  *((void **) dyingitem) = deaditemstack;
3925  deaditemstack = dyingitem;
3926  items--;
3927 }
3928 
3930 // //
3931 // traversalinit() Prepare to traverse the entire list of items. //
3932 // //
3933 // This routine is used in conjunction with traverse(). //
3934 // //
3936 
3937 void tetgenmesh::memorypool::traversalinit()
3938 {
3939  unsigned long alignptr;
3940 
3941  // Begin the traversal in the first block.
3942  pathblock = firstblock;
3943  // Find the first item in the block. Increment by the size of (void *).
3944  alignptr = (unsigned long) (pathblock + 1);
3945  // Align with item on an `alignbytes'-byte boundary.
3946  pathitem = (void *)
3947  (alignptr + (unsigned long) alignbytes -
3948  (alignptr % (unsigned long) alignbytes));
3949  // Set the number of items left in the current block.
3950  pathitemsleft = itemsperblock;
3951 }
3952 
3954 // //
3955 // traverse() Find the next item in the list. //
3956 // //
3957 // This routine is used in conjunction with traversalinit(). Be forewarned //
3958 // that this routine successively returns all items in the list, including //
3959 // deallocated ones on the deaditemqueue. It's up to you to figure out which //
3960 // ones are actually dead. It can usually be done more space-efficiently by //
3961 // a routine that knows something about the structure of the item. //
3962 // //
3964 
3965 void* tetgenmesh::memorypool::traverse()
3966 {
3967  void *newitem;
3968  unsigned long alignptr;
3969 
3970  // Stop upon exhausting the list of items.
3971  if (pathitem == nextitem) {
3972  return (void *) NULL;
3973  }
3974  // Check whether any untraversed items remain in the current block.
3975  if (pathitemsleft == 0) {
3976  // Find the next block.
3977  pathblock = (void **) *pathblock;
3978  // Find the first item in the block. Increment by the size of (void *).
3979  alignptr = (unsigned long) (pathblock + 1);
3980  // Align with item on an `alignbytes'-byte boundary.
3981  pathitem = (void *)
3982  (alignptr + (unsigned long) alignbytes -
3983  (alignptr % (unsigned long) alignbytes));
3984  // Set the number of items left in the current block.
3985  pathitemsleft = itemsperblock;
3986  }
3987  newitem = pathitem;
3988  // Find the next item in the block.
3989  if (itemwordtype == POINTER) {
3990  pathitem = (void *) ((void **) pathitem + itemwords);
3991  } else {
3992  pathitem = (void *) ((REAL *) pathitem + itemwords);
3993  }
3994  pathitemsleft--;
3995  return newitem;
3996 }
3997 
3999 // //
4000 // linkinit() Initialize a link for storing items. //
4001 // //
4002 // The input parameters are the size of each item, a pointer of a linear //
4003 // order function and the number of items allocating in one memory bulk. //
4004 // //
4006 
4007 void tetgenmesh::link::linkinit(int bytecount, compfunc pcomp, int itemcount)
4008 {
4009 #ifdef SELF_CHECK
4010  assert(bytecount > 0 && itemcount > 0);
4011 #endif
4012  // Remember the real size of each item.
4013  linkitembytes = bytecount;
4014  // Set the linear order function for this link.
4015  comp = pcomp;
4016 
4017  // Call the constructor of 'memorypool' to initialize its variables.
4018  // like: itembytes, itemwords, items, ... Each node has size
4019  // bytecount + 2 * sizeof(void **), and total 'itemcount + 2' (because
4020  // link has additional two nodes 'head' and 'tail').
4021  poolinit(bytecount + 2 * sizeof(void **), itemcount + 2, POINTER, 0);
4022 
4023  // Initial state of this link.
4024  head = (void **) alloc();
4025  tail = (void **) alloc();
4026  *head = (void *) tail;
4027  *(head + 1) = NULL;
4028  *tail = NULL;
4029  *(tail + 1) = (void *) head;
4030  nextlinkitem = *head;
4031  curpos = 1;
4032  linkitems = 0;
4033 }
4034 
4036 // //
4037 // clear() Deallocate all nodes in this link. //
4038 // //
4039 // The link is returned to its starting state, except that no memory is //
4040 // freed to the operating system. Rather, the previously allocated blocks //
4041 // are ready to be reused. //
4042 // //
4044 
4045 void tetgenmesh::link::clear()
4046 {
4047  // Reset the pool.
4048  restart();
4049 
4050  // Initial state of this link.
4051  head = (void **) alloc();
4052  tail = (void **) alloc();
4053  *head = (void *) tail;
4054  *(head + 1) = NULL;
4055  *tail = NULL;
4056  *(tail + 1) = (void *) head;
4057  nextlinkitem = *head;
4058  curpos = 1;
4059  linkitems = 0;
4060 }
4061 
4063 // //
4064 // move() Causes 'nextlinkitem' to traverse the specified number of nodes,//
4065 // updates 'curpos' to be the node to which 'nextlinkitem' points. //
4066 // //
4067 // 'numberofnodes' is a number indicating how many nodes need be traversed //
4068 // (not counter the current node) need be traversed. It may be positive(move //
4069 // forward) or negative (move backward). Return TRUE if it is successful. //
4070 // //
4072 
4073 bool tetgenmesh::link::move(int numberofnodes)
4074 {
4075  void **nownode;
4076  int i;
4077 
4078  nownode = (void **) nextlinkitem;
4079  if (numberofnodes > 0) {
4080  // Move forward.
4081  i = 0;
4082  while ((i < numberofnodes) && *nownode) {
4083  nownode = (void **) *nownode;
4084  i++;
4085  }
4086  if (*nownode == NULL) return false;
4087  nextlinkitem = (void *) nownode;
4088  curpos += numberofnodes;
4089  } else if (numberofnodes < 0) {
4090  // Move backward.
4091  i = 0;
4092  numberofnodes = -numberofnodes;
4093  while ((i < numberofnodes) && *(nownode + 1)) {
4094  nownode = (void **) *(nownode + 1);
4095  i++;
4096  }
4097  if (*(nownode + 1) == NULL) return false;
4098  nextlinkitem = (void *) nownode;
4099  curpos -= numberofnodes;
4100  }
4101  return true;
4102 }
4103 
4105 // //
4106 // locate() Locates the node at the specified position. //
4107 // //
4108 // The number 'pos' (between 1 and 'linkitems') indicates the location. This //
4109 // routine first decides the shortest path traversing from 'curpos' to 'pos',//
4110 // i.e., from head, tail or 'curpos'. Routine 'move()' is called to really //
4111 // traverse the link. If success, 'nextlinkitem' points to the node, 'curpos'//
4112 // and 'pos' are equal. Otherwise, return FALSE. //
4113 // //
4115 
4116 bool tetgenmesh::link::locate(int pos)
4117 {
4118  int headdist, taildist, curdist;
4119  int abscurdist, mindist;
4120 
4121  if (pos < 1 || pos > linkitems) return false;
4122 
4123  headdist = pos - 1;
4124  taildist = linkitems - pos;
4125  curdist = pos - curpos;
4126  abscurdist = curdist >= 0 ? curdist : -curdist;
4127 
4128  if (headdist > taildist) {
4129  if (taildist > abscurdist) {
4130  mindist = curdist;
4131  } else {
4132  // taildist <= abs(curdist)
4133  mindist = -taildist;
4134  goend();
4135  }
4136  } else {
4137  // headdist <= taildist
4138  if (headdist > abscurdist) {
4139  mindist = curdist;
4140  } else {
4141  // headdist <= abs(curdist)
4142  mindist = headdist;
4143  rewind();
4144  }
4145  }
4146 
4147  return move(mindist);
4148 }
4149 
4151 // //
4152 // add() Add a node at the end of this link. //
4153 // //
4154 // A new node is appended to the end of the link. If 'newitem' is not NULL, //
4155 // its conents will be copied to the data slot of the new node. Returns the //
4156 // pointer to the newest added node. //
4157 // //
4159 
4160 void* tetgenmesh::link::add(void* newitem)
4161 {
4162  void **newnode = tail;
4163  if (newitem != (void *) NULL) {
4164  memcpy((void *)(newnode + 2), newitem, linkitembytes);
4165  }
4166  tail = (void **) alloc();
4167  *tail = NULL;
4168  *newnode = (void*) tail;
4169  *(tail + 1) = (void*) newnode;
4170  linkitems++;
4171  return (void *)(newnode + 2);
4172 }
4173 
4175 // //
4176 // insert() Inserts a node before the specified position. //
4177 // //
4178 // 'pos' (between 1 and 'linkitems') indicates the inserting position. This //
4179 // routine inserts a new node before the node of 'pos'. If 'newitem' is not //
4180 // NULL, its conents will be copied into the data slot of the new node. If //
4181 // 'pos' is larger than 'linkitems', it is equal as 'add()'. A pointer to //
4182 // the newest inserted item is returned. //
4183 // //
4185 
4186 void* tetgenmesh::link::insert(int pos, void* insitem)
4187 {
4188  if (!locate(pos)) {
4189  return add(insitem);
4190  }
4191 
4192  void **nownode = (void **) nextlinkitem;
4193 
4194  // Insert a node before 'nownode'.
4195  void **newnode = (void **) alloc();
4196  if (insitem != (void *) NULL) {
4197  memcpy((void *)(newnode + 2), insitem, linkitembytes);
4198  }
4199 
4200  *(void **)(*(nownode + 1)) = (void *) newnode;
4201  *newnode = (void *) nownode;
4202  *(newnode + 1) = *(nownode + 1);
4203  *(nownode + 1) = (void *) newnode;
4204 
4205  linkitems++;
4206 
4207  nextlinkitem = (void *) newnode;
4208  return (void *)(newnode + 2);
4209 }
4210 
4212 // //
4213 // del() Delete a node. //
4214 // //
4215 // Returns a pointer of the deleted data. If you try to delete a non-existed //
4216 // node (e.g. link is empty or a wrong index is given) return NULL. //
4217 // //
4219 
4220 void* tetgenmesh::link::deletenode(void** deadnode)
4221 {
4222  void **nextnode = (void **) *deadnode;
4223  void **prevnode = (void **) *(deadnode + 1);
4224  *prevnode = (void *) nextnode;
4225  *(nextnode + 1) = (void *) prevnode;
4226 
4227  dealloc((void *) deadnode);
4228  linkitems--;
4229 
4230  nextlinkitem = (void *) nextnode;
4231  return (void *)(deadnode + 2);
4232 }
4233 
4235 // //
4236 // del() Delete a node at the specified position. //
4237 // //
4238 // 'pos' between 1 and 'linkitems'. Returns a pointer of the deleted data. //
4239 // If you try to delete a non-existed node (e.g. link is empty or a wrong //
4240 // index is given) return NULL. //
4241 // //
4243 
4244 void* tetgenmesh::link::del(int pos)
4245 {
4246  if (!locate(pos) || (linkitems == 0)) {
4247  return (void *) NULL;
4248  }
4249  return deletenode((void **) nextlinkitem);
4250 }
4251 
4253 // //
4254 // getitem() The link traversal routine. //
4255 // //
4256 // Returns the node to which 'nextlinkitem' points. Returns a 'NULL' if the //
4257 // end of the link is reaching. Both 'nextlinkitem' and 'curpos' will be //
4258 // updated after this operation. //
4259 // //
4261 
4262 void* tetgenmesh::link::getitem()
4263 {
4264  if (nextlinkitem == (void *) tail) return NULL;
4265  void **nownode = (void **) nextlinkitem;
4266  nextlinkitem = *nownode;
4267  curpos += 1;
4268  return (void *)(nownode + 2);
4269 }
4270 
4272 // //
4273 // getnitem() Returns the node at a specified position. //
4274 // //
4275 // 'pos' between 1 and 'linkitems'. After this operation, 'nextlinkitem' and //
4276 // 'curpos' will be updated to indicate this node. //
4277 // //
4279 
4280 void* tetgenmesh::link::getnitem(int pos)
4281 {
4282  if (!locate(pos)) return NULL;
4283  return (void *)((void **) nextlinkitem + 2);
4284 }
4285 
4287 // //
4288 // hasitem() Search in this link to find if 'checkitem' exists. //
4289 // //
4290 // If 'checkitem' exists, return its position (between 1 to 'linkitems'), //
4291 // otherwise, return -1. This routine requires the linear order function has //
4292 // been set. //
4293 // //
4295 
4296 int tetgenmesh::link::hasitem(void* checkitem)
4297 {
4298  void *pathitem;
4299  int count;
4300 
4301  rewind();
4302  pathitem = getitem();
4303  count = 0;
4304  while (pathitem) {
4305  count ++;
4306  if (comp) {
4307  if ((* comp)(pathitem, checkitem) == 0) {
4308  return count;
4309  }
4310  }
4311  pathitem = getitem();
4312  }
4313  return -1;
4314 }
4315 
4316 //
4317 // End of class 'list', 'memorypool' and 'link' implementation
4318 //
4319 
4320 //
4321 // Begin of mesh manipulation primitives
4322 //
4323 
4324 //
4325 // Begin of tables initialization.
4326 //
4327 
4328 // For enumerating three edges of a triangle.
4329 
4330 int tetgenmesh::plus1mod3[3] = {1, 2, 0};
4331 int tetgenmesh::minus1mod3[3] = {2, 0, 1};
4332 
4333 // Table 've' takes an edge version as input, returns the next edge version
4334 // in the same edge ring.
4335 
4336 int tetgenmesh::ve[6] = { 2, 5, 4, 1, 0, 3 };
4337 
4338 // Tables 'vo', 'vd' and 'va' take an edge version, return the positions of
4339 // the origin, destination and apex in the triangle.
4340 
4341 int tetgenmesh::vo[6] = { 0, 1, 1, 2, 2, 0 };
4342 int tetgenmesh::vd[6] = { 1, 0, 2, 1, 0, 2 };
4343 int tetgenmesh::va[6] = { 2, 2, 0, 0, 1, 1 };
4344 
4345 // The following tables are for tetrahedron primitives (operate on trifaces).
4346 
4347 // For 'org()', 'dest()' and 'apex()'. Use 'loc' as the first index and
4348 // 'ver' as the second index.
4349 
4350 int tetgenmesh::locver2org[4][6] = {
4351  {0, 1, 1, 2, 2, 0},
4352  {0, 3, 3, 1, 1, 0},
4353  {1, 3, 3, 2, 2, 1},
4354  {2, 3, 3, 0, 0, 2}
4355 };
4356 int tetgenmesh::locver2dest[4][6] = {
4357  {1, 0, 2, 1, 0, 2},
4358  {3, 0, 1, 3, 0, 1},
4359  {3, 1, 2, 3, 1, 2},
4360  {3, 2, 0, 3, 2, 0}
4361 };
4362 int tetgenmesh::locver2apex[4][6] = {
4363  {2, 2, 0, 0, 1, 1},
4364  {1, 1, 0, 0, 3, 3},
4365  {2, 2, 1, 1, 3, 3},
4366  {0, 0, 2, 2, 3, 3}
4367 };
4368 
4369 // For oppo() primitives, use 'loc' as the index.
4370 
4371 int tetgenmesh::loc2oppo[4] = { 3, 2, 0, 1 };
4372 
4373 // For fnext() primitive. Use 'loc' as the first index and 'ver' as the
4374 // second index. Returns a new 'loc' and new 'ver' in an array. (It is
4375 // only valid for edge version equals one of {0, 2, 4}.)
4376 
4377 int tetgenmesh::locver2nextf[4][6][2] = {
4378  { {1, 5}, {-1, -1}, {2, 5}, {-1, -1}, {3, 5}, {-1, -1} },
4379  { {3, 3}, {-1, -1}, {2, 1}, {-1, -1}, {0, 1}, {-1, -1} },
4380  { {1, 3}, {-1, -1}, {3, 1}, {-1, -1}, {0, 3}, {-1, -1} },
4381  { {2, 3}, {-1, -1}, {1, 1}, {-1, -1}, {0, 5}, {-1, -1} }
4382 };
4383 
4384 // The edge number (from 0 to 5) of a tet is defined as follows:
4385 // 0 - (v0, v1), 1 - (v1, v2), 2 - (v2, v0)
4386 // 3 - (v3, v0), 4 - (v3, v1), 5 - (v3, v2).
4387 
4388 int tetgenmesh::locver2edge[4][6] = {
4389  {0, 0, 1, 1, 2, 2},
4390  {3, 3, 4, 4, 0, 0},
4391  {4, 4, 5, 5, 1, 1},
4392  {5, 5, 3, 3, 2, 2}
4393 };
4394 
4395 int tetgenmesh::edge2locver[6][2] = {
4396  {0, 0}, // 0 v0 -> v1
4397  {0, 2}, // 1 v1 -> v2
4398  {0, 4}, // 2 v2 -> v1
4399  {1, 0}, // 3 v0 -> v3
4400  {1, 2}, // 4 v1 -> v3
4401  {2, 2} // 5 v2 -> v3
4402 };
4403 
4404 //
4405 // End of tables initialization.
4406 //
4407 
4408 // Some macros for convenience
4409 
4410 #define Div2 >> 1
4411 #define Mod2 & 01
4412 
4413 // NOTE: These bit operators should only be used in macros below.
4414 
4415 // Get orient(Range from 0 to 2) from face version(Range from 0 to 5).
4416 
4417 #define Orient(V) ((V) Div2)
4418 
4419 // Determine edge ring(0 or 1) from face version(Range from 0 to 5).
4420 
4421 #define EdgeRing(V) ((V) Mod2)
4422 
4423 //
4424 // Begin of primitives for tetrahedra
4425 //
4426 
4427 // Each tetrahedron contains four pointers to its neighboring tetrahedra,
4428 // with face indices. To save memory, both information are kept in a
4429 // single pointer. To make this possible, all tetrahedra are aligned to
4430 // eight-byte boundaries, so that the last three bits of each pointer are
4431 // zeros. A face index (in the range 0 to 3) is compressed into the last
4432 // two bits of each pointer by the function 'encode()'. The function
4433 // 'decode()' decodes a pointer, extracting a face index and a pointer to
4434 // the beginning of a tetrahedron.
4435 
4436 inline void tetgenmesh::decode(tetrahedron ptr, triface& t) {
4437  t.loc = (int) ((unsigned long) (ptr) & (unsigned long) 3l);
4438  t.tet = (tetrahedron *) ((unsigned long) (ptr) & ~(unsigned long) 7l);
4439 }
4440 
4441 inline tetgenmesh::tetrahedron tetgenmesh::encode(triface& t) {
4442  return (tetrahedron) ((unsigned long) t.tet | (unsigned long) t.loc);
4443 }
4444 
4445 // sym() finds the abutting tetrahedron on the same face.
4446 
4447 inline void tetgenmesh::sym(triface& t1, triface& t2) {
4448  tetrahedron ptr = t1.tet[t1.loc];
4449  decode(ptr, t2);
4450 }
4451 
4452 inline void tetgenmesh::symself(triface& t) {
4453  tetrahedron ptr = t.tet[t.loc];
4454  decode(ptr, t);
4455 }
4456 
4457 // Bond two tetrahedra together at their faces.
4458 
4459 inline void tetgenmesh::bond(triface& t1, triface& t2) {
4460  t1.tet[t1.loc] = encode(t2);
4461  t2.tet[t2.loc] = encode(t1);
4462 }
4463 
4464 // Dissolve a bond (from one side). Note that the other tetrahedron will
4465 // still think it is connected to this tetrahedron. Usually, however,
4466 // the other tetrahedron is being deleted entirely, or bonded to another
4467 // tetrahedron, so it doesn't matter.
4468 
4469 inline void tetgenmesh::dissolve(triface& t) {
4470  t.tet[t.loc] = (tetrahedron) dummytet;
4471 }
4472 
4473 // These primitives determine or set the origin, destination, apex or
4474 // opposition of a tetrahedron with respect to 'loc' and 'ver'.
4475 
4476 inline tetgenmesh::point tetgenmesh::org(triface& t) {
4477  return (point) t.tet[locver2org[t.loc][t.ver] + 4];
4478 }
4479 
4480 inline tetgenmesh::point tetgenmesh::dest(triface& t) {
4481  return (point) t.tet[locver2dest[t.loc][t.ver] + 4];
4482 }
4483 
4484 inline tetgenmesh::point tetgenmesh::apex(triface& t) {
4485  return (point) t.tet[locver2apex[t.loc][t.ver] + 4];
4486 }
4487 
4488 inline tetgenmesh::point tetgenmesh::oppo(triface& t) {
4489  return (point) t.tet[loc2oppo[t.loc] + 4];
4490 }
4491 
4492 inline void tetgenmesh::setorg(triface& t, point pointptr) {
4493  t.tet[locver2org[t.loc][t.ver] + 4] = (tetrahedron) pointptr;
4494 }
4495 
4496 inline void tetgenmesh::setdest(triface& t, point pointptr) {
4497  t.tet[locver2dest[t.loc][t.ver] + 4] = (tetrahedron) pointptr;
4498 }
4499 
4500 inline void tetgenmesh::setapex(triface& t, point pointptr) {
4501  t.tet[locver2apex[t.loc][t.ver] + 4] = (tetrahedron) pointptr;
4502 }
4503 
4504 inline void tetgenmesh::setoppo(triface& t, point pointptr) {
4505  t.tet[loc2oppo[t.loc] + 4] = (tetrahedron) pointptr;
4506 }
4507 
4508 // These primitives were drived from Mucke's triangle-edge data structure
4509 // to change face-edge relation in a tetrahedron (esym, enext and enext2)
4510 // or between two tetrahedra (fnext).
4511 
4512 // If e0 = e(i, j), e1 = e(j, i), that is e0 and e1 are the two directions
4513 // of the same undirected edge of a face. e0.sym() = e1 and vice versa.
4514 
4515 inline void tetgenmesh::esym(triface& t1, triface& t2) {
4516  t2.tet = t1.tet;
4517  t2.loc = t1.loc;
4518  t2.ver = t1.ver + (EdgeRing(t1.ver) ? -1 : 1);
4519 }
4520 
4521 inline void tetgenmesh::esymself(triface& t) {
4522  t.ver += (EdgeRing(t.ver) ? -1 : 1);
4523 }
4524 
4525 // If e0 and e1 are both in the same edge ring of a face, e1 = e0.enext().
4526 
4527 inline void tetgenmesh::enext(triface& t1, triface& t2) {
4528  t2.tet = t1.tet;
4529  t2.loc = t1.loc;
4530  t2.ver = ve[t1.ver];
4531 }
4532 
4533 inline void tetgenmesh::enextself(triface& t) {
4534  t.ver = ve[t.ver];
4535 }
4536 
4537 // enext2() is equal to e2 = e0.enext().enext()
4538 
4539 inline void tetgenmesh::enext2(triface& t1, triface& t2) {
4540  t2.tet = t1.tet;
4541  t2.loc = t1.loc;
4542  t2.ver = ve[ve[t1.ver]];
4543 }
4544 
4545 inline void tetgenmesh::enext2self(triface& t) {
4546  t.ver = ve[ve[t.ver]];
4547 }
4548 
4549 // If f0 and f1 are both in the same face ring of a face, f1 = f0.fnext().
4550 // If f1 exists, return true. Otherwise, return false, i.e., f0 is a
4551 // boundary or hull face.
4552 
4553 inline bool tetgenmesh::fnext(triface& t1, triface& t2)
4554 {
4555  // Get the next face.
4556  t2.loc = locver2nextf[t1.loc][t1.ver][0];
4557  // Is the next face in the same tet?
4558  if (t2.loc != -1) {
4559  // It's in the same tet. Get the edge version.
4560  t2.ver = locver2nextf[t1.loc][t1.ver][1];
4561  t2.tet = t1.tet;
4562  } else {
4563  // The next face is in the neigbhour of 't1'.
4564  sym(t1, t2);
4565  if (t2.tet != dummytet) {
4566  // Find the corresponding edge in t2.
4567  point torg;
4568  int tloc, tver, i;
4569  t2.ver = 0;
4570  torg = org(t1);
4571  for (i = 0; (i < 3) && (org(t2) != torg); i++) {
4572  enextself(t2);
4573  }
4574  // Go to the next face in t2.
4575  tloc = t2.loc;
4576  tver = t2.ver;
4577  t2.loc = locver2nextf[tloc][tver][0];
4578  t2.ver = locver2nextf[tloc][tver][1];
4579  }
4580  }
4581  return t2.tet != dummytet;
4582 }
4583 
4584 inline bool tetgenmesh::fnextself(triface& t1)
4585 {
4586  triface t2;
4587 
4588  // Get the next face.
4589  t2.loc = locver2nextf[t1.loc][t1.ver][0];
4590  // Is the next face in the same tet?
4591  if (t2.loc != -1) {
4592  // It's in the same tet. Get the edge version.
4593  t2.ver = locver2nextf[t1.loc][t1.ver][1];
4594  t1.loc = t2.loc;
4595  t1.ver = t2.ver;
4596  } else {
4597  // The next face is in the neigbhour of 't1'.
4598  sym(t1, t2);
4599  if (t2.tet != dummytet) {
4600  // Find the corresponding edge in t2.
4601  point torg;
4602  int i;
4603  t2.ver = 0;
4604  torg = org(t1);
4605  for (i = 0; (i < 3) && (org(t2) != torg); i++) {
4606  enextself(t2);
4607  }
4608  t1.loc = locver2nextf[t2.loc][t2.ver][0];
4609  t1.ver = locver2nextf[t2.loc][t2.ver][1];
4610  t1.tet = t2.tet;
4611  }
4612  }
4613  return t2.tet != dummytet;
4614 }
4615 
4616 // enextfnext() and enext2fnext() are combination primitives of enext(),
4617 // enext2() and fnext().
4618 
4619 inline void tetgenmesh::enextfnext(triface& t1, triface& t2) {
4620  enext(t1, t2);
4621  fnextself(t2);
4622 }
4623 
4624 inline void tetgenmesh::enextfnextself(triface& t) {
4625  enextself(t);
4626  fnextself(t);
4627 }
4628 
4629 inline void tetgenmesh::enext2fnext(triface& t1, triface& t2) {
4630  enext2(t1, t2);
4631  fnextself(t2);
4632 }
4633 
4634 inline void tetgenmesh::enext2fnextself(triface& t) {
4635  enext2self(t);
4636  fnextself(t);
4637 }
4638 
4639 // Primitives to infect or cure a tetrahedron with the virus. The last
4640 // third bit of the pointer is marked for infection. These rely on the
4641 // assumption that all tetrahedron are aligned to eight-byte boundaries.
4642 
4643 inline void tetgenmesh::infect(triface& t) {
4644  t.tet[0] = (tetrahedron) ((unsigned long) t.tet[0] | (unsigned long) 4l);
4645 }
4646 
4647 inline void tetgenmesh::uninfect(triface& t) {
4648  t.tet[0] = (tetrahedron) ((unsigned long) t.tet[0] & ~ (unsigned long) 4l);
4649 }
4650 
4651 // Test a tetrahedron for viral infection.
4652 
4653 inline bool tetgenmesh::infected(triface& t) {
4654  return (((unsigned long) t.tet[0] & (unsigned long) 4l) != 0);
4655 }
4656 
4657 // Check or set a tetrahedron's attributes.
4658 
4659 inline REAL tetgenmesh::elemattribute(tetrahedron* ptr, int attnum) {
4660  return ((REAL *) (ptr))[elemattribindex + attnum];
4661 }
4662 
4663 inline void tetgenmesh::
4664 setelemattribute(tetrahedron* ptr, int attnum, REAL value){
4665  ((REAL *) (ptr))[elemattribindex + attnum] = value;
4666 }
4667 
4668 // Check or set a tetrahedron's maximum volume bound.
4669 
4670 inline REAL tetgenmesh::volumebound(tetrahedron* ptr) {
4671  return ((REAL *) (ptr))[volumeboundindex];
4672 }
4673 
4674 inline void tetgenmesh::setvolumebound(tetrahedron* ptr, REAL value) {
4675  ((REAL *) (ptr))[volumeboundindex] = value;
4676 }
4677 
4678 //
4679 // End of primitives for tetrahedra
4680 //
4681 
4682 //
4683 // Begin of primitives for subfaces/subsegments
4684 //
4685 
4686 // Each subface contains three pointers to its neighboring subfaces, with
4687 // edge versions. To save memory, both information are kept in a single
4688 // pointer. To make this possible, all subfaces are aligned to eight-byte
4689 // boundaries, so that the last three bits of each pointer are zeros. An
4690 // edge version (in the range 0 to 5) is compressed into the last three
4691 // bits of each pointer by 'sencode()'. 'sdecode()' decodes a pointer,
4692 // extracting an edge version and a pointer to the beginning of a subface.
4693 
4694 inline void tetgenmesh::sdecode(shellface sptr, face& s) {
4695  s.shver = (int) ((unsigned long) (sptr) & (unsigned long) 7l);
4696  s.sh = (shellface *) ((unsigned long) (sptr) & ~ (unsigned long) 7l);
4697 }
4698 
4699 inline tetgenmesh::shellface tetgenmesh::sencode(face& s) {
4700  return (shellface) ((unsigned long) s.sh | (unsigned long) s.shver);
4701 }
4702 
4703 // spivot() finds the other subface (from this subface) that shares the
4704 // same edge.
4705 
4706 inline void tetgenmesh::spivot(face& s1, face& s2) {
4707  shellface sptr = s1.sh[Orient(s1.shver)];
4708  sdecode(sptr, s2);
4709 }
4710 
4711 inline void tetgenmesh::spivotself(face& s) {
4712  shellface sptr = s.sh[Orient(s.shver)];
4713  sdecode(sptr, s);
4714 }
4715 
4716 // sbond() bonds two subfaces together, i.e., after bonding, both faces
4717 // are pointing to each other.
4718 
4719 inline void tetgenmesh::sbond(face& s1, face& s2) {
4720  s1.sh[Orient(s1.shver)] = sencode(s2);
4721  s2.sh[Orient(s2.shver)] = sencode(s1);
4722 }
4723 
4724 // sbond1() only bonds s2 to s1, i.e., after bonding, s1 is pointing to s2,
4725 // but s2 is not pointing to s1.
4726 
4727 inline void tetgenmesh::sbond1(face& s1, face& s2) {
4728  s1.sh[Orient(s1.shver)] = sencode(s2);
4729 }
4730 
4731 // Dissolve a subface bond (from one side). Note that the other subface
4732 // will still think it's connected to this subface.
4733 
4734 inline void tetgenmesh::sdissolve(face& s) {
4735  s.sh[Orient(s.shver)] = (shellface) dummysh;
4736 }
4737 
4738 // These primitives determine or set the origin, destination, or apex
4739 // of a subface with respect to the edge version.
4740 
4741 inline tetgenmesh::point tetgenmesh::sorg(face& s) {
4742  return (point) s.sh[3 + vo[s.shver]];
4743 }
4744 
4745 inline tetgenmesh::point tetgenmesh::sdest(face& s) {
4746  return (point) s.sh[3 + vd[s.shver]];
4747 }
4748 
4749 inline tetgenmesh::point tetgenmesh::sapex(face& s) {
4750  return (point) s.sh[3 + va[s.shver]];
4751 }
4752 
4753 inline void tetgenmesh::setsorg(face& s, point pointptr) {
4754  s.sh[3 + vo[s.shver]] = (shellface) pointptr;
4755 }
4756 
4757 inline void tetgenmesh::setsdest(face& s, point pointptr) {
4758  s.sh[3 + vd[s.shver]] = (shellface) pointptr;
4759 }
4760 
4761 inline void tetgenmesh::setsapex(face& s, point pointptr) {
4762  s.sh[3 + va[s.shver]] = (shellface) pointptr;
4763 }
4764 
4765 // These primitives were drived from Mucke[2]'s triangle-edge data structure
4766 // to change face-edge relation in a subface (sesym, senext and senext2).
4767 
4768 inline void tetgenmesh::sesym(face& s1, face& s2) {
4769  s2.sh = s1.sh;
4770  s2.shver = s1.shver + (EdgeRing(s1.shver) ? -1 : 1);
4771 }
4772 
4773 inline void tetgenmesh::sesymself(face& s) {
4774  s.shver += (EdgeRing(s.shver) ? -1 : 1);
4775 }
4776 
4777 inline void tetgenmesh::senext(face& s1, face& s2) {
4778  s2.sh = s1.sh;
4779  s2.shver = ve[s1.shver];
4780 }
4781 
4782 inline void tetgenmesh::senextself(face& s) {
4783  s.shver = ve[s.shver];
4784 }
4785 
4786 inline void tetgenmesh::senext2(face& s1, face& s2) {
4787  s2.sh = s1.sh;
4788  s2.shver = ve[ve[s1.shver]];
4789 }
4790 
4791 inline void tetgenmesh::senext2self(face& s) {
4792  s.shver = ve[ve[s.shver]];
4793 }
4794 
4795 // If f0 and f1 are both in the same face ring, then f1 = f0.fnext(),
4796 
4797 inline void tetgenmesh::sfnext(face& s1, face& s2) {
4798  getnextsface(&s1, &s2);
4799 }
4800 
4801 inline void tetgenmesh::sfnextself(face& s) {
4802  getnextsface(&s, NULL);
4803 }
4804 
4805 // These primitives read or set a pointer of the badface structure. The
4806 // pointer is stored sh[11].
4807 
4808 inline tetgenmesh::badface* tetgenmesh::shell2badface(face& s) {
4809  return (badface*) s.sh[11];
4810 }
4811 
4812 inline void tetgenmesh::setshell2badface(face& s, badface* value) {
4813  s.sh[11] = (shellface) value;
4814 }
4815 
4816 // Check or set a subface's maximum area bound.
4817 
4818 inline REAL tetgenmesh::areabound(face& s) {
4819  return ((REAL *) (s.sh))[areaboundindex];
4820 }
4821 
4822 inline void tetgenmesh::setareabound(face& s, REAL value) {
4823  ((REAL *) (s.sh))[areaboundindex] = value;
4824 }
4825 
4826 // These two primitives read or set a shell marker. Shell markers are used
4827 // to hold user boundary information.
4828 
4829 inline int tetgenmesh::shellmark(face& s) {
4830  return ((int *) (s.sh))[shmarkindex];
4831 }
4832 
4833 inline void tetgenmesh::setshellmark(face& s, int value) {
4834  ((int *) (s.sh))[shmarkindex] = value;
4835 }
4836 
4837 // These two primitives set or read the type of the subface or subsegment.
4838 
4839 inline enum tetgenmesh::shestype tetgenmesh::shelltype(face& s) {
4840  return (enum shestype) ((int *) (s.sh))[shmarkindex + 1];
4841 }
4842 
4843 inline void tetgenmesh::setshelltype(face& s, enum shestype value) {
4844  ((int *) (s.sh))[shmarkindex + 1] = (int) value;
4845 }
4846 
4847 // These two primitives set or read the pbc group of the subface.
4848 
4849 inline int tetgenmesh::shellpbcgroup(face& s) {
4850  return ((int *) (s.sh))[shmarkindex + 2];
4851 }
4852 
4853 inline void tetgenmesh::setshellpbcgroup(face& s, int value) {
4854  ((int *) (s.sh))[shmarkindex + 2] = value;
4855 }
4856 
4857 // Primitives to infect or cure a subface with the virus. These rely on the
4858 // assumption that all tetrahedra are aligned to eight-byte boundaries.
4859 
4860 inline void tetgenmesh::sinfect(face& s) {
4861  s.sh[6] = (shellface) ((unsigned long) s.sh[6] | (unsigned long) 4l);
4862 }
4863 
4864 inline void tetgenmesh::suninfect(face& s) {
4865  s.sh[6] = (shellface)((unsigned long) s.sh[6] & ~(unsigned long) 4l);
4866 }
4867 
4868 // Test a subface for viral infection.
4869 
4870 inline bool tetgenmesh::sinfected(face& s) {
4871  return (((unsigned long) s.sh[6] & (unsigned long) 4l) != 0);
4872 }
4873 
4874 //
4875 // End of primitives for subfaces/subsegments
4876 //
4877 
4878 //
4879 // Begin of primitives for interacting between tetrahedra and subfaces
4880 //
4881 
4882 // tspivot() finds a subface abutting on this tetrahdera.
4883 
4884 inline void tetgenmesh::tspivot(triface& t, face& s) {
4885  shellface sptr = (shellface) t.tet[8 + t.loc];
4886  sdecode(sptr, s);
4887 }
4888 
4889 // stpivot() finds a tetrahedron abutting a subface.
4890 
4891 inline void tetgenmesh::stpivot(face& s, triface& t) {
4892  tetrahedron ptr = (tetrahedron) s.sh[6 + EdgeRing(s.shver)];
4893  decode(ptr, t);
4894 }
4895 
4896 // tsbond() bond a tetrahedron to a subface.
4897 
4898 inline void tetgenmesh::tsbond(triface& t, face& s) {
4899  t.tet[8 + t.loc] = (tetrahedron) sencode(s);
4900  s.sh[6 + EdgeRing(s.shver)] = (shellface) encode(t);
4901 }
4902 
4903 // tsdissolve() dissolve a bond (from the tetrahedron side).
4904 
4905 inline void tetgenmesh::tsdissolve(triface& t) {
4906  t.tet[8 + t.loc] = (tetrahedron) dummysh;
4907 }
4908 
4909 // stdissolve() dissolve a bond (from the subface side).
4910 
4911 inline void tetgenmesh::stdissolve(face& s) {
4912  s.sh[6 + EdgeRing(s.shver)] = (shellface) dummytet;
4913 }
4914 
4915 //
4916 // End of primitives for interacting between tetrahedra and subfaces
4917 //
4918 
4919 //
4920 // Begin of primitives for interacting between subfaces and subsegs
4921 //
4922 
4923 // sspivot() finds a subsegment abutting a subface.
4924 
4925 inline void tetgenmesh::sspivot(face& s, face& edge) {
4926  shellface sptr = (shellface) s.sh[8 + Orient(s.shver)];
4927  sdecode(sptr, edge);
4928 }
4929 
4930 // ssbond() bond a subface to a subsegment.
4931 
4932 inline void tetgenmesh::ssbond(face& s, face& edge) {
4933  s.sh[8 + Orient(s.shver)] = sencode(edge);
4934  edge.sh[0] = sencode(s);
4935 }
4936 
4937 // ssdisolve() dissolve a bond (from the subface side)
4938 
4939 inline void tetgenmesh::ssdissolve(face& s) {
4940  s.sh[8 + Orient(s.shver)] = (shellface) dummysh;
4941 }
4942 
4943 //
4944 // End of primitives for interacting between subfaces and subsegs
4945 //
4946 
4947 //
4948 // Begin of primitives for interacting between tet and subsegs.
4949 //
4950 
4951 inline void tetgenmesh::tsspivot1(triface& t, face& seg)
4952 {
4953  shellface sptr = (shellface) t.tet[8 + locver2edge[t.loc][t.ver]];
4954  sdecode(sptr, seg);
4955 }
4956 
4957 // Only bond/dissolve at tet's side, but not vice versa.
4958 
4959 inline void tetgenmesh::tssbond1(triface& t, face& seg)
4960 {
4961  t.tet[8 + locver2edge[t.loc][t.ver]] = (tetrahedron) sencode(seg);
4962 }
4963 
4964 inline void tetgenmesh::tssdissolve1(triface& t)
4965 {
4966  t.tet[8 + locver2edge[t.loc][t.ver]] = (tetrahedron) dummysh;
4967 }
4968 
4969 //
4970 // End of primitives for interacting between tet and subsegs.
4971 //
4972 
4973 //
4974 // Begin of primitives for points
4975 //
4976 
4977 inline int tetgenmesh::pointmark(point pt) {
4978  return ((int *) (pt))[pointmarkindex];
4979 }
4980 
4981 inline void tetgenmesh::setpointmark(point pt, int value) {
4982  ((int *) (pt))[pointmarkindex] = value;
4983 }
4984 
4985 // These two primitives set and read the type of the point.
4986 
4987 inline enum tetgenmesh::verttype tetgenmesh::pointtype(point pt) {
4988  return (enum verttype) ((int *) (pt))[pointmarkindex + 1];
4989 }
4990 
4991 inline void tetgenmesh::setpointtype(point pt, enum verttype value) {
4992  ((int *) (pt))[pointmarkindex + 1] = (int) value;
4993 }
4994 
4995 // These following primitives set and read a pointer to a tetrahedron
4996 // a subface/subsegment, a point, or a tet of background mesh.
4997 
4998 inline tetgenmesh::tetrahedron tetgenmesh::point2tet(point pt) {
4999  return ((tetrahedron *) (pt))[point2simindex];
5000 }
5001 
5002 inline void tetgenmesh::setpoint2tet(point pt, tetrahedron value) {
5003  ((tetrahedron *) (pt))[point2simindex] = value;
5004 }
5005 
5006 inline tetgenmesh::shellface tetgenmesh::point2sh(point pt) {
5007  return (shellface) ((tetrahedron *) (pt))[point2simindex + 1];
5008 }
5009 
5010 inline void tetgenmesh::setpoint2sh(point pt, shellface value) {
5011  ((tetrahedron *) (pt))[point2simindex + 1] = (tetrahedron) value;
5012 }
5013 
5014 inline tetgenmesh::point tetgenmesh::point2ppt(point pt) {
5015  return (point) ((tetrahedron *) (pt))[point2simindex + 2];
5016 }
5017 
5018 inline void tetgenmesh::setpoint2ppt(point pt, point value) {
5019  ((tetrahedron *) (pt))[point2simindex + 2] = (tetrahedron) value;
5020 }
5021 
5022 inline tetgenmesh::tetrahedron tetgenmesh::point2bgmtet(point pt) {
5023  return ((tetrahedron *) (pt))[point2simindex + 3];
5024 }
5025 
5026 inline void tetgenmesh::setpoint2bgmtet(point pt, tetrahedron value) {
5027  ((tetrahedron *) (pt))[point2simindex + 3] = value;
5028 }
5029 
5030 // These primitives set and read a pointer to its pbc point.
5031 
5032 inline tetgenmesh::point tetgenmesh::point2pbcpt(point pt) {
5033  return (point) ((tetrahedron *) (pt))[point2pbcptindex];
5034 }
5035 
5036 inline void tetgenmesh::setpoint2pbcpt(point pt, point value) {
5037  ((tetrahedron *) (pt))[point2pbcptindex] = (tetrahedron) value;
5038 }
5039 
5040 //
5041 // End of primitives for points
5042 //
5043 
5044 //
5045 // Begin of advanced primitives
5046 //
5047 
5048 // adjustedgering() adjusts the edge version so that it belongs to the
5049 // indicated edge ring. The 'direction' only can be 0(CCW) or 1(CW).
5050 // If the edge is not in the wanted edge ring, reverse it.
5051 
5052 inline void tetgenmesh::adjustedgering(triface& t, int direction) {
5053  if (EdgeRing(t.ver) != direction) {
5054  esymself(t);
5055  }
5056 }
5057 
5058 inline void tetgenmesh::adjustedgering(face& s, int direction) {
5059  if (EdgeRing(s.shver) != direction) {
5060  sesymself(s);
5061  }
5062 }
5063 
5064 // isdead() returns TRUE if the tetrahedron or subface has been dealloced.
5065 
5066 inline bool tetgenmesh::isdead(triface* t) {
5067  if (t->tet == (tetrahedron *) NULL) return true;
5068  else return t->tet[4] == (tetrahedron) NULL;
5069 }
5070 
5071 inline bool tetgenmesh::isdead(face* s) {
5072  if (s->sh == (shellface *) NULL) return true;
5073  else return s->sh[3] == (shellface) NULL;
5074 }
5075 
5076 // isfacehaspoint() returns TRUE if the 'testpoint' is one of the vertices
5077 // of the tetface 't' subface 's'.
5078 
5079 inline bool tetgenmesh::isfacehaspoint(triface* t, point testpoint) {
5080  return ((org(*t) == testpoint) || (dest(*t) == testpoint) ||
5081  (apex(*t) == testpoint));
5082 }
5083 
5084 inline bool tetgenmesh::isfacehaspoint(face* s, point testpoint) {
5085  return (s->sh[3] == (shellface) testpoint) ||
5086  (s->sh[4] == (shellface) testpoint) ||
5087  (s->sh[5] == (shellface) testpoint);
5088 }
5089 
5090 // isfacehasedge() returns TRUE if the edge (given by its two endpoints) is
5091 // one of the three edges of the subface 's'.
5092 
5093 inline bool tetgenmesh::isfacehasedge(face* s, point tend1, point tend2) {
5094  return (isfacehaspoint(s, tend1) && isfacehaspoint(s, tend2));
5095 }
5096 
5097 // issymexist() returns TRUE if the adjoining tetrahedron is not 'duumytet'.
5098 
5099 inline bool tetgenmesh::issymexist(triface* t) {
5100  tetrahedron *ptr = (tetrahedron *)
5101  ((unsigned long)(t->tet[t->loc]) & ~(unsigned long)7l);
5102  return ptr != dummytet;
5103 }
5104 
5106 // //
5107 // getnextsface() Finds the next subface in the face ring. //
5108 // //
5109 // For saving space in the data structure of subface, there only exists one //
5110 // face ring around a segment (see programming manual). This routine imple- //
5111 // ments the double face ring as desired in Muecke's data structure. //
5112 // //
5114 
5115 void tetgenmesh::getnextsface(face* s1, face* s2)
5116 {
5117  face neighsh, spinsh;
5118  face testseg;
5119 
5120  sspivot(*s1, testseg);
5121  if (testseg.sh != dummysh) {
5122  testseg.shver = 0;
5123  if (sorg(testseg) == sorg(*s1)) {
5124  spivot(*s1, neighsh);
5125  } else {
5126  spinsh = *s1;
5127  do {
5128  neighsh = spinsh;
5129  spivotself(spinsh);
5130  } while (spinsh.sh != s1->sh);
5131  }
5132  } else {
5133  spivot(*s1, neighsh);
5134  }
5135  if (sorg(neighsh) != sorg(*s1)) {
5136  sesymself(neighsh);
5137  }
5138  if (s2 != (face *) NULL) {
5139  *s2 = neighsh;
5140  } else {
5141  *s1 = neighsh;
5142  }
5143 }
5144 
5146 // //
5147 // tsspivot() Finds a subsegment abutting on a tetrahderon's edge. //
5148 // //
5149 // The edge is represented in the primary edge of 'checkedge'. If there is a //
5150 // subsegment bonded at this edge, it is returned in handle 'checkseg', the //
5151 // edge direction of 'checkseg' is conformed to 'checkedge'. If there isn't, //
5152 // set 'checkseg.sh = dummysh' to indicate it is not a subsegment. //
5153 // //
5154 // To find whether an edge of a tetrahedron is a subsegment or not. First we //
5155 // need find a subface around this edge to see if it contains a subsegment. //
5156 // The reason is there is no direct connection between a tetrahedron and its //
5157 // adjoining subsegments. //
5158 // //
5160 
5161 void tetgenmesh::tsspivot(triface* checkedge, face* checkseg)
5162 {
5163  triface spintet;
5164  face parentsh;
5165  point tapex;
5166  int hitbdry;
5167 
5168  spintet = *checkedge;
5169  tapex = apex(*checkedge);
5170  hitbdry = 0;
5171  do {
5172  tspivot(spintet, parentsh);
5173  // Does spintet have a (non-fake) subface attached?
5174  if ((parentsh.sh != dummysh) && (sapex(parentsh) != NULL)) {
5175  // Find a subface! Find the edge in it.
5176  findedge(&parentsh, org(*checkedge), dest(*checkedge));
5177  sspivot(parentsh, *checkseg);
5178  if (checkseg->sh != dummysh) {
5179  // Find a subsegment! Correct its edge direction before return.
5180  if (sorg(*checkseg) != org(*checkedge)) {
5181  sesymself(*checkseg);
5182  }
5183  }
5184  return;
5185  }
5186  if (!fnextself(spintet)) {
5187  hitbdry++;
5188  if (hitbdry < 2) {
5189  esym(*checkedge, spintet);
5190  if (!fnextself(spintet)) {
5191  hitbdry++;
5192  }
5193  }
5194  }
5195  } while ((apex(spintet) != tapex) && (hitbdry < 2));
5196  // Not find.
5197  checkseg->sh = dummysh;
5198 }
5199 
5201 // //
5202 // sstpivot() Finds a tetrahedron abutting a subsegment. //
5203 // //
5204 // This is the inverse operation of 'tsspivot()'. One subsegment shared by //
5205 // arbitrary number of tetrahedron, the returned tetrahedron is not unique. //
5206 // The edge direction of the returned tetrahedron is conformed to the given //
5207 // subsegment. //
5208 // //
5210 
5211 void tetgenmesh::sstpivot(face* checkseg, triface* retedge)
5212 {
5213  face parentsh;
5214 
5215  // Get the subface which holds the subsegment.
5216  sdecode(checkseg->sh[0], parentsh);
5217 #ifdef SELF_CHECK
5218  assert(parentsh.sh != dummysh);
5219 #endif
5220  // Get a tetraheron to which the subface attches.
5221  stpivot(parentsh, *retedge);
5222  if (retedge->tet == dummytet) {
5223  sesymself(parentsh);
5224  stpivot(parentsh, *retedge);
5225 #ifdef SELF_CHECK
5226  assert(retedge->tet != dummytet);
5227 #endif
5228  }
5229  // Correct the edge direction before return.
5230  findedge(retedge, sorg(*checkseg), sdest(*checkseg));
5231 }
5232 
5234 // //
5235 // findorg() Finds a point in the given handle (tetrahedron or subface). //
5236 // //
5237 // If 'dorg' is a one of vertices of the given handle, set the origin of //
5238 // this handle be that point and return TRUE. Otherwise, return FALSE and //
5239 // 'tface' remains unchanged. //
5240 // //
5242 
5243 bool tetgenmesh::findorg(triface* tface, point dorg)
5244 {
5245  if (org(*tface) == dorg) {
5246  return true;
5247  } else {
5248  if (dest(*tface) == dorg) {
5249  enextself(*tface);
5250  return true;
5251  } else {
5252  if (apex(*tface) == dorg) {
5253  enext2self(*tface);
5254  return true;
5255  } else {
5256  if (oppo(*tface) == dorg) {
5257  // Keep 'tface' referring to the same tet after fnext().
5258  adjustedgering(*tface, CCW);
5259  fnextself(*tface);
5260  enext2self(*tface);
5261  return true;
5262  }
5263  }
5264  }
5265  }
5266  return false;
5267 }
5268 
5269 bool tetgenmesh::findorg(face* sface, point dorg)
5270 {
5271  if (sorg(*sface) == dorg) {
5272  return true;
5273  } else {
5274  if (sdest(*sface) == dorg) {
5275  senextself(*sface);
5276  return true;
5277  } else {
5278  if (sapex(*sface) == dorg) {
5279  senext2self(*sface);
5280  return true;
5281  }
5282  }
5283  }
5284  return false;
5285 }
5286 
5288 // //
5289 // findedge() Find an edge in the given handle (tetrahedron or subface). //
5290 // //
5291 // The edge is given in two points 'eorg' and 'edest'. It is assumed that //
5292 // the edge must exist in the given handle (tetrahedron or subface). This //
5293 // routine sets the right edge version for the input handle. //
5294 // //
5296 
5297 void tetgenmesh::findedge(triface* tface, point eorg, point edest)
5298 {
5299  int i;
5300 
5301  for (i = 0; i < 3; i++) {
5302  if (org(*tface) == eorg) {
5303  if (dest(*tface) == edest) {
5304  // Edge is found, return.
5305  return;
5306  }
5307  } else {
5308  if (org(*tface) == edest) {
5309  if (dest(*tface) == eorg) {
5310  // Edge is found, inverse the direction and return.
5311  esymself(*tface);
5312  return;
5313  }
5314  }
5315  }
5316  enextself(*tface);
5317  }
5318  // It should never be here.
5319  printf("Internalerror in findedge(): Unable to find an edge in tet.\n");
5320  internalerror();
5321 }
5322 
5323 void tetgenmesh::findedge(face* sface, point eorg, point edest)
5324 {
5325  int i;
5326 
5327  for (i = 0; i < 3; i++) {
5328  if (sorg(*sface) == eorg) {
5329  if (sdest(*sface) == edest) {
5330  // Edge is found, return.
5331  return;
5332  }
5333  } else {
5334  if (sorg(*sface) == edest) {
5335  if (sdest(*sface) == eorg) {
5336  // Edge is found, inverse the direction and return.
5337  sesymself(*sface);
5338  return;
5339  }
5340  }
5341  }
5342  senextself(*sface);
5343  }
5344  printf("Internalerror in findedge(): Unable to find an edge in subface.\n");
5345  internalerror();
5346 }
5347 
5349 // //
5350 // findface() Find the face has the given origin, destination and apex. //
5351 // //
5352 // On input, 'fface' is a handle which may contain the three corners or may //
5353 // not or may be dead. On return, it represents exactly the face with the //
5354 // given origin, destination and apex. //
5355 // //
5357 
5358 void tetgenmesh::findface(triface *fface, point forg, point fdest, point fapex)
5359 {
5360  triface spintet;
5361  enum finddirectionresult collinear;
5362  int hitbdry;
5363 
5364  if (!isdead(fface)) {
5365  // First check the easiest case, that 'fface' is just the right one.
5366  if (org(*fface) == forg && dest(*fface) == fdest &&
5367  apex(*fface) == fapex) return;
5368  } else {
5369  // The input handle is dead, use the 'recenttet' if it is alive.
5370  if (!isdead(&recenttet)) *fface = recenttet;
5371  }
5372 
5373  if (!isdead(fface)) {
5374  if (!findorg(fface, forg)) {
5375  // 'forg' is not a corner of 'fface', locate it.
5376  preciselocate(forg, fface, tetrahedrons->items);
5377  }
5378  // It is possible that forg is not found in a non-convex mesh.
5379  if (org(*fface) == forg) {
5380  collinear = finddirection(fface, fdest, tetrahedrons->items);
5381  if (collinear == RIGHTCOLLINEAR) {
5382  // fdest is just the destination.
5383  } else if (collinear == LEFTCOLLINEAR) {
5384  enext2self(*fface);
5385  esymself(*fface);
5386  } else if (collinear == TOPCOLLINEAR) {
5387  fnextself(*fface);
5388  enext2self(*fface);
5389  esymself(*fface);
5390  }
5391  }
5392  // It is possible taht fdest is not found in a non-convex mesh.
5393  if ((org(*fface) == forg) && (dest(*fface) == fdest)) {
5394  // Find the apex of 'fapex'.
5395  spintet = *fface;
5396  hitbdry = 0;
5397  do {
5398  if (apex(spintet) == fapex) {
5399  // We have done. Be careful the edge direction of 'spintet',
5400  // it may reversed because of hitting boundary once.
5401  if (org(spintet) != org(*fface)) {
5402  esymself(spintet);
5403  }
5404  *fface = spintet;
5405  return;
5406  }
5407  if (!fnextself(spintet)) {
5408  hitbdry ++;
5409  if (hitbdry < 2) {
5410  esym(*fface, spintet);
5411  if (!fnextself(spintet)) {
5412  hitbdry ++;
5413  }
5414  }
5415  }
5416  } while (hitbdry < 2 && apex(spintet) != apex(*fface));
5417  // It is possible that fapex is not found in a non-convex mesh.
5418  }
5419  }
5420 
5421  if (isdead(fface) || (org(*fface) != forg) || (dest(*fface) != fdest) ||
5422  (apex(*fface) != fapex)) {
5423  // Too bad, the input handle is useless. We have to find a handle
5424  // for 'fface' contains the 'forg' and 'fdest'. Here a brute force
5425  // search is performed.
5426  if (b->verbose > 1) {
5427  printf("Warning in findface(): Perform a brute-force searching.\n");
5428  }
5429  enum verttype forgty, fdestty, fapexty;
5430  int share, i;
5431  forgty = pointtype(forg);
5432  fdestty = pointtype(fdest);
5433  fapexty = pointtype(fapex);
5434  setpointtype(forg, DEADVERTEX);
5435  setpointtype(fdest, DEADVERTEX);
5436  setpointtype(fapex, DEADVERTEX);
5437  tetrahedrons->traversalinit();
5438  fface->tet = tetrahedrontraverse();
5439  while (fface->tet != (tetrahedron *) NULL) {
5440  share = 0;
5441  for (i = 0; i < 4; i++) {
5442  if (pointtype((point) fface->tet[4 + i]) == DEADVERTEX) share ++;
5443  }
5444  if (share == 3) {
5445  // Found! Set the correct face and desired corners.
5446  if (pointtype((point) fface->tet[4]) != DEADVERTEX) {
5447  fface->loc = 2;
5448  } else if (pointtype((point) fface->tet[5]) != DEADVERTEX) {
5449  fface->loc = 3;
5450  } else if (pointtype((point) fface->tet[6]) != DEADVERTEX) {
5451  fface->loc = 1;
5452  } else { // pointtype((point) fface->tet[7]) != DEADVERTEX
5453  fface->loc = 0;
5454  }
5455  findedge(fface, forg, fdest);
5456  break;
5457  }
5458  fface->tet = tetrahedrontraverse();
5459  }
5460  setpointtype(forg, forgty);
5461  setpointtype(fdest, fdestty);
5462  setpointtype(fapex, fapexty);
5463  if (fface->tet == (tetrahedron *) NULL) {
5464  // It is impossible to reach here.
5465  printf("Internal error: Fail to find the indicated face.\n");
5466  internalerror();
5467  }
5468  }
5469 }
5470 
5472 // //
5473 // getonextseg() Get the next SEGMENT counterclockwise with the same org. //
5474 // //
5475 // 's' is a subface. This routine reteuns the segment which is counterclock- //
5476 // wise with the origin of s. //
5477 // //
5479 
5480 void tetgenmesh::getonextseg(face* s, face* lseg)
5481 {
5482  face checksh, checkseg;
5483  point forg;
5484 
5485  forg = sorg(*s);
5486  checksh = *s;
5487  do {
5488  // Go to the edge at forg's left side.
5489  senext2self(checksh);
5490  // Check if there is a segment attaching this edge.
5491  sspivot(checksh, checkseg);
5492  if (checkseg.sh != dummysh) break;
5493  // No segment! Go to the neighbor of this subface.
5494  spivotself(checksh);
5495 #ifdef SELF_CHECK
5496  // It should always meet a segment before come back.
5497  assert(checksh.sh != s->sh);
5498 #endif
5499  if (sorg(checksh) != forg) {
5500  sesymself(checksh);
5501 #ifdef SELF_CHECK
5502  assert(sorg(checksh) == forg);
5503 #endif
5504  }
5505  } while (true);
5506  if (sorg(checkseg) != forg) sesymself(checkseg);
5507  *lseg = checkseg;
5508 }
5509 
5511 // //
5512 // getseghasorg() Get the segment containing the given point. //
5513 // //
5514 // 'dorg' is an endpoint of a segment S. 'sseg' is a subsegment of S. This //
5515 // routine search a subsegment (along sseg) of S containing dorg. On return, //
5516 // 'sseg' contains 'dorg' as its origin. //
5517 // //
5519 
5520 void tetgenmesh::getseghasorg(face* sseg, point dorg)
5521 {
5522  face nextseg;
5523  point checkpt;
5524 
5525  nextseg = *sseg;
5526  checkpt = sorg(nextseg);
5527  while ((checkpt != dorg) && (pointtype(checkpt) == FREESEGVERTEX)) {
5528  // Search dorg along the original direction of sseg.
5529  senext2self(nextseg);
5530  spivotself(nextseg);
5531  nextseg.shver = 0;
5532  if (sdest(nextseg) != checkpt) sesymself(nextseg);
5533  checkpt = sorg(nextseg);
5534  }
5535  if (checkpt == dorg) {
5536  *sseg = nextseg;
5537  return;
5538  }
5539  nextseg = *sseg;
5540  checkpt = sdest(nextseg);
5541  while ((checkpt != dorg) && (pointtype(checkpt) == FREESEGVERTEX)) {
5542  // Search dorg along the destinational direction of sseg.
5543  senextself(nextseg);
5544  spivotself(nextseg);
5545  nextseg.shver = 0;
5546  if (sorg(nextseg) != checkpt) sesymself(nextseg);
5547  checkpt = sdest(nextseg);
5548  }
5549  if (checkpt == dorg) {
5550  sesym(nextseg, *sseg);
5551  return;
5552  }
5553  // Should never be here.
5554  printf("Internalerror in getseghasorg(): Unable to find the subseg.\n");
5555  internalerror();
5556 }
5557 
5559 // //
5560 // getsubsegfarorg() Get the origin of the parent segment of a subseg. //
5561 // //
5563 
5564 tetgenmesh::point tetgenmesh::getsubsegfarorg(face* sseg)
5565 {
5566  face prevseg;
5567  point checkpt;
5568 
5569  checkpt = sorg(*sseg);
5570  senext2(*sseg, prevseg);
5571  spivotself(prevseg);
5572  // Search dorg along the original direction of sseg.
5573  while (prevseg.sh != dummysh) {
5574  prevseg.shver = 0;
5575  if (sdest(prevseg) != checkpt) sesymself(prevseg);
5576  checkpt = sorg(prevseg);
5577  senext2self(prevseg);
5578  spivotself(prevseg);
5579  }
5580  return checkpt;
5581 }
5582 
5584 // //
5585 // getsubsegfardest() Get the dest. of the parent segment of a subseg. //
5586 // //
5588 
5589 tetgenmesh::point tetgenmesh::getsubsegfardest(face* sseg)
5590 {
5591  face nextseg;
5592  point checkpt;
5593 
5594  checkpt = sdest(*sseg);
5595  senext(*sseg, nextseg);
5596  spivotself(nextseg);
5597  // Search dorg along the destinational direction of sseg.
5598  while (nextseg.sh != dummysh) {
5599  nextseg.shver = 0;
5600  if (sorg(nextseg) != checkpt) sesymself(nextseg);
5601  checkpt = sdest(nextseg);
5602  senextself(nextseg);
5603  spivotself(nextseg);
5604  }
5605  return checkpt;
5606 }
5607 
5609 // //
5610 // printtet() Print out the details of a tetrahedron on screen. //
5611 // //
5612 // It's also used when the highest level of verbosity (`-VVV') is specified. //
5613 // //
5615 
5616 void tetgenmesh::printtet(triface* tface)
5617 {
5618  triface tmpface, prtface;
5619  point tmppt;
5620  face tmpsh;
5621  int facecount;
5622 
5623  printf("Tetra x%lx with loc(%i) and ver(%i):",
5624  (unsigned long)(tface->tet), tface->loc, tface->ver);
5625  if (infected(*tface)) {
5626  printf(" (infected)");
5627  }
5628  printf("\n");
5629 
5630  tmpface = *tface;
5631  facecount = 0;
5632  while(facecount < 4) {
5633  tmpface.loc = facecount;
5634  sym(tmpface, prtface);
5635  if(prtface.tet == dummytet) {
5636  printf(" [%i] Outer space.\n", facecount);
5637  } else {
5638  printf(" [%i] x%lx loc(%i).", facecount,
5639  (unsigned long)(prtface.tet), prtface.loc);
5640  if (infected(prtface)) {
5641  printf(" (infected)");
5642  }
5643  printf("\n");
5644  }
5645  facecount ++;
5646  }
5647 
5648  tmppt = org(*tface);
5649  if(tmppt == (point) NULL) {
5650  printf(" Org [%i] NULL\n", locver2org[tface->loc][tface->ver]);
5651  } else {
5652  printf(" Org [%i] x%lx (%.12g,%.12g,%.12g) %d\n",
5653  locver2org[tface->loc][tface->ver], (unsigned long)(tmppt),
5654  tmppt[0], tmppt[1], tmppt[2], pointmark(tmppt));
5655  }
5656  tmppt = dest(*tface);
5657  if(tmppt == (point) NULL) {
5658  printf(" Dest[%i] NULL\n", locver2dest[tface->loc][tface->ver]);
5659  } else {
5660  printf(" Dest[%i] x%lx (%.12g,%.12g,%.12g) %d\n",
5661  locver2dest[tface->loc][tface->ver], (unsigned long)(tmppt),
5662  tmppt[0], tmppt[1], tmppt[2], pointmark(tmppt));
5663  }
5664  tmppt = apex(*tface);
5665  if(tmppt == (point) NULL) {
5666  printf(" Apex[%i] NULL\n", locver2apex[tface->loc][tface->ver]);
5667  } else {
5668  printf(" Apex[%i] x%lx (%.12g,%.12g,%.12g) %d\n",
5669  locver2apex[tface->loc][tface->ver], (unsigned long)(tmppt),
5670  tmppt[0], tmppt[1], tmppt[2], pointmark(tmppt));
5671  }
5672  tmppt = oppo(*tface);
5673  if(tmppt == (point) NULL) {
5674  printf(" Oppo[%i] NULL\n", loc2oppo[tface->loc]);
5675  } else {
5676  printf(" Oppo[%i] x%lx (%.12g,%.12g,%.12g) %d\n",
5677  loc2oppo[tface->loc], (unsigned long)(tmppt),
5678  tmppt[0], tmppt[1], tmppt[2], pointmark(tmppt));
5679  }
5680 
5681  if (b->useshelles) {
5682  tmpface = *tface;
5683  facecount = 0;
5684  while(facecount < 6) {
5685  tmpface.loc = facecount;
5686  tspivot(tmpface, tmpsh);
5687  if(tmpsh.sh != dummysh) {
5688  printf(" [%i] x%lx ID(%i) ", facecount,
5689  (unsigned long)(tmpsh.sh), shellmark(tmpsh));
5690  if (sorg(tmpsh) == (point) NULL) {
5691  printf("(fake)");
5692  }
5693  printf("\n");
5694  }
5695  facecount ++;
5696  }
5697  }
5698 }
5699 
5701 // //
5702 // printsh() Print out the details of a subface or subsegment on screen. //
5703 // //
5704 // It's also used when the highest level of verbosity (`-VVV') is specified. //
5705 // //
5707 
5708 void tetgenmesh::printsh(face* sface)
5709 {
5710  face prtsh;
5711  triface prttet;
5712  point printpoint;
5713 
5714  if (sapex(*sface) != NULL) {
5715  printf("subface x%lx, ver %d, mark %d:",
5716  (unsigned long)(sface->sh), sface->shver, shellmark(*sface));
5717  } else {
5718  printf("Subsegment x%lx, ver %d, mark %d:",
5719  (unsigned long)(sface->sh), sface->shver, shellmark(*sface));
5720  }
5721  if (sinfected(*sface)) {
5722  printf(" (infected)");
5723  }
5724  if (shell2badface(*sface)) {
5725  printf(" (queued)");
5726  }
5727  if (sapex(*sface) != NULL) {
5728  if (shelltype(*sface) == SHARP) {
5729  printf(" (sharp)");
5730  }
5731  } else {
5732  if (shelltype(*sface) == SHARP) {
5733  printf(" (sharp)");
5734  }
5735  }
5736  if (checkpbcs) {
5737  if (shellpbcgroup(*sface) >= 0) {
5738  printf(" (pbc %d)", shellpbcgroup(*sface));
5739  }
5740  }
5741  printf("\n");
5742 
5743  sdecode(sface->sh[0], prtsh);
5744  if (prtsh.sh == dummysh) {
5745  printf(" [0] = No shell\n");
5746  } else {
5747  printf(" [0] = x%lx %d\n", (unsigned long)(prtsh.sh), prtsh.shver);
5748  }
5749  sdecode(sface->sh[1], prtsh);
5750  if (prtsh.sh == dummysh) {
5751  printf(" [1] = No shell\n");
5752  } else {
5753  printf(" [1] = x%lx %d\n", (unsigned long)(prtsh.sh), prtsh.shver);
5754  }
5755  sdecode(sface->sh[2], prtsh);
5756  if (prtsh.sh == dummysh) {
5757  printf(" [2] = No shell\n");
5758  } else {
5759  printf(" [2] = x%lx %d\n", (unsigned long)(prtsh.sh), prtsh.shver);
5760  }
5761 
5762  printpoint = sorg(*sface);
5763  if (printpoint == (point) NULL)
5764  printf(" Org [%d] = NULL\n", vo[sface->shver]);
5765  else
5766  printf(" Org [%d] = x%lx (%.12g,%.12g,%.12g) %d\n",
5767  vo[sface->shver], (unsigned long)(printpoint), printpoint[0],
5768  printpoint[1], printpoint[2], pointmark(printpoint));
5769  printpoint = sdest(*sface);
5770  if (printpoint == (point) NULL)
5771  printf(" Dest[%d] = NULL\n", vd[sface->shver]);
5772  else
5773  printf(" Dest[%d] = x%lx (%.12g,%.12g,%.12g) %d\n",
5774  vd[sface->shver], (unsigned long)(printpoint), printpoint[0],
5775  printpoint[1], printpoint[2], pointmark(printpoint));
5776 
5777  if (sapex(*sface) != NULL) {
5778  printpoint = sapex(*sface);
5779  if (printpoint == (point) NULL)
5780  printf(" Apex[%d] = NULL\n", va[sface->shver]);
5781  else
5782  printf(" Apex[%d] = x%lx (%.12g,%.12g,%.12g) %d\n",
5783  va[sface->shver], (unsigned long)(printpoint), printpoint[0],
5784  printpoint[1], printpoint[2], pointmark(printpoint));
5785 
5786  decode(sface->sh[6], prttet);
5787  if (prttet.tet == dummytet) {
5788  printf(" [6] = Outer space\n");
5789  } else {
5790  printf(" [6] = x%lx %d\n",
5791  (unsigned long)(prttet.tet), prttet.loc);
5792  }
5793  decode(sface->sh[7], prttet);
5794  if (prttet.tet == dummytet) {
5795  printf(" [7] = Outer space\n");
5796  } else {
5797  printf(" [7] = x%lx %d\n",
5798  (unsigned long)(prttet.tet), prttet.loc);
5799  }
5800 
5801  sdecode(sface->sh[8], prtsh);
5802  if (prtsh.sh == dummysh) {
5803  printf(" [8] = No subsegment\n");
5804  } else {
5805  printf(" [8] = x%lx %d\n",
5806  (unsigned long)(prtsh.sh), prtsh.shver);
5807  }
5808  sdecode(sface->sh[9], prtsh);
5809  if (prtsh.sh == dummysh) {
5810  printf(" [9] = No subsegment\n");
5811  } else {
5812  printf(" [9] = x%lx %d\n",
5813  (unsigned long)(prtsh.sh), prtsh.shver);
5814  }
5815  sdecode(sface->sh[10], prtsh);
5816  if (prtsh.sh == dummysh) {
5817  printf(" [10]= No subsegment\n");
5818  } else {
5819  printf(" [10]= x%lx %d\n",
5820  (unsigned long)(prtsh.sh), prtsh.shver);
5821  }
5822  }
5823 }
5824 
5825 //
5826 // End of advanced primitives
5827 //
5828 
5829 //
5830 // End of mesh manipulation primitives
5831 //
5832 
5833 //
5834 // Begin of mesh items searching routines
5835 //
5836 
5838 // //
5839 // makepoint2tetmap() Construct a mapping from points to tetrahedra. //
5840 // //
5841 // Traverses all the tetrahedra, provides each corner of each tetrahedron //
5842 // with a pointer to that tetrahedera. Some pointers will be overwritten by //
5843 // other pointers because each point may be a corner of several tetrahedra, //
5844 // but in the end every point will point to a tetrahedron that contains it. //
5845 // //
5847 
5848 void tetgenmesh::makepoint2tetmap()
5849 {
5850  triface tetloop;
5851  point pointptr;
5852 
5853  if (b->verbose > 0) {
5854  printf(" Constructing mapping from points to tetrahedra.\n");
5855  }
5856 
5857  // Initialize the point2tet field of each point.
5858  points->traversalinit();
5859  pointptr = pointtraverse();
5860  while (pointptr != (point) NULL) {
5861  setpoint2tet(pointptr, (tetrahedron) NULL);
5862  pointptr = pointtraverse();
5863  }
5864 
5865  tetrahedrons->traversalinit();
5866  tetloop.tet = tetrahedrontraverse();
5867  while (tetloop.tet != (tetrahedron *) NULL) {
5868  // Check all four points of the tetrahedron.
5869  tetloop.loc = 0;
5870  pointptr = org(tetloop);
5871  setpoint2tet(pointptr, encode(tetloop));
5872  pointptr = dest(tetloop);
5873  setpoint2tet(pointptr, encode(tetloop));
5874  pointptr = apex(tetloop);
5875  setpoint2tet(pointptr, encode(tetloop));
5876  pointptr = oppo(tetloop);
5877  setpoint2tet(pointptr, encode(tetloop));
5878  // Get the next tetrahedron in the list.
5879  tetloop.tet = tetrahedrontraverse();
5880  }
5881 }
5882 
5884 // //
5885 // makeindex2pointmap() Create a map from index to vertices. //
5886 // //
5887 // 'idx2verlist' returns the created map. Traverse all vertices, a pointer //
5888 // to each vertex is set into the array. The pointer to the first vertex is //
5889 // saved in 'idx2verlist[0]'. Don't forget to minus 'in->firstnumber' when //
5890 // to get the vertex form its index. //
5891 // //
5893 
5894 void tetgenmesh::makeindex2pointmap(point*& idx2verlist)
5895 {
5896  point pointloop;
5897  int idx;
5898 
5899  if (b->verbose > 0) {
5900  printf(" Constructing mapping from indices to points.\n");
5901  }
5902 
5903  idx2verlist = new point[points->items];
5904 
5905  points->traversalinit();
5906  pointloop = pointtraverse();
5907  idx = 0;
5908  while (pointloop != (point) NULL) {
5909  idx2verlist[idx] = pointloop;
5910  idx++;
5911  pointloop = pointtraverse();
5912  }
5913 }
5914 
5916 // //
5917 // makesegmentmap() Create a map from vertices (their indices) to //
5918 // segments incident at the same vertices. //
5919 // //
5920 // Two arrays 'idx2seglist' and 'segsperverlist' together return the map. //
5921 // They form a sparse matrix structure with size (n + 1) x (n + 1), n is the //
5922 // number of segments. idx2seglist contains row information and //
5923 // segsperverlist contains all (non-zero) elements. The i-th entry of //
5924 // idx2seglist is the starting position of i-th row's (non-zero) elements in //
5925 // segsperverlist. The number of elements of i-th row is calculated by the //
5926 // (i+1)-th entry minus i-th entry of idx2seglist. //
5927 // //
5928 // NOTE: These two arrays will be created inside this routine, don't forget //
5929 // to free them after using. //
5930 // //
5932 
5933 void tetgenmesh::makesegmentmap(int*& idx2seglist, shellface**& segsperverlist)
5934 {
5935  shellface *shloop;
5936  int i, j, k;
5937 
5938  if (b->verbose > 0) {
5939  printf(" Constructing mapping from points to segments.\n");
5940  }
5941 
5942  // Create and initialize 'idx2seglist'.
5943  idx2seglist = new int[points->items + 1];
5944  for (i = 0; i < points->items + 1; i++) idx2seglist[i] = 0;
5945 
5946  // Loop the set of segments once, counter the number of segments sharing
5947  // each vertex.
5948  subsegs->traversalinit();
5949  shloop = shellfacetraverse(subsegs);
5950  while (shloop != (shellface *) NULL) {
5951  // Increment the number of sharing segments for each endpoint.
5952  for (i = 0; i < 2; i++) {
5953  j = pointmark((point) shloop[3 + i]) - in->firstnumber;
5954  idx2seglist[j]++;
5955  }
5956  shloop = shellfacetraverse(subsegs);
5957  }
5958 
5959  // Calculate the total length of array 'facesperverlist'.
5960  j = idx2seglist[0];
5961  idx2seglist[0] = 0; // Array starts from 0 element.
5962  for (i = 0; i < points->items; i++) {
5963  k = idx2seglist[i + 1];
5964  idx2seglist[i + 1] = idx2seglist[i] + j;
5965  j = k;
5966  }
5967  // The total length is in the last unit of idx2seglist.
5968  segsperverlist = new shellface*[idx2seglist[i]];
5969  // Loop the set of segments again, set the info. of segments per vertex.
5970  subsegs->traversalinit();
5971  shloop = shellfacetraverse(subsegs);
5972  while (shloop != (shellface *) NULL) {
5973  for (i = 0; i < 2; i++) {
5974  j = pointmark((point) shloop[3 + i]) - in->firstnumber;
5975  segsperverlist[idx2seglist[j]] = shloop;
5976  idx2seglist[j]++;
5977  }
5978  shloop = shellfacetraverse(subsegs);
5979  }
5980  // Contents in 'idx2seglist' are shifted, now shift them back.
5981  for (i = points->items - 1; i >= 0; i--) {
5982  idx2seglist[i + 1] = idx2seglist[i];
5983  }
5984  idx2seglist[0] = 0;
5985 }
5986 
5988 // //
5989 // makesubfacemap() Create a map from vertices (their indices) to //
5990 // subfaces incident at the same vertices. //
5991 // //
5992 // Two arrays 'idx2facelist' and 'facesperverlist' together return the map. //
5993 // They form a sparse matrix structure with size (n + 1) x (n + 1), n is the //
5994 // number of subfaces. idx2facelist contains row information and //
5995 // facesperverlist contains all (non-zero) elements. The i-th entry of //
5996 // idx2facelist is the starting position of i-th row's(non-zero) elements in //
5997 // facesperverlist. The number of elements of i-th row is calculated by the //
5998 // (i+1)-th entry minus i-th entry of idx2facelist. //
5999 // //
6000 // NOTE: These two arrays will be created inside this routine, don't forget //
6001 // to free them after using. //
6002 // //
6004 
6005 void tetgenmesh::
6006 makesubfacemap(int*& idx2facelist, shellface**& facesperverlist)
6007 {
6008  shellface *shloop;
6009  int i, j, k;
6010 
6011  if (b->verbose > 0) {
6012  printf(" Constructing mapping from points to subfaces.\n");
6013  }
6014 
6015  // Create and initialize 'idx2facelist'.
6016  idx2facelist = new int[points->items + 1];
6017  for (i = 0; i < points->items + 1; i++) idx2facelist[i] = 0;
6018 
6019  // Loop the set of subfaces once, counter the number of subfaces sharing
6020  // each vertex.
6021  subfaces->traversalinit();
6022  shloop = shellfacetraverse(subfaces);
6023  while (shloop != (shellface *) NULL) {
6024  // Increment the number of sharing segments for each endpoint.
6025  for (i = 0; i < 3; i++) {
6026  j = pointmark((point) shloop[3 + i]) - in->firstnumber;
6027  idx2facelist[j]++;
6028  }
6029  shloop = shellfacetraverse(subfaces);
6030  }
6031 
6032  // Calculate the total length of array 'facesperverlist'.
6033  j = idx2facelist[0];
6034  idx2facelist[0] = 0; // Array starts from 0 element.
6035  for (i = 0; i < points->items; i++) {
6036  k = idx2facelist[i + 1];
6037  idx2facelist[i + 1] = idx2facelist[i] + j;
6038  j = k;
6039  }
6040  // The total length is in the last unit of idx2facelist.
6041  facesperverlist = new shellface*[idx2facelist[i]];
6042  // Loop the set of segments again, set the info. of segments per vertex.
6043  subfaces->traversalinit();
6044  shloop = shellfacetraverse(subfaces);
6045  while (shloop != (shellface *) NULL) {
6046  for (i = 0; i < 3; i++) {
6047  j = pointmark((point) shloop[3 + i]) - in->firstnumber;
6048  facesperverlist[idx2facelist[j]] = shloop;
6049  idx2facelist[j]++;
6050  }
6051  shloop = shellfacetraverse(subfaces);
6052  }
6053  // Contents in 'idx2facelist' are shifted, now shift them back.
6054  for (i = points->items - 1; i >= 0; i--) {
6055  idx2facelist[i + 1] = idx2facelist[i];
6056  }
6057  idx2facelist[0] = 0;
6058 }
6059 
6061 // //
6062 // maketetrahedronmap() Create a map from vertices (their indices) to //
6063 // tetrahedra incident at the same vertices. //
6064 // //
6065 // Two arrays 'idx2tetlist' and 'tetsperverlist' together return the map. //
6066 // They form a sparse matrix structure with size (n + 1) x (n + 1), n is the //
6067 // number of tetrahedra. idx2tetlist contains row information and //
6068 // tetsperverlist contains all (non-zero) elements. The i-th entry of //
6069 // idx2tetlist is the starting position of i-th row's (non-zero) elements in //
6070 // tetsperverlist. The number of elements of i-th row is calculated by the //
6071 // (i+1)-th entry minus i-th entry of idx2tetlist. //
6072 // //
6073 // NOTE: These two arrays will be created inside this routine, don't forget //
6074 // to free them after using. //
6075 // //
6077 
6078 void tetgenmesh::
6079 maketetrahedronmap(int*& idx2tetlist, tetrahedron**& tetsperverlist)
6080 {
6081  tetrahedron *tetloop;
6082  int i, j, k;
6083 
6084  if (b->verbose > 0) {
6085  printf(" Constructing mapping from points to tetrahedra.\n");
6086  }
6087 
6088  // Create and initialize 'idx2tetlist'.
6089  idx2tetlist = new int[points->items + 1];
6090  for (i = 0; i < points->items + 1; i++) idx2tetlist[i] = 0;
6091 
6092  // Loop the set of tetrahedra once, counter the number of tetrahedra
6093  // sharing each vertex.
6094  tetrahedrons->traversalinit();
6095  tetloop = tetrahedrontraverse();
6096  while (tetloop != (tetrahedron *) NULL) {
6097  // Increment the number of sharing tetrahedra for each endpoint.
6098  for (i = 0; i < 4; i++) {
6099  j = pointmark((point) tetloop[4 + i]) - in->firstnumber;
6100  idx2tetlist[j]++;
6101  }
6102  tetloop = tetrahedrontraverse();
6103  }
6104 
6105  // Calculate the total length of array 'tetsperverlist'.
6106  j = idx2tetlist[0];
6107  idx2tetlist[0] = 0; // Array starts from 0 element.
6108  for (i = 0; i < points->items; i++) {
6109  k = idx2tetlist[i + 1];
6110  idx2tetlist[i + 1] = idx2tetlist[i] + j;
6111  j = k;
6112  }
6113  // The total length is in the last unit of idx2tetlist.
6114  tetsperverlist = new tetrahedron*[idx2tetlist[i]];
6115  // Loop the set of tetrahedra again, set the info. of tet. per vertex.
6116  tetrahedrons->traversalinit();
6117  tetloop = tetrahedrontraverse();
6118  while (tetloop != (tetrahedron *) NULL) {
6119  for (i = 0; i < 4; i++) {
6120  j = pointmark((point) tetloop[4 + i]) - in->firstnumber;
6121  tetsperverlist[idx2tetlist[j]] = tetloop;
6122  idx2tetlist[j]++;
6123  }
6124  tetloop = tetrahedrontraverse();
6125  }
6126  // Contents in 'idx2tetlist' are shifted, now shift them back.
6127  for (i = points->items - 1; i >= 0; i--) {
6128  idx2tetlist[i + 1] = idx2tetlist[i];
6129  }
6130  idx2tetlist[0] = 0;
6131 }
6132 
6133 //
6134 // End of mesh items searching routines
6135 //
6136 
6137 //
6138 // Begin of linear algebra functions
6139 //
6140 
6141 // dot() returns the dot product: v1 dot v2.
6142 
6143 inline REAL tetgenmesh::dot(REAL* v1, REAL* v2)
6144 {
6145  return v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2];
6146 }
6147 
6148 // cross() computes the cross product: n = v1 cross v2.
6149 
6150 inline void tetgenmesh::cross(REAL* v1, REAL* v2, REAL* n)
6151 {
6152  n[0] = v1[1] * v2[2] - v2[1] * v1[2];
6153  n[1] = -(v1[0] * v2[2] - v2[0] * v1[2]);
6154  n[2] = v1[0] * v2[1] - v2[0] * v1[1];
6155 }
6156 
6157 // initm44() initializes a 4x4 matrix.
6158 static void initm44(REAL a00, REAL a01, REAL a02, REAL a03,
6159  REAL a10, REAL a11, REAL a12, REAL a13,
6160  REAL a20, REAL a21, REAL a22, REAL a23,
6161  REAL a30, REAL a31, REAL a32, REAL a33,
6162  REAL M[4][4])
6163 {
6164  M[0][0] = a00; M[0][1] = a01; M[0][2] = a02; M[0][3] = a03;
6165  M[1][0] = a10; M[1][1] = a11; M[1][2] = a12; M[1][3] = a13;
6166  M[2][0] = a20; M[2][1] = a21; M[2][2] = a22; M[2][3] = a23;
6167  M[3][0] = a30; M[3][1] = a31; M[3][2] = a32; M[3][3] = a33;
6168 }
6169 
6170 // m4xm4() multiplies 2 4x4 matrics: m1 = m1 * m2.
6171 static void m4xm4(REAL m1[4][4], REAL m2[4][4])
6172 {
6173  REAL tmp[4];
6174  int i, j;
6175 
6176  for (i = 0; i < 4; i++) { // i-th row
6177  for (j = 0; j < 4; j++) { // j-th col
6178  tmp[j] = m1[i][0] * m2[0][j] + m1[i][1] * m2[1][j]
6179  + m1[i][2] * m2[2][j] + m1[i][3] * m2[3][j];
6180  }
6181  for (j = 0; j < 4; j++)
6182  m1[i][j] = tmp[j];
6183  }
6184 }
6185 
6186 // m4xv4() multiplies a 4x4 matrix and 4x1 vector: v2 = m * v1
6187 static void m4xv4(REAL v2[4], REAL m[4][4], REAL v1[4])
6188 {
6189  v2[0] = m[0][0]*v1[0] + m[0][1]*v1[1] + m[0][2]*v1[2] + m[0][3]*v1[3];
6190  v2[1] = m[1][0]*v1[0] + m[1][1]*v1[1] + m[1][2]*v1[2] + m[1][3]*v1[3];
6191  v2[2] = m[2][0]*v1[0] + m[2][1]*v1[1] + m[2][2]*v1[2] + m[2][3]*v1[3];
6192  v2[3] = m[3][0]*v1[0] + m[3][1]*v1[1] + m[3][2]*v1[2] + m[3][3]*v1[3];
6193 }
6194 
6196 // //
6197 // lu_decmp() Compute the LU decomposition of a matrix. //
6198 // //
6199 // Compute the LU decomposition of a (non-singular) square matrix A using //
6200 // partial pivoting and implicit row exchanges. The result is: //
6201 // A = P * L * U, //
6202 // where P is a permutation matrix, L is unit lower triangular, and U is //
6203 // upper triangular. The factored form of A is used in combination with //
6204 // 'lu_solve()' to solve linear equations: Ax = b, or invert a matrix. //
6205 // //
6206 // The inputs are a square matrix 'lu[N..n+N-1][N..n+N-1]', it's size is 'n'.//
6207 // On output, 'lu' is replaced by the LU decomposition of a rowwise permuta- //
6208 // tion of itself, 'ps[N..n+N-1]' is an output vector that records the row //
6209 // permutation effected by the partial pivoting, effectively, 'ps' array //
6210 // tells the user what the permutation matrix P is; 'd' is output as +1/-1 //
6211 // depending on whether the number of row interchanges was even or odd, //
6212 // respectively. //
6213 // //
6214 // Return true if the LU decomposition is successfully computed, otherwise, //
6215 // return false in case that A is a singular matrix. //
6216 // //
6218 
6219 bool tetgenmesh::lu_decmp(REAL lu[4][4], int n, int* ps, REAL* d, int N)
6220 {
6221  REAL scales[4];
6222  REAL pivot, biggest, mult, tempf;
6223  int pivotindex = 0;
6224  int i, j, k;
6225 
6226  *d = 1.0; // No row interchanges yet.
6227 
6228  for (i = N; i < n + N; i++) { // For each row.
6229  // Find the largest element in each row for row equilibration
6230  biggest = 0.0;
6231  for (j = N; j < n + N; j++)
6232  if (biggest < (tempf = fabs(lu[i][j])))
6233  biggest = tempf;
6234  if (biggest != 0.0)
6235  scales[i] = 1.0 / biggest;
6236  else {
6237  scales[i] = 0.0;
6238  return false; // Zero row: singular matrix.
6239  }
6240  ps[i] = i; // Initialize pivot sequence.
6241  }
6242 
6243  for (k = N; k < n + N - 1; k++) { // For each column.
6244  // Find the largest element in each column to pivot around.
6245  biggest = 0.0;
6246  for (i = k; i < n + N; i++) {
6247  if (biggest < (tempf = fabs(lu[ps[i]][k]) * scales[ps[i]])) {
6248  biggest = tempf;
6249  pivotindex = i;
6250  }
6251  }
6252  if (biggest == 0.0) {
6253  return false; // Zero column: singular matrix.
6254  }
6255  if (pivotindex != k) { // Update pivot sequence.
6256  j = ps[k];
6257  ps[k] = ps[pivotindex];
6258  ps[pivotindex] = j;
6259  *d = -(*d); // ...and change the parity of d.
6260  }
6261 
6262  // Pivot, eliminating an extra variable each time
6263  pivot = lu[ps[k]][k];
6264  for (i = k + 1; i < n + N; i++) {
6265  lu[ps[i]][k] = mult = lu[ps[i]][k] / pivot;
6266  if (mult != 0.0) {
6267  for (j = k + 1; j < n + N; j++)
6268  lu[ps[i]][j] -= mult * lu[ps[k]][j];
6269  }
6270  }
6271  }
6272 
6273  // (lu[ps[n + N - 1]][n + N - 1] == 0.0) ==> A is singular.
6274  return lu[ps[n + N - 1]][n + N - 1] != 0.0;
6275 }
6276 
6278 // //
6279 // lu_solve() Solves the linear equation: Ax = b, after the matrix A //
6280 // has been decomposed into the lower and upper triangular //
6281 // matrices L and U, where A = LU. //
6282 // //
6283 // 'lu[N..n+N-1][N..n+N-1]' is input, not as the matrix 'A' but rather as //
6284 // its LU decomposition, computed by the routine 'lu_decmp'; 'ps[N..n+N-1]' //
6285 // is input as the permutation vector returned by 'lu_decmp'; 'b[N..n+N-1]' //
6286 // is input as the right-hand side vector, and returns with the solution //
6287 // vector. 'lu', 'n', and 'ps' are not modified by this routine and can be //
6288 // left in place for successive calls with different right-hand sides 'b'. //
6289 // //
6291 
6292 void tetgenmesh::lu_solve(REAL lu[4][4], int n, int* ps, REAL* b, int N)
6293 {
6294  int i, j;
6295  REAL X[4], dot;
6296 
6297  for (i = N; i < n + N; i++) X[i] = 0.0;
6298 
6299  // Vector reduction using U triangular matrix.
6300  for (i = N; i < n + N; i++) {
6301  dot = 0.0;
6302  for (j = N; j < i + N; j++)
6303  dot += lu[ps[i]][j] * X[j];
6304  X[i] = b[ps[i]] - dot;
6305  }
6306 
6307  // Back substitution, in L triangular matrix.
6308  for (i = n + N - 1; i >= N; i--) {
6309  dot = 0.0;
6310  for (j = i + 1; j < n + N; j++)
6311  dot += lu[ps[i]][j] * X[j];
6312  X[i] = (X[i] - dot) / lu[ps[i]][i];
6313  }
6314 
6315  for (i = N; i < n + N; i++) b[i] = X[i];
6316 }
6317 
6318 //
6319 // End of linear algebra functions
6320 //
6321 
6322 //
6323 // Begin of geometric tests
6324 //
6325 
6326 // All the following routines require the input objects are not degenerate.
6327 // i.e., a triangle must has three non-collinear corners; an edge must
6328 // has two identical endpoints. Degenerate cases should have to detect
6329 // first and then handled as special cases.
6330 
6332 // //
6333 // edge_vert_col_inter() Test whether an edge (ab) and a collinear vertex //
6334 // (p) are intersecting or not. //
6335 // //
6336 // Possible cases are p is coincident to a (p = a), or to b (p = b), or p is //
6337 // inside ab (a < p < b), or outside ab (p < a or p > b). These cases can be //
6338 // quickly determined by comparing the corresponding coords of a, b, and p //
6339 // (which are not all equal). //
6340 // //
6341 // The return value indicates one of the three cases: DISJOINT, SHAREVERTEX //
6342 // (p = a or p = b), and INTERSECT (a < p < b). //
6343 // //
6345 
6346 enum tetgenmesh::interresult tetgenmesh::edge_vert_col_inter(REAL* A, REAL* B,
6347  REAL* P)
6348 {
6349  int i = 0;
6350  do {
6351  if (A[i] < B[i]) {
6352  if (P[i] < A[i]) {
6353  return DISJOINT;
6354  } else if (P[i] > A[i]) {
6355  if (P[i] < B[i]) {
6356  return INTERSECT;
6357  } else if (P[i] > B[i]) {
6358  return DISJOINT;
6359  } else {
6360  // assert(P[i] == B[i]);
6361  return SHAREVERTEX;
6362  }
6363  } else {
6364  // assert(P[i] == A[i]);
6365  return SHAREVERTEX;
6366  }
6367  } else if (A[i] > B[i]) {
6368  if (P[i] < B[i]) {
6369  return DISJOINT;
6370  } else if (P[i] > B[i]) {
6371  if (P[i] < A[i]) {
6372  return INTERSECT;
6373  } else if (P[i] > A[i]) {
6374  return DISJOINT;
6375  } else {
6376  // assert(P[i] == A[i]);
6377  return SHAREVERTEX;
6378  }
6379  } else {
6380  // assert(P[i] == B[i]);
6381  return SHAREVERTEX;
6382  }
6383  }
6384  // i-th coordinates are equal, try i+1-th;
6385  i++;
6386  } while (i < 3);
6387  // Should never be here.
6388  return DISJOINT;
6389 }
6390 
6392 // //
6393 // edge_edge_cop_inter() Test whether two coplanar edges (ab, and pq) are //
6394 // intersecting or not. //
6395 // //
6396 // Possible cases are ab and pq are disjointed, or proper intersecting (int- //
6397 // ersect at a point other than their endpoints), or both collinear and int- //
6398 // ersecting, or sharing at a common endpoint, or are coincident. //
6399 // //
6400 // A reference point R is required, which is exactly not coplanar with these //
6401 // two edges. Since the caller knows these two edges are coplanar, it must //
6402 // be able to provide (or calculate) such a point. //
6403 // //
6404 // The return value indicates one of the four cases: DISJOINT, SHAREVERTEX, //
6405 // SHAREEDGE, and INTERSECT. //
6406 // //
6408 
6409 enum tetgenmesh::interresult tetgenmesh:: edge_edge_cop_inter(REAL* A, REAL* B,
6410  REAL* P, REAL* Q, REAL* R)
6411 {
6412  REAL s1, s2, s3, s4;
6413 
6414 #ifdef SELF_CHECK
6415  assert(R != NULL);
6416 #endif
6417  s1 = orient3d(A, B, R, P);
6418  s2 = orient3d(A, B, R, Q);
6419  if (s1 * s2 > 0.0) {
6420  // Both p and q are at the same side of ab.
6421  return DISJOINT;
6422  }
6423  s3 = orient3d(P, Q, R, A);
6424  s4 = orient3d(P, Q, R, B);
6425  if (s3 * s4 > 0.0) {
6426  // Both a and b are at the same side of pq.
6427  return DISJOINT;
6428  }
6429 
6430  // Possible degenerate cases are:
6431  // (1) Only one of p and q is collinear with ab;
6432  // (2) Both p and q are collinear with ab;
6433  // (3) Only one of a and b is collinear with pq.
6434  enum interresult abp, abq;
6435  enum interresult pqa, pqb;
6436 
6437  if (s1 == 0.0) {
6438  // p is collinear with ab.
6439  abp = edge_vert_col_inter(A, B, P);
6440  if (abp == INTERSECT) {
6441  // p is inside ab.
6442  return INTERSECT;
6443  }
6444  if (s2 == 0.0) {
6445  // q is collinear with ab. Case (2).
6446  abq = edge_vert_col_inter(A, B, Q);
6447  if (abq == INTERSECT) {
6448  // q is inside ab.
6449  return INTERSECT;
6450  }
6451  if (abp == SHAREVERTEX && abq == SHAREVERTEX) {
6452  // ab and pq are identical.
6453  return SHAREEDGE;
6454  }
6455  pqa = edge_vert_col_inter(P, Q, A);
6456  if (pqa == INTERSECT) {
6457  // a is inside pq.
6458  return INTERSECT;
6459  }
6460  pqb = edge_vert_col_inter(P, Q, B);
6461  if (pqb == INTERSECT) {
6462  // b is inside pq.
6463  return INTERSECT;
6464  }
6465  if (abp == SHAREVERTEX || abq == SHAREVERTEX) {
6466  // either p or q is coincident with a or b.
6467 #ifdef SELF_CHECK
6468  // ONLY one case is possible, otherwise, shoule be SHAREEDGE.
6469  assert(abp ^ abq);
6470 #endif
6471  return SHAREVERTEX;
6472  }
6473  // The last case. They are disjointed.
6474 #ifdef SELF_CHECK
6475  assert((abp == DISJOINT) && (abp == abq && abq == pqa && pqa == pqb));
6476 #endif
6477  return DISJOINT;
6478  } else {
6479  // p is collinear with ab. Case (1).
6480 #ifdef SELF_CHECK
6481  assert(abp == SHAREVERTEX || abp == DISJOINT);
6482 #endif
6483  return abp;
6484  }
6485  }
6486  // p is NOT collinear with ab.
6487  if (s2 == 0.0) {
6488  // q is collinear with ab. Case (1).
6489  abq = edge_vert_col_inter(A, B, Q);
6490 #ifdef SELF_CHECK
6491  assert(abq == SHAREVERTEX || abq == DISJOINT || abq == INTERSECT);
6492 #endif
6493  return abq;
6494  }
6495 
6496  // We have found p and q are not collinear with ab. However, it is still
6497  // possible that a or b is collinear with pq (ONLY one of a and b).
6498  if (s3 == 0.0) {
6499  // a is collinear with pq. Case (3).
6500 #ifdef SELF_CHECK
6501  assert(s4 != 0.0);
6502 #endif
6503  pqa = edge_vert_col_inter(P, Q, A);
6504 #ifdef SELF_CHECK
6505  // This case should have been detected in above.
6506  assert(pqa != SHAREVERTEX);
6507  assert(pqa == INTERSECT || pqa == DISJOINT);
6508 #endif
6509  return pqa;
6510  }
6511  if (s4 == 0.0) {
6512  // b is collinear with pq. Case (3).
6513 #ifdef SELF_CHECK
6514  assert(s3 != 0.0);
6515 #endif
6516  pqb = edge_vert_col_inter(P, Q, B);
6517 #ifdef SELF_CHECK
6518  // This case should have been detected in above.
6519  assert(pqb != SHAREVERTEX);
6520  assert(pqb == INTERSECT || pqb == DISJOINT);
6521 #endif
6522  return pqb;
6523  }
6524 
6525  // ab and pq are intersecting properly.
6526  return INTERSECT;
6527 }
6528 
6530 // //
6531 // Notations //
6532 // //
6533 // Let ABC be the plane passes through a, b, and c; ABC+ be the halfspace //
6534 // including the set of all points x, such that orient3d(a, b, c, x) > 0; //
6535 // ABC- be the other halfspace, such that for each point x in ABC-, //
6536 // orient3d(a, b, c, x) < 0. For the set of x which are on ABC, orient3d(a, //
6537 // b, c, x) = 0. //
6538 // //
6540 
6542 // //
6543 // tri_vert_copl_inter() Test whether a triangle (abc) and a coplanar //
6544 // point (p) are intersecting or not. //
6545 // //
6546 // Possible cases are p is inside abc, or on an edge of, or coincident with //
6547 // a vertex of, or outside abc. //
6548 // //
6549 // A reference point R is required. R is exactly not coplanar with abc and p.//
6550 // Since the caller knows they are coplanar, it must be able to provide (or //
6551 // calculate) such a point. //
6552 // //
6553 // The return value indicates one of the four cases: DISJOINT, SHAREVERTEX, //
6554 // and INTERSECT. //
6555 // //
6557 
6558 enum tetgenmesh::interresult tetgenmesh::tri_vert_cop_inter(REAL* A, REAL* B,
6559  REAL* C, REAL* P, REAL* R)
6560 {
6561  REAL s1, s2, s3;
6562  int sign;
6563 
6564 #ifdef SELF_CHECK
6565  assert(R != (REAL *) NULL);
6566 #endif
6567  // Adjust the orientation of a, b, c and r, so that we can assume that
6568  // r is strictly in ABC- (i.e., r is above ABC wrt. right-hand rule).
6569  s1 = orient3d(A, B, C, R);
6570 #ifdef SELF_CHECK
6571  assert(s1 != 0.0);
6572 #endif
6573  sign = s1 < 0.0 ? 1 : -1;
6574 
6575  // Test starts from here.
6576  s1 = orient3d(A, B, R, P) * sign;
6577  if (s1 < 0.0) {
6578  // p is in ABR-.
6579  return DISJOINT;
6580  }
6581  s2 = orient3d(B, C, R, P) * sign;
6582  if (s2 < 0.0) {
6583  // p is in BCR-.
6584  return DISJOINT;
6585  }
6586  s3 = orient3d(C, A, R, P) * sign;
6587  if (s3 < 0.0) {
6588  // p is in CAR-.
6589  return DISJOINT;
6590  }
6591  if (s1 == 0.0) {
6592  // p is on ABR.
6593  if (s2 == 0.0) {
6594  // p is on BCR.
6595 #ifdef SELF_CHECK
6596  assert(s3 > 0.0);
6597 #endif
6598  // p is coincident with b.
6599  return SHAREVERTEX;
6600  }
6601  if (s3 == 0.0) {
6602  // p is on CAR.
6603  // p is coincident with a.
6604  return SHAREVERTEX;
6605  }
6606  // p is on edge ab.
6607  return INTERSECT;
6608  }
6609  // p is in ABR+.
6610  if (s2 == 0.0) {
6611  // p is on BCR.
6612  if (s3 == 0.0) {
6613  // p is on CAR.
6614  // p is coincident with c.
6615  return SHAREVERTEX;
6616  }
6617  // p is on edge bc.
6618  return INTERSECT;
6619  }
6620  if (s3 == 0.0) {
6621  // p is on CAR.
6622  // p is on edge ca.
6623  return INTERSECT;
6624  }
6625 
6626  // p is strictly inside abc.
6627  return INTERSECT;
6628 }
6629 
6631 // //
6632 // tri_edge_cop_inter() Test whether a triangle (abc) and a coplanar edge //
6633 // (pq) are intersecting or not. //
6634 // //
6635 // A reference point R is required. R is exactly not coplanar with abc and //
6636 // pq. Since the caller knows they are coplanar, it must be able to provide //
6637 // (or calculate) such a point. //
6638 // //
6639 // The return value indicates one of the four cases: DISJOINT, SHAREVERTEX, //
6640 // SHAREEDGE, and INTERSECT. //
6641 // //
6643 
6644 enum tetgenmesh::interresult tetgenmesh::tri_edge_cop_inter(REAL* A, REAL* B,
6645  REAL* C, REAL* P, REAL* Q, REAL* R)
6646 {
6647  enum interresult abpq, bcpq, capq;
6648  enum interresult abcp, abcq;
6649 
6650  // Test if pq is intersecting one of edges of abc.
6651  abpq = edge_edge_cop_inter(A, B, P, Q, R);
6652  if (abpq == INTERSECT || abpq == SHAREEDGE) {
6653  return abpq;
6654  }
6655  bcpq = edge_edge_cop_inter(B, C, P, Q, R);
6656  if (bcpq == INTERSECT || bcpq == SHAREEDGE) {
6657  return bcpq;
6658  }
6659  capq = edge_edge_cop_inter(C, A, P, Q, R);
6660  if (capq == INTERSECT || capq == SHAREEDGE) {
6661  return capq;
6662  }
6663 
6664  // Test if p and q is inside abc.
6665  abcp = tri_vert_cop_inter(A, B, C, P, R);
6666  if (abcp == INTERSECT) {
6667  return INTERSECT;
6668  }
6669  abcq = tri_vert_cop_inter(A, B, C, Q, R);
6670  if (abcq == INTERSECT) {
6671  return INTERSECT;
6672  }
6673 
6674  // Combine the test results of edge intersectings and triangle insides
6675  // to detect whether abc and pq are sharing vertex or disjointed.
6676  if (abpq == SHAREVERTEX) {
6677  // p or q is coincident with a or b.
6678 #ifdef SELF_CHECK
6679  assert(abcp ^ abcq);
6680 #endif
6681  return SHAREVERTEX;
6682  }
6683  if (bcpq == SHAREVERTEX) {
6684  // p or q is coincident with b or c.
6685 #ifdef SELF_CHECK
6686  assert(abcp ^ abcq);
6687 #endif
6688  return SHAREVERTEX;
6689  }
6690  if (capq == SHAREVERTEX) {
6691  // p or q is coincident with c or a.
6692 #ifdef SELF_CHECK
6693  assert(abcp ^ abcq);
6694 #endif
6695  return SHAREVERTEX;
6696  }
6697 
6698  // They are disjointed.
6699  return DISJOINT;
6700 }
6701 
6703 // //
6704 // tri_edge_inter_tail() Test whether a triangle (abc) and an edge (pq) //
6705 // are intersecting or not. //
6706 // //
6707 // s1 and s2 are results of pre-performed orientation tests. s1 = orient3d( //
6708 // a, b, c, p); s2 = orient3d(a, b, c, q). To separate this routine from //
6709 // tri_edge_inter() can save two orientation tests in tri_tri_inter(). //
6710 // //
6711 // The return value indicates one of the four cases: DISJOINT, SHAREVERTEX, //
6712 // SHAREEDGE, and INTERSECT. //
6713 // //
6715 
6716 enum tetgenmesh::interresult tetgenmesh::tri_edge_inter_tail(REAL* A, REAL* B,
6717  REAL* C, REAL* P, REAL* Q, REAL s1, REAL s2)
6718 {
6719  REAL s3, s4, s5;
6720  int sign;
6721 
6722  if (s1 * s2 > 0.0) {
6723  // p, q are at the same halfspace of ABC, no intersection.
6724  return DISJOINT;
6725  }
6726 
6727  if (s1 * s2 < 0.0) {
6728  // p, q are both not on ABC (and not sharing vertices, edges of abc).
6729  // Adjust the orientation of a, b, c and p, so that we can assume that
6730  // p is strictly in ABC-, and q is strictly in ABC+.
6731  sign = s1 < 0.0 ? 1 : -1;
6732  s3 = orient3d(A, B, P, Q) * sign;
6733  if (s3 < 0.0) {
6734  // q is at ABP-.
6735  return DISJOINT;
6736  }
6737  s4 = orient3d(B, C, P, Q) * sign;
6738  if (s4 < 0.0) {
6739  // q is at BCP-.
6740  return DISJOINT;
6741  }
6742  s5 = orient3d(C, A, P, Q) * sign;
6743  if (s5 < 0.0) {
6744  // q is at CAP-.
6745  return DISJOINT;
6746  }
6747  if (s3 == 0.0) {
6748  // q is on ABP.
6749  if (s4 == 0.0) {
6750  // q is on BCP (and q must in CAP+).
6751 #ifdef SELF_CHECK
6752  assert(s5 > 0.0);
6753 #endif
6754  // pq intersects abc at vertex b.
6755  return SHAREVERTEX;
6756  }
6757  if (s5 == 0.0) {
6758  // q is on CAP (and q must in BCP+).
6759  // pq intersects abc at vertex a.
6760  return SHAREVERTEX;
6761  }
6762  // q in both BCP+ and CAP+.
6763  // pq crosses ab properly.
6764  return INTERSECT;
6765  }
6766  // q is in ABP+;
6767  if (s4 == 0.0) {
6768  // q is on BCP.
6769  if (s5 == 0.0) {
6770  // q is on CAP.
6771  // pq intersects abc at vertex c.
6772  return SHAREVERTEX;
6773  }
6774  // pq crosses bc properly.
6775  return INTERSECT;
6776  }
6777  // q is in BCP+;
6778  if (s5 == 0.0) {
6779  // q is on CAP.
6780  // pq crosses ca properly.
6781  return INTERSECT;
6782  }
6783  // q is in CAP+;
6784  // pq crosses abc properly.
6785  return INTERSECT;
6786  }
6787 
6788  if (s1 != 0.0 || s2 != 0.0) {
6789  // Either p or q is coplanar with abc. ONLY one of them is possible.
6790  if (s1 == 0.0) {
6791  // p is coplanar with abc, q can be used as reference point.
6792 #ifdef SELF_CHECK
6793  assert(s2 != 0.0);
6794 #endif
6795  return tri_vert_cop_inter(A, B, C, P, Q);
6796  } else {
6797  // q is coplanar with abc, p can be used as reference point.
6798 #ifdef SELF_CHECK
6799  assert(s2 == 0.0);
6800 #endif
6801  return tri_vert_cop_inter(A, B, C, Q, P);
6802  }
6803  }
6804 
6805  // pq is coplanar with abc. Calculate a point which is exactly not
6806  // coplanar with a, b, and c.
6807  REAL R[3], N[3];
6808  REAL ax, ay, az, bx, by, bz;
6809 
6810  ax = A[0] - B[0];
6811  ay = A[1] - B[1];
6812  az = A[2] - B[2];
6813  bx = A[0] - C[0];
6814  by = A[1] - C[1];
6815  bz = A[2] - C[2];
6816  N[0] = ay * bz - by * az;
6817  N[1] = az * bx - bz * ax;
6818  N[2] = ax * by - bx * ay;
6819  // The normal should not be a zero vector (otherwise, abc are collinear).
6820 #ifdef SELF_CHECK
6821  assert((fabs(N[0]) + fabs(N[1]) + fabs(N[2])) > 0.0);
6822 #endif
6823  // The reference point R is lifted from A to the normal direction with
6824  // a distance d = average edge length of the triangle abc.
6825  R[0] = N[0] + A[0];
6826  R[1] = N[1] + A[1];
6827  R[2] = N[2] + A[2];
6828  // Becareful the case: if the non-zero component(s) in N is smaller than
6829  // the machine epsilon (i.e., 2^(-16) for double), R will exactly equal
6830  // to A due to the round-off error. Do check if it is.
6831  if (R[0] == A[0] && R[1] == A[1] && R[2] == A[2]) {
6832  int i, j;
6833  for (i = 0; i < 3; i++) {
6834 #ifdef SELF_CHECK
6835  assert (R[i] == A[i]);
6836 #endif
6837  j = 2;
6838  do {
6839  if (N[i] > 0.0) {
6840  N[i] += (j * macheps);
6841  } else {
6842  N[i] -= (j * macheps);
6843  }
6844  R[i] = N[i] + A[i];
6845  j *= 2;
6846  } while (R[i] == A[i]);
6847  }
6848  }
6849 
6850  return tri_edge_cop_inter(A, B, C, P, Q, R);
6851 }
6852 
6854 // //
6855 // tri_edge_inter() Test whether a triangle (abc) and an edge (pq) are //
6856 // intersecting or not. //
6857 // //
6858 // The return value indicates one of the four cases: DISJOINT, SHAREVERTEX, //
6859 // SHAREEDGE, and INTERSECT. //
6860 // //
6862 
6863 enum tetgenmesh::interresult tetgenmesh::tri_edge_inter(REAL* A, REAL* B,
6864  REAL* C, REAL* P, REAL* Q)
6865 {
6866  REAL s1, s2;
6867 
6868  // Test the locations of p and q with respect to ABC.
6869  s1 = orient3d(A, B, C, P);
6870  s2 = orient3d(A, B, C, Q);
6871 
6872  return tri_edge_inter_tail(A, B, C, P, Q, s1, s2);
6873 }
6874 
6876 // //
6877 // tri_tri_inter() Test whether two triangle (abc) and (opq) are //
6878 // intersecting or not. //
6879 // //
6880 // The return value indicates one of the five cases: DISJOINT, SHAREVERTEX, //
6881 // SHAREEDGE, SHAREFACE, and INTERSECT. //
6882 // //
6884 
6885 enum tetgenmesh::interresult tetgenmesh::tri_tri_inter(REAL* A, REAL* B,
6886  REAL* C, REAL* O, REAL* P, REAL* Q)
6887 {
6888  REAL s_o, s_p, s_q;
6889  REAL s_a, s_b, s_c;
6890 
6891  s_o = orient3d(A, B, C, O);
6892  s_p = orient3d(A, B, C, P);
6893  s_q = orient3d(A, B, C, Q);
6894  if ((s_o * s_p > 0.0) && (s_o * s_q > 0.0)) {
6895  // o, p, q are all in the same halfspace of ABC.
6896  return DISJOINT;
6897  }
6898 
6899  s_a = orient3d(O, P, Q, A);
6900  s_b = orient3d(O, P, Q, B);
6901  s_c = orient3d(O, P, Q, C);
6902  if ((s_a * s_b > 0.0) && (s_a * s_c > 0.0)) {
6903  // a, b, c are all in the same halfspace of OPQ.
6904  return DISJOINT;
6905  }
6906 
6907  enum interresult abcop, abcpq, abcqo;
6908  int shareedge = 0;
6909 
6910  abcop = tri_edge_inter_tail(A, B, C, O, P, s_o, s_p);
6911  if (abcop == INTERSECT) {
6912  return INTERSECT;
6913  } else if (abcop == SHAREEDGE) {
6914  shareedge++;
6915  }
6916  abcpq = tri_edge_inter_tail(A, B, C, P, Q, s_p, s_q);
6917  if (abcpq == INTERSECT) {
6918  return INTERSECT;
6919  } else if (abcpq == SHAREEDGE) {
6920  shareedge++;
6921  }
6922  abcqo = tri_edge_inter_tail(A, B, C, Q, O, s_q, s_o);
6923  if (abcqo == INTERSECT) {
6924  return INTERSECT;
6925  } else if (abcqo == SHAREEDGE) {
6926  shareedge++;
6927  }
6928  if (shareedge == 3) {
6929  // opq are coincident with abc.
6930  return SHAREFACE;
6931  }
6932 #ifdef SELF_CHECK
6933  // It is only possible either no share edge or one.
6934  assert(shareedge == 0 || shareedge == 1);
6935 #endif
6936 
6937  // Continue to detect whether opq and abc are intersecting or not.
6938  enum interresult opqab, opqbc, opqca;
6939 
6940  opqab = tri_edge_inter_tail(O, P, Q, A, B, s_a, s_b);
6941  if (opqab == INTERSECT) {
6942  return INTERSECT;
6943  }
6944  opqbc = tri_edge_inter_tail(O, P, Q, B, C, s_b, s_c);
6945  if (opqbc == INTERSECT) {
6946  return INTERSECT;
6947  }
6948  opqca = tri_edge_inter_tail(O, P, Q, C, A, s_c, s_a);
6949  if (opqca == INTERSECT) {
6950  return INTERSECT;
6951  }
6952 
6953  // At this point, two triangles are not intersecting and not coincident.
6954  // They may be share an edge, or share a vertex, or disjoint.
6955  if (abcop == SHAREEDGE) {
6956 #ifdef SELF_CHECK
6957  assert(abcpq == SHAREVERTEX && abcqo == SHAREVERTEX);
6958 #endif
6959  // op is coincident with an edge of abc.
6960  return SHAREEDGE;
6961  }
6962  if (abcpq == SHAREEDGE) {
6963 #ifdef SELF_CHECK
6964  assert(abcop == SHAREVERTEX && abcqo == SHAREVERTEX);
6965 #endif
6966  // pq is coincident with an edge of abc.
6967  return SHAREEDGE;
6968  }
6969  if (abcqo == SHAREEDGE) {
6970 #ifdef SELF_CHECK
6971  assert(abcop == SHAREVERTEX && abcpq == SHAREVERTEX);
6972 #endif
6973  // qo is coincident with an edge of abc.
6974  return SHAREEDGE;
6975  }
6976 
6977  // They may share a vertex or disjoint.
6978  if (abcop == SHAREVERTEX) {
6979  // o or p is coincident with a vertex of abc.
6980  if (abcpq == SHAREVERTEX) {
6981  // p is the coincident vertex.
6982 #ifdef SELF_CHECK
6983  assert(abcqo != SHAREVERTEX);
6984 #endif
6985  } else {
6986  // o is the coincident vertex.
6987 #ifdef SELF_CHECK
6988  assert(abcqo == SHAREVERTEX);
6989 #endif
6990  }
6991  return SHAREVERTEX;
6992  }
6993  if (abcpq == SHAREVERTEX) {
6994  // q is the coincident vertex.
6995 #ifdef SELF_CHECK
6996  assert(abcqo == SHAREVERTEX);
6997 #endif
6998  return SHAREVERTEX;
6999  }
7000 
7001  // They are disjoint.
7002  return DISJOINT;
7003 }
7004 
7006 // //
7007 // insphere_sos() Insphere test with symbolic perturbation. //
7008 // //
7009 // The input points a, b, c, and d should be non-coplanar. They must be ord- //
7010 // ered so that they have a positive orientation (as defined by orient3d()), //
7011 // or the sign of the result will be reversed. //
7012 // //
7013 // Return a positive value if the point e lies inside the circumsphere of a, //
7014 // b, c, and d; a negative value if it lies outside. //
7015 // //
7017 
7018 REAL tetgenmesh::insphere_sos(REAL* pa, REAL* pb, REAL* pc, REAL* pd, REAL* pe,
7019  int ia, int ib, int ic, int id, int ie)
7020 {
7021  REAL det;
7022 
7023  det = insphere(pa, pb, pc, pd, pe);
7024  if (det != 0.0) {
7025  return det;
7026  }
7027 
7028  // det = 0.0, use symbolic perturbation.
7029  REAL *p[5], *tmpp;
7030  REAL sign, det_c, det_d;
7031  int idx[5], perm, tmp;
7032  int n, i, j;
7033 
7034  p[0] = pa; idx[0] = ia;
7035  p[1] = pb; idx[1] = ib;
7036  p[2] = pc; idx[2] = ic;
7037  p[3] = pd; idx[3] = id;
7038  p[4] = pe; idx[4] = ie;
7039 
7040  // Bubble sort the points by the increasing order of the indices.
7041  n = 5;
7042  perm = 0; // The number of total swaps.
7043  for (i = 0; i < n - 1; i++) {
7044  for (j = 0; j < n - 1 - i; j++) {
7045  if (idx[j + 1] < idx[j]) { // compare the two neighbors.
7046  tmp = idx[j]; // swap idx[j] and idx[j + 1]
7047  idx[j] = idx[j + 1];
7048  idx[j + 1] = tmp;
7049  tmpp = p[j]; // swap p[j] and p[j + 1]
7050  p[j] = p[j + 1];
7051  p[j + 1] = tmpp;
7052  perm++;
7053  }
7054  }
7055  }
7056 
7057  sign = (perm % 2 == 0) ? 1.0 : -1.0;
7058  det_c = orient3d(p[1], p[2], p[3], p[4]); // orient3d(b, c, d, e)
7059  if (det_c != 0.0) {
7060  return sign * det_c;
7061  }
7062  det_d = orient3d(p[0], p[2], p[3], p[4]); // orient3d(a, c, d, e)
7063  return -sign * det_d;
7064 }
7065 
7067 // //
7068 // iscollinear() Check if three points are approximately collinear. //
7069 // //
7070 // 'eps' is a relative error tolerance. The collinearity is determined by //
7071 // the value q = cos(theta), where theta is the angle between two vectors //
7072 // A->B and A->C. They're collinear if 1.0 - q <= epspp. //
7073 // //
7075 
7076 bool tetgenmesh::iscollinear(REAL* A, REAL* B, REAL* C, REAL eps)
7077 {
7078  REAL abx, aby, abz;
7079  REAL acx, acy, acz;
7080  REAL Lv, Lw, dd;
7081  REAL d, q;
7082 
7083  // Limit of two closed points.
7084  q = longest * eps;
7085  q *= q;
7086 
7087  abx = A[0] - B[0];
7088  aby = A[1] - B[1];
7089  abz = A[2] - B[2];
7090  acx = A[0] - C[0];
7091  acy = A[1] - C[1];
7092  acz = A[2] - C[2];
7093  Lv = abx * abx + aby * aby + abz * abz;
7094  // Is AB (nearly) indentical?
7095  if (Lv < q) return true;
7096  Lw = acx * acx + acy * acy + acz * acz;
7097  // Is AC (nearly) indentical?
7098  if (Lw < q) return true;
7099  dd = abx * acx + aby * acy + abz * acz;
7100 
7101  d = (dd * dd) / (Lv * Lw);
7102  if (d > 1.0) d = 1.0; // Rounding.
7103  q = 1.0 - sqrt(d); // Notice 0 < q < 1.0.
7104 
7105  return q <= eps;
7106 }
7107 
7109 // //
7110 // iscoplanar() Check if four points are approximately coplanar. //
7111 // //
7112 // 'vol6' is six times of the signed volume of the tetrahedron formed by the //
7113 // four points. 'eps' is the relative error tolerance. The coplanarity is //
7114 // determined by the value: q = fabs(vol6) / L^3, where L is the average //
7115 // edge length of the tet. They're coplanar if q <= eps. //
7116 // //
7118 
7119 bool tetgenmesh::
7120 iscoplanar(REAL* k, REAL* l, REAL* m, REAL* n, REAL vol6, REAL eps)
7121 {
7122  REAL L, q;
7123  REAL x, y, z;
7124 
7125  if (vol6 == 0.0) return true;
7126 
7127  x = k[0] - l[0];
7128  y = k[1] - l[1];
7129  z = k[2] - l[2];
7130  L = sqrt(x * x + y * y + z * z);
7131  x = l[0] - m[0];
7132  y = l[1] - m[1];
7133  z = l[2] - m[2];
7134  L += sqrt(x * x + y * y + z * z);
7135  x = m[0] - k[0];
7136  y = m[1] - k[1];
7137  z = m[2] - k[2];
7138  L += sqrt(x * x + y * y + z * z);
7139  x = k[0] - n[0];
7140  y = k[1] - n[1];
7141  z = k[2] - n[2];
7142  L += sqrt(x * x + y * y + z * z);
7143  x = l[0] - n[0];
7144  y = l[1] - n[1];
7145  z = l[2] - n[2];
7146  L += sqrt(x * x + y * y + z * z);
7147  x = m[0] - n[0];
7148  y = m[1] - n[1];
7149  z = m[2] - n[2];
7150  L += sqrt(x * x + y * y + z * z);
7151 #ifdef SELF_CHECK
7152  assert(L > 0.0);
7153 #endif
7154  L /= 6.0;
7155  q = fabs(vol6) / (L * L * L);
7156 
7157  return q <= eps;
7158 }
7159 
7161 // //
7162 // iscospheric() Check if five points are approximately coplanar. //
7163 // //
7164 // 'vol24' is the 24 times of the signed volume of the 4-dimensional simplex //
7165 // formed by the five points. 'eps' is the relative tolerance. The cosphere //
7166 // case is determined by the value: q = fabs(vol24) / L^4, where L is the //
7167 // average edge length of the simplex. They're cosphere if q <= eps. //
7168 // //
7170 
7171 bool tetgenmesh::
7172 iscospheric(REAL* k, REAL* l, REAL* m, REAL* n, REAL* o, REAL vol24, REAL eps)
7173 {
7174  REAL L, q;
7175 
7176  // A 4D simplex has 10 edges.
7177  L = distance(k, l);
7178  L += distance(l, m);
7179  L += distance(m, k);
7180  L += distance(k, n);
7181  L += distance(l, n);
7182  L += distance(m, n);
7183  L += distance(k, o);
7184  L += distance(l, o);
7185  L += distance(m, o);
7186  L += distance(n, o);
7187 #ifdef SELF_CHECK
7188  assert(L > 0.0);
7189 #endif
7190  L /= 10.0;
7191  q = fabs(vol24) / (L * L * L * L);
7192 
7193  return q < eps;
7194 }
7195 
7196 //
7197 // End of geometric tests
7198 //
7199 
7200 //
7201 // Begin of Geometric quantities calculators
7202 //
7203 
7204 // distance() computs the Euclidean distance between two points.
7205 inline REAL tetgenmesh::distance(REAL* p1, REAL* p2)
7206 {
7207  return sqrt((p2[0] - p1[0]) * (p2[0] - p1[0]) +
7208  (p2[1] - p1[1]) * (p2[1] - p1[1]) +
7209  (p2[2] - p1[2]) * (p2[2] - p1[2]));
7210 }
7211 
7213 // //
7214 // shortdistance() Returns the shortest distance from point p to a line //
7215 // defined by two points e1 and e2. //
7216 // //
7217 // First compute the projection length l_p of the vector v1 = p - e1 along //
7218 // the vector v2 = e2 - e1. Then Pythagoras' Theorem is used to compute the //
7219 // shortest distance. //
7220 // //
7221 // This routine allows that p is collinear with the line. In this case, the //
7222 // return value is zero. The two points e1 and e2 should not be identical. //
7223 // //
7225 
7226 REAL tetgenmesh::shortdistance(REAL* p, REAL* e1, REAL* e2)
7227 {
7228  REAL v1[3], v2[3];
7229  REAL len, l_p;
7230 
7231  v1[0] = e2[0] - e1[0];
7232  v1[1] = e2[1] - e1[1];
7233  v1[2] = e2[2] - e1[2];
7234  v2[0] = p[0] - e1[0];
7235  v2[1] = p[1] - e1[1];
7236  v2[2] = p[2] - e1[2];
7237 
7238  len = sqrt(dot(v1, v1));
7239 #ifdef SELF_CHECK
7240  assert(len != 0.0);
7241 #endif
7242  v1[0] /= len;
7243  v1[1] /= len;
7244  v1[2] /= len;
7245  l_p = dot(v1, v2);
7246 
7247  return sqrt(dot(v2, v2) - l_p * l_p);
7248 }
7249 
7251 // //
7252 // shortdistance() Returns the shortest distance from point p to a face. //
7253 // //
7255 
7256 REAL tetgenmesh::shortdistance(REAL* p, REAL* e1, REAL* e2, REAL* e3)
7257 {
7258  REAL prj[3];
7259 
7260  projpt2face(p, e1, e2, e3, prj);
7261  return distance(p, prj);
7262 }
7263 
7265 // //
7266 // interiorangle() Return the interior angle (0 - 2 * PI) between vectors //
7267 // o->p1 and o->p2. //
7268 // //
7269 // 'n' is the normal of the plane containing face (o, p1, p2). The interior //
7270 // angle is the total angle rotating from o->p1 around n to o->p2. Exchange //
7271 // the position of p1 and p2 will get the complement angle of the other one. //
7272 // i.e., interiorangle(o, p1, p2) = 2 * PI - interiorangle(o, p2, p1). Set //
7273 // 'n' be NULL if you only want the interior angle between 0 - PI. //
7274 // //
7276 
7277 REAL tetgenmesh::interiorangle(REAL* o, REAL* p1, REAL* p2, REAL* n)
7278 {
7279  REAL v1[3], v2[3], np[3];
7280  REAL theta, costheta, lenlen;
7281  REAL ori, len1, len2;
7282 
7283  // Get the interior angle (0 - PI) between o->p1, and o->p2.
7284  v1[0] = p1[0] - o[0];
7285  v1[1] = p1[1] - o[1];
7286  v1[2] = p1[2] - o[2];
7287  v2[0] = p2[0] - o[0];
7288  v2[1] = p2[1] - o[1];
7289  v2[2] = p2[2] - o[2];
7290  len1 = sqrt(dot(v1, v1));
7291  len2 = sqrt(dot(v2, v2));
7292  lenlen = len1 * len2;
7293 #ifdef SELF_CHECK
7294  assert(lenlen != 0.0);
7295 #endif
7296  costheta = dot(v1, v2) / lenlen;
7297  if (costheta > 1.0) {
7298  costheta = 1.0; // Roundoff.
7299  } else if (costheta < -1.0) {
7300  costheta = -1.0; // Roundoff.
7301  }
7302  theta = acos(costheta);
7303  if (n != NULL) {
7304  // Get a point above the face (o, p1, p2);
7305  np[0] = o[0] + n[0];
7306  np[1] = o[1] + n[1];
7307  np[2] = o[2] + n[2];
7308  // Adjust theta (0 - 2 * PI).
7309  ori = orient3d(p1, o, np, p2);
7310  if (ori > 0.0) {
7311  theta = 2 * PI - theta;
7312  }
7313  }
7314 
7315  return theta;
7316 }
7317 
7319 // //
7320 // projpt2edge() Return the projection point from a point to an edge. //
7321 // //
7323 
7324 void tetgenmesh::projpt2edge(REAL* p, REAL* e1, REAL* e2, REAL* prj)
7325 {
7326  REAL v1[3], v2[3];
7327  REAL len, l_p;
7328 
7329  v1[0] = e2[0] - e1[0];
7330  v1[1] = e2[1] - e1[1];
7331  v1[2] = e2[2] - e1[2];
7332  v2[0] = p[0] - e1[0];
7333  v2[1] = p[1] - e1[1];
7334  v2[2] = p[2] - e1[2];
7335 
7336  len = sqrt(dot(v1, v1));
7337 #ifdef SELF_CHECK
7338  assert(len != 0.0);
7339 #endif
7340  v1[0] /= len;
7341  v1[1] /= len;
7342  v1[2] /= len;
7343  l_p = dot(v1, v2);
7344 
7345  prj[0] = e1[0] + l_p * v1[0];
7346  prj[1] = e1[1] + l_p * v1[1];
7347  prj[2] = e1[2] + l_p * v1[2];
7348 }
7349 
7351 // //
7352 // projpt2face() Return the projection point from a point to a face. //
7353 // //
7355 
7356 void tetgenmesh::projpt2face(REAL* p, REAL* f1, REAL* f2, REAL* f3, REAL* prj)
7357 {
7358  REAL fnormal[3], v1[3];
7359  REAL len, dist;
7360 
7361  // Get the unit face normal.
7362  facenormal(f1, f2, f3, fnormal, &len);
7363 #ifdef SELF_CHECK
7364  assert(len > 0.0);
7365 #endif
7366  fnormal[0] /= len;
7367  fnormal[1] /= len;
7368  fnormal[2] /= len;
7369  // Get the vector v1 = |p - f1|.
7370  v1[0] = p[0] - f1[0];
7371  v1[1] = p[1] - f1[1];
7372  v1[2] = p[2] - f1[2];
7373  // Get the project distance.
7374  dist = dot(fnormal, v1);
7375 
7376  // Get the project point.
7377  prj[0] = p[0] - dist * fnormal[0];
7378  prj[1] = p[1] - dist * fnormal[1];
7379  prj[2] = p[2] - dist * fnormal[2];
7380 }
7381 
7383 // //
7384 // facenormal() Calculate the normal of a face given by three points. //
7385 // //
7386 // In general, the face normal can be calculate by the cross product of any //
7387 // pair of the three edge vectors. However, if the three points are nearly //
7388 // collinear, the rounding error may harm the result. To choose a good pair //
7389 // of vectors is helpful to reduce the error. //
7390 // //
7392 
7393 void tetgenmesh::facenormal(REAL* pa, REAL* pb, REAL* pc, REAL* n, REAL* nlen)
7394 {
7395  REAL v1[3], v2[3];
7396 
7397  v1[0] = pb[0] - pa[0];
7398  v1[1] = pb[1] - pa[1];
7399  v1[2] = pb[2] - pa[2];
7400  v2[0] = pc[0] - pa[0];
7401  v2[1] = pc[1] - pa[1];
7402  v2[2] = pc[2] - pa[2];
7403 
7404  cross(v1, v2, n);
7405  if (nlen != (REAL *) NULL) {
7406  *nlen = sqrt(dot(n, n));
7407  }
7408 }
7409 
7411 // //
7412 // edgeorthonormal() Return the unit normal of an edge in a given plane. //
7413 // //
7414 // The edge is from e1 to e2, the plane is defined by given an additional //
7415 // point op, which is non-collinear with the edge. In addition, the side of //
7416 // the edge in which op lies defines the positive position of the normal. //
7417 // //
7418 // Let v1 be the unit vector from e1 to e2, v2 be the unit edge vector from //
7419 // e1 to op, fn be the unit face normal calculated by fn = v1 x v2. Then the //
7420 // unit edge normal of e1e2 pointing to op is n = fn x v1. Note, we should //
7421 // not change the position of fn and v1, otherwise, we get the edge normal //
7422 // pointing to the other side of op. //
7423 // //
7425 
7426 void tetgenmesh::edgeorthonormal(REAL* e1, REAL* e2, REAL* op, REAL* n)
7427 {
7428  REAL v1[3], v2[3], fn[3];
7429  REAL len;
7430 
7431  // Get the edge vector v1.
7432  v1[0] = e2[0] - e1[0];
7433  v1[1] = e2[1] - e1[1];
7434  v1[2] = e2[2] - e1[2];
7435  // Get the edge vector v2.
7436  v2[0] = op[0] - e1[0];
7437  v2[1] = op[1] - e1[1];
7438  v2[2] = op[2] - e1[2];
7439  // Get the face normal fn = v1 x v2.
7440  cross(v1, v2, fn);
7441  // Get the edge normal n pointing to op. n = fn x v1.
7442  cross(fn, v1, n);
7443  // Normalize the vector.
7444  len = sqrt(dot(n, n));
7445  n[0] /= len;
7446  n[1] /= len;
7447  n[2] /= len;
7448 }
7449 
7451 // //
7452 // facedihedral() Return the dihedral angle (in radian) between two //
7453 // adjoining faces. //
7454 // //
7455 // 'pa', 'pb' are the shared edge of these two faces, 'pc1', and 'pc2' are //
7456 // apexes of these two faces. Return the angle (between 0 to 2*pi) between //
7457 // the normal of face (pa, pb, pc1) and normal of face (pa, pb, pc2). //
7458 // //
7460 
7461 REAL tetgenmesh::facedihedral(REAL* pa, REAL* pb, REAL* pc1, REAL* pc2)
7462 {
7463  REAL n1[3], n2[3];
7464  REAL n1len, n2len;
7465  REAL costheta, ori;
7466  REAL theta;
7467 
7468  facenormal(pa, pb, pc1, n1, &n1len);
7469  facenormal(pa, pb, pc2, n2, &n2len);
7470  costheta = dot(n1, n2) / (n1len * n2len);
7471  // Be careful rounding error!
7472  if (costheta > 1.0) {
7473  costheta = 1.0;
7474  } else if (costheta < -1.0) {
7475  costheta = -1.0;
7476  }
7477  theta = acos(costheta);
7478  ori = orient3d(pa, pb, pc1, pc2);
7479  if (ori > 0.0) {
7480  theta = 2 * PI - theta;
7481  }
7482 
7483  return theta;
7484 }
7485 
7487 // //
7488 // tetalldihedral() Get all (six) dihedral angles of a tet. //
7489 // //
7490 // The tet is given by its four corners a, b, c, and d. If 'cosdd' is not //
7491 // NULL, it returns the cosines of the 6 dihedral angles, the corresponding //
7492 // edges are: ab, bc, ca, ad, bd, and cd. If 'cosmaxd' (or 'cosmind') is not //
7493 // NULL, it returns the cosine of the maximal (or minimal) dihedral angle. //
7494 // //
7496 
7497 void tetgenmesh::tetalldihedral(point pa, point pb, point pc, point pd,
7498  REAL* cosdd, REAL* cosmaxd, REAL* cosmind)
7499 {
7500  REAL N[4][3], cosd, len;
7501  int f1, f2, i, j;
7502 
7503  f1=0;
7504  f2=0;
7505 
7506  // Get four normals of faces of the tet.
7507  tetallnormal(pa, pb, pc, pd, N, NULL);
7508  // Normalize the normals.
7509  for (i = 0; i < 4; i++) {
7510  len = sqrt(dot(N[i], N[i]));
7511  if (len != 0.0) {
7512  for (j = 0; j < 3; j++) N[i][j] /= len;
7513  }
7514  }
7515 
7516  for (i = 0; i < 6; i++) {
7517  switch (i) {
7518  case 0: f1 = 2; f2 = 3; break; // edge ab.
7519  case 1: f1 = 0; f2 = 3; break; // edge bc.
7520  case 2: f1 = 1; f2 = 3; break; // edge ca.
7521  case 3: f1 = 1; f2 = 2; break; // edge ad.
7522  case 4: f1 = 2; f2 = 0; break; // edge bd.
7523  case 5: f1 = 0; f2 = 1; break; // edge cd.
7524  }
7525  cosd = -dot(N[f1], N[f2]);
7526  if (cosdd) cosdd[i] = cosd;
7527  if (i == 0) {
7528  if (cosmaxd) *cosmaxd = cosd;
7529  if (cosmind) *cosmind = cosd;
7530  } else {
7531  if (cosmaxd) *cosmaxd = cosd < *cosmaxd ? cosd : *cosmaxd;
7532  if (cosmind) *cosmind = cosd > *cosmind ? cosd : *cosmind;
7533  }
7534  }
7535 }
7536 
7538 // //
7539 // tetallnormal() Get the in-noramls of the four faces of a given tet. //
7540 // //
7541 // Let tet be abcd. N[4][3] returns the four normals, which are: N[0] cbd, //
7542 // N[1] acd, N[2] bad, N[3] abc. These normals are unnormalized. //
7543 // //
7545 
7546 void tetgenmesh::tetallnormal(point pa, point pb, point pc, point pd,
7547  REAL N[4][3], REAL* volume)
7548 {
7549  REAL A[4][4], rhs[4], D;
7550  int indx[4];
7551  int i, j;
7552 
7553  // get the entries of A[3][3].
7554  for (i = 0; i < 3; i++) A[0][i] = pa[i] - pd[i]; // d->a vec
7555  for (i = 0; i < 3; i++) A[1][i] = pb[i] - pd[i]; // d->b vec
7556  for (i = 0; i < 3; i++) A[2][i] = pc[i] - pd[i]; // d->c vec
7557  // Compute the inverse of matrix A, to get 3 normals of the 4 faces.
7558  lu_decmp(A, 3, indx, &D, 0); // Decompose the matrix just once.
7559  if (volume != NULL) {
7560  // Get the volume of the tet.
7561  *volume = fabs((A[indx[0]][0] * A[indx[1]][1] * A[indx[2]][2])) / 6.0;
7562  }
7563  for (j = 0; j < 3; j++) {
7564  for (i = 0; i < 3; i++) rhs[i] = 0.0;
7565  rhs[j] = 1.0; // Positive means the inside direction
7566  lu_solve(A, 3, indx, rhs, 0);
7567  for (i = 0; i < 3; i++) N[j][i] = rhs[i];
7568  }
7569  // Get the fourth normal by summing up the first three.
7570  for (i = 0; i < 3; i++) N[3][i] = - N[0][i] - N[1][i] - N[2][i];