9 #define UNUSED_OPT(var) var=var
11 #define UNUSED_OPT(var)
58 void terminatetetgen(
int x)
64 #endif // #ifdef TETLIBRARY
84 void tetgenio::initialize()
90 pointlist = (REAL *) NULL;
91 pointattributelist = (REAL *) NULL;
92 pointmtrlist = (REAL *) NULL;
93 pointmarkerlist = (
int *) NULL;
95 numberofpointattributes = 0;
96 numberofpointmtrs = 0;
98 tetrahedronlist = (
int *) NULL;
99 tetrahedronattributelist = (REAL *) NULL;
100 tetrahedronvolumelist = (REAL *) NULL;
101 neighborlist = (
int *) NULL;
102 numberoftetrahedra = 0;
104 numberoftetrahedronattributes = 0;
106 trifacelist = (
int *) NULL;
107 adjtetlist = (
int *) NULL;
108 trifacemarkerlist = (
int *) NULL;
109 numberoftrifaces = 0;
111 facetlist = (facet *) NULL;
112 facetmarkerlist = (
int *) NULL;
115 edgelist = (
int *) NULL;
116 edgemarkerlist = (
int *) NULL;
119 holelist = (REAL *) NULL;
122 regionlist = (REAL *) NULL;
125 facetconstraintlist = (REAL *) NULL;
126 numberoffacetconstraints = 0;
127 segmentconstraintlist = (REAL *) NULL;
128 numberofsegmentconstraints = 0;
130 pbcgrouplist = (pbcgroup *) NULL;
131 numberofpbcgroups = 0;
133 vpointlist = (REAL *) NULL;
134 vedgelist = (voroedge *) NULL;
135 vfacetlist = (vorofacet *) NULL;
136 vcelllist = (
int **) NULL;
158 void tetgenio::deinitialize()
165 if (pointlist != (REAL *) NULL) {
168 if (pointattributelist != (REAL *) NULL) {
169 delete [] pointattributelist;
171 if (pointmtrlist != (REAL *) NULL) {
172 delete [] pointmtrlist;
174 if (pointmarkerlist != (
int *) NULL) {
175 delete [] pointmarkerlist;
178 if (tetrahedronlist != (
int *) NULL) {
179 delete [] tetrahedronlist;
181 if (tetrahedronattributelist != (REAL *) NULL) {
182 delete [] tetrahedronattributelist;
184 if (tetrahedronvolumelist != (REAL *) NULL) {
185 delete [] tetrahedronvolumelist;
187 if (neighborlist != (
int *) NULL) {
188 delete [] neighborlist;
191 if (trifacelist != (
int *) NULL) {
192 delete [] trifacelist;
194 if (adjtetlist != (
int *) NULL) {
195 delete [] adjtetlist;
197 if (trifacemarkerlist != (
int *) NULL) {
198 delete [] trifacemarkerlist;
201 if (edgelist != (
int *) NULL) {
204 if (edgemarkerlist != (
int *) NULL) {
205 delete [] edgemarkerlist;
208 if (facetlist != (facet *) NULL) {
209 for (i = 0; i < numberoffacets; i++) {
211 for (j = 0; j < f->numberofpolygons; j++) {
212 p = &f->polygonlist[j];
213 delete [] p->vertexlist;
215 delete [] f->polygonlist;
216 if (f->holelist != (REAL *) NULL) {
217 delete [] f->holelist;
222 if (facetmarkerlist != (
int *) NULL) {
223 delete [] facetmarkerlist;
226 if (holelist != (REAL *) NULL) {
229 if (regionlist != (REAL *) NULL) {
230 delete [] regionlist;
232 if (facetconstraintlist != (REAL *) NULL) {
233 delete [] facetconstraintlist;
235 if (segmentconstraintlist != (REAL *) NULL) {
236 delete [] segmentconstraintlist;
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;
245 delete [] pbcgrouplist;
247 if (vpointlist != (REAL *) NULL) {
248 delete [] vpointlist;
250 if (vedgelist != (voroedge *) NULL) {
253 if (vfacetlist != (vorofacet *) NULL) {
254 for (i = 0; i < numberofvfacets; i++) {
255 delete [] vfacetlist[i].elist;
257 delete [] vfacetlist;
259 if (vcelllist != (
int **) NULL) {
260 for (i = 0; i < numberofvcells; i++) {
261 delete [] vcelllist[i];
282 bool tetgenio::load_node_call(FILE* infile,
int markers,
char* infilename)
284 char inputline[INPUTLINESIZE];
286 REAL x, y, z, attrib;
287 int firstnode, currentmarker;
288 int index, attribindex;
292 pointlist =
new REAL[numberofpoints * 3];
293 if (pointlist == (REAL *) NULL) {
294 printf(
"Error: Out of memory.\n");
297 if (numberofpointattributes > 0) {
298 pointattributelist =
new REAL[numberofpoints * numberofpointattributes];
299 if (pointattributelist == (REAL *) NULL) {
300 printf(
"Error: Out of memory.\n");
305 pointmarkerlist =
new int[numberofpoints];
306 if (pointmarkerlist == (
int *) NULL) {
307 printf(
"Error: Out of memory.\n");
315 for (i = 0; i < numberofpoints; i++) {
316 stringptr = readnumberline(inputline, infile, infilename);
319 firstnode = (int) strtol (stringptr, &stringptr, 0);
320 if ((firstnode == 0) || (firstnode == 1)) {
321 firstnumber = firstnode;
324 stringptr = findnextnumber(stringptr);
326 if (*stringptr ==
'\0') {
327 printf(
"Error: Point %d has no x coordinate.\n", firstnumber + i);
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);
336 y = (REAL) strtod(stringptr, &stringptr);
338 stringptr = findnextnumber(stringptr);
339 if (*stringptr ==
'\0') {
340 printf(
"Error: Point %d has no z coordinate.\n", firstnumber + i);
343 z = (REAL) strtod(stringptr, &stringptr);
347 pointlist[index++] = x;
348 pointlist[index++] = y;
349 pointlist[index++] = z;
351 for (j = 0; j < numberofpointattributes; j++) {
352 stringptr = findnextnumber(stringptr);
353 if (*stringptr ==
'\0') {
356 attrib = (REAL) strtod(stringptr, &stringptr);
358 pointattributelist[attribindex++] = attrib;
362 stringptr = findnextnumber(stringptr);
363 if (*stringptr ==
'\0') {
366 currentmarker = (int) strtol (stringptr, &stringptr, 0);
368 pointmarkerlist[i] = currentmarker;
371 if (i < numberofpoints) {
374 pointlist = (REAL *) NULL;
376 delete [] pointmarkerlist;
377 pointmarkerlist = (
int *) NULL;
379 if (numberofpointattributes > 0) {
380 delete [] pointattributelist;
381 pointattributelist = (REAL *) NULL;
398 bool tetgenio::load_node(
char* filename)
401 char innodefilename[FILENAMESIZE];
402 char inputline[INPUTLINESIZE];
408 strcpy(innodefilename, filename);
409 strcat(innodefilename,
".node");
412 infile = fopen(innodefilename,
"r");
413 if (infile == (FILE *) NULL) {
414 printf(
"File I/O Error: Cannot access file %s.\n", innodefilename);
417 printf(
"Opening %s.\n", innodefilename);
419 stringptr = readnumberline(inputline, infile, innodefilename);
421 stringptr = strstr(inputline,
"rbox");
422 if (stringptr == NULL) {
425 stringptr = inputline;
426 numberofpoints = (int) strtol (stringptr, &stringptr, 0);
427 stringptr = findnextnumber(stringptr);
428 if (*stringptr ==
'\0') {
431 mesh_dim = (int) strtol (stringptr, &stringptr, 0);
433 stringptr = findnextnumber(stringptr);
434 if (*stringptr ==
'\0') {
435 numberofpointattributes = 0;
437 numberofpointattributes = (int) strtol (stringptr, &stringptr, 0);
439 stringptr = findnextnumber(stringptr);
440 if (*stringptr ==
'\0') {
443 markers = (int) strtol (stringptr, &stringptr, 0);
447 stringptr = inputline;
449 mesh_dim = (int) strtol (stringptr, &stringptr, 0);
451 stringptr = readnumberline(inputline, infile, innodefilename);
452 numberofpoints = (int) strtol (stringptr, &stringptr, 0);
462 if (numberofpoints < (mesh_dim + 1)) {
463 printf(
"Input error: TetGen needs at least %d points.\n", mesh_dim + 1);
469 if (!load_node_call(infile, markers, innodefilename)) {
488 bool tetgenio::load_pbc(
char* filename)
491 char pbcfilename[FILENAMESIZE];
492 char inputline[INPUTLINESIZE];
499 strcpy(pbcfilename, filename);
500 strcat(pbcfilename,
".pbc");
501 infile = fopen(pbcfilename,
"r");
502 if (infile != (FILE *) NULL) {
503 printf(
"Opening %s.\n", pbcfilename);
510 stringptr = readnumberline(inputline, infile, pbcfilename);
511 numberofpbcgroups = (int) strtol (stringptr, &stringptr, 0);
512 if (numberofpbcgroups == 0) {
518 pbcgrouplist =
new pbcgroup[numberofpbcgroups];
521 for (i = 0; i < numberofpbcgroups; i++) {
522 pg = &(pbcgrouplist[i]);
524 pg->numberofpointpairs = 0;
525 pg->pointpairlist = (
int *) NULL;
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);
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++) {
541 stringptr = findnextnumber(stringptr);
542 if (*stringptr ==
'\0') {
544 stringptr = readnumberline(inputline, infile, pbcfilename);
545 if (*stringptr ==
'\0')
break;
547 pg->transmat[j][k] = (REAL) strtod(stringptr, &stringptr);
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];
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;
572 if (i < numberofpbcgroups) {
574 delete [] pbcgrouplist;
575 pbcgrouplist = (pbcgroup *) NULL;
576 numberofpbcgroups = 0;
591 bool tetgenio::load_var(
char* filename)
594 char varfilename[FILENAMESIZE];
595 char inputline[INPUTLINESIZE];
601 strcpy(varfilename, filename);
602 strcat(varfilename,
".var");
603 infile = fopen(varfilename,
"r");
604 if (infile != (FILE *) NULL) {
605 printf(
"Opening %s.\n", varfilename);
612 stringptr = readnumberline(inputline, infile, varfilename);
613 if (*stringptr !=
'\0') {
614 numberoffacetconstraints = (int) strtol (stringptr, &stringptr, 0);
616 numberoffacetconstraints = 0;
618 if (numberoffacetconstraints > 0) {
620 facetconstraintlist =
new REAL[numberoffacetconstraints * 2];
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",
630 facetconstraintlist[index++] = (REAL) strtod(stringptr, &stringptr);
632 stringptr = findnextnumber(stringptr);
633 if (*stringptr ==
'\0') {
634 printf(
"Error: facet constraint %d has no maximum area bound.\n",
638 facetconstraintlist[index++] = (REAL) strtod(stringptr, &stringptr);
641 if (i < numberoffacetconstraints) {
649 stringptr = readnumberline(inputline, infile, varfilename);
650 if (*stringptr !=
'\0') {
651 numberofsegmentconstraints = (int) strtol (stringptr, &stringptr, 0);
653 numberofsegmentconstraints = 0;
655 if (numberofsegmentconstraints > 0) {
657 segmentconstraintlist =
new REAL[numberofsegmentconstraints * 3];
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",
667 segmentconstraintlist[index++] = (REAL) strtod(stringptr, &stringptr);
669 stringptr = findnextnumber(stringptr);
670 if (*stringptr ==
'\0') {
671 printf(
"Error: segment constraint %d has no second endpoint.\n",
675 segmentconstraintlist[index++] = (REAL) strtod(stringptr, &stringptr);
677 stringptr = findnextnumber(stringptr);
678 if (*stringptr ==
'\0') {
679 printf(
"Error: segment constraint %d has no maximum length bound.\n",
683 segmentconstraintlist[index++] = (REAL) strtod(stringptr, &stringptr);
686 if (i < numberofsegmentconstraints) {
706 bool tetgenio::load_mtr(
char* filename)
709 char mtrfilename[FILENAMESIZE];
710 char inputline[INPUTLINESIZE];
716 strcpy(mtrfilename, filename);
717 strcat(mtrfilename,
".mtr");
718 infile = fopen(mtrfilename,
"r");
719 if (infile != (FILE *) NULL) {
720 printf(
"Opening %s.\n", mtrfilename);
727 stringptr = readnumberline(inputline, infile, mtrfilename);
728 stringptr = findnextnumber(stringptr);
729 if (*stringptr !=
'\0') {
730 numberofpointmtrs = (int) strtol (stringptr, &stringptr, 0);
732 if (numberofpointmtrs == 0) {
734 numberofpointmtrs = 1;
738 pointmtrlist =
new REAL[numberofpoints * numberofpointmtrs];
739 if (pointmtrlist == (REAL *) NULL) {
740 printf(
"Error: Out of memory.\n");
744 for (i = 0; i < numberofpoints; i++) {
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);
753 mtr = (REAL) strtod(stringptr, &stringptr);
754 pointmtrlist[mtrindex++] = mtr;
755 stringptr = findnextnumber(stringptr);
773 bool tetgenio::load_poly(
char* filename)
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;
786 strcpy(innodefilename, filename);
787 strcpy(inpolyfilename, filename);
788 strcpy(insmeshfilename, filename);
789 strcat(innodefilename,
".node");
790 strcat(inpolyfilename,
".poly");
791 strcat(insmeshfilename,
".smesh");
796 polyfile = fopen(inpolyfilename,
"r");
797 if (polyfile == (FILE *) NULL) {
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);
805 printf(
"Opening %s.\n", insmeshfilename);
809 printf(
"Opening %s.\n", inpolyfilename);
813 numberofpointattributes = 0;
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);
823 stringptr = findnextnumber(stringptr);
824 if (*stringptr !=
'\0') {
825 numberofpointattributes = (int) strtol (stringptr, &stringptr, 0);
827 stringptr = findnextnumber(stringptr);
828 if (*stringptr !=
'\0') {
829 markers = (int) strtol (stringptr, &stringptr, 0);
831 if (numberofpoints > 0) {
834 infilename = insmeshfilename;
836 infilename = inpolyfilename;
843 infilename = innodefilename;
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);
856 numberofpointattributes = 0;
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);
866 stringptr = findnextnumber(stringptr);
867 if (*stringptr !=
'\0') {
868 numberofpointattributes = (int) strtol (stringptr, &stringptr, 0);
870 stringptr = findnextnumber(stringptr);
871 if (*stringptr !=
'\0') {
872 markers = (int) strtol (stringptr, &stringptr, 0);
876 if ((mesh_dim != 3) && (mesh_dim != 2)) {
877 printf(
"Input error: TetGen only works for 2D & 3D point sets.\n");
881 if (numberofpoints < (mesh_dim + 1)) {
882 printf(
"Input error: TetGen needs at least %d points.\n", mesh_dim + 1);
888 if (!load_node_call(infile, markers, infilename)) {
903 stringptr = readnumberline(inputline, polyfile, inpolyfilename);
904 numberoffacets = (int) strtol (stringptr, &stringptr, 0);
905 if (numberoffacets <= 0) {
910 stringptr = findnextnumber(stringptr);
911 if (*stringptr ==
'\0') {
914 markers = (int) strtol (stringptr, &stringptr, 0);
918 facetlist =
new facet[numberoffacets];
920 facetmarkerlist =
new int[numberoffacets];
926 for (i = 1; i <= numberoffacets; i++) {
927 f = &(facetlist[i - 1]);
929 f->numberofholes = 0;
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);
938 stringptr = findnextnumber(stringptr);
939 if (*stringptr !=
'\0') {
940 currentmarker = (int) strtol(stringptr, &stringptr, 0);
946 facetmarkerlist[i - 1] = currentmarker;
949 if (f->numberofpolygons <= 0) {
950 printf(
"Error: Wrong number of polygon in %d facet.\n", i);
954 f->polygonlist =
new polygon[f->numberofpolygons];
956 for (j = 1; j <= f->numberofpolygons; j++) {
957 p = &(f->polygonlist[j - 1]);
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);
967 p->vertexlist =
new int[p->numberofvertices];
969 for (k = 1; k <= p->numberofvertices; k++) {
970 stringptr = findnextnumber(stringptr);
971 if (*stringptr ==
'\0') {
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);
981 p->vertexlist[k - 1] = (int) strtol (stringptr, &stringptr, 0);
984 if (j <= f->numberofpolygons) {
989 delete [] f->polygonlist;
991 f->numberofpolygons = j - 1;
993 f->numberofholes = 0;
997 if (f->numberofholes > 0) {
999 f->holelist =
new REAL[f->numberofholes * 3];
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);
1010 f->holelist[index++] = (REAL) strtod (stringptr, &stringptr);
1017 if (j <= f->numberofholes) {
1023 if (i <= numberoffacets) {
1025 numberoffacets = i - 1;
1031 for (i = 1; i <= numberoffacets; i++) {
1032 f = &(facetlist[i - 1]);
1036 f->numberofpolygons = 1;
1037 f->polygonlist =
new polygon[f->numberofpolygons];
1038 p = &(f->polygonlist[0]);
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);
1048 p->vertexlist =
new int[p->numberofvertices];
1049 for (k = 1; k <= p->numberofvertices; k++) {
1050 stringptr = findnextnumber(stringptr);
1051 if (*stringptr ==
'\0') {
1054 stringptr = readnumberline(inputline, polyfile, inpolyfilename);
1055 if (*stringptr ==
'\0') {
1056 printf(
"Error: Missing %d endpoints in facet %d",
1057 p->numberofvertices - k, i);
1061 p->vertexlist[k - 1] = (int) strtol (stringptr, &stringptr, 0);
1063 if (k <= p->numberofvertices) {
1069 stringptr = findnextnumber(stringptr);
1070 if (*stringptr ==
'\0') {
1073 currentmarker = (int) strtol(stringptr, &stringptr, 0);
1075 facetmarkerlist[i - 1] = currentmarker;
1078 if (i <= numberoffacets) {
1080 numberoffacets = i - 1;
1087 stringptr = readnumberline(inputline, polyfile, inpolyfilename);
1088 if (*stringptr !=
'\0') {
1089 numberofholes = (int) strtol (stringptr, &stringptr, 0);
1093 if (numberofholes > 0) {
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));
1103 holelist[i] = (REAL) strtod(stringptr, &stringptr);
1105 stringptr = findnextnumber(stringptr);
1106 if (*stringptr ==
'\0') {
1107 printf(
"Error: Hole %d has no y coord.\n", firstnumber + (i / 3));
1110 holelist[i + 1] = (REAL) strtod(stringptr, &stringptr);
1112 stringptr = findnextnumber(stringptr);
1113 if (*stringptr ==
'\0') {
1114 printf(
"Error: Hole %d has no z coord.\n", firstnumber + (i / 3));
1117 holelist[i + 2] = (REAL) strtod(stringptr, &stringptr);
1120 if (i < 3 * numberofholes) {
1129 stringptr = readnumberline(inputline, polyfile, NULL);
1130 if (stringptr != (
char *) NULL && *stringptr !=
'\0') {
1131 numberofregions = (int) strtol (stringptr, &stringptr, 0);
1133 numberofregions = 0;
1135 if (numberofregions > 0) {
1137 regionlist =
new REAL[numberofregions * 5];
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);
1146 regionlist[index++] = (REAL) strtod(stringptr, &stringptr);
1148 stringptr = findnextnumber(stringptr);
1149 if (*stringptr ==
'\0') {
1150 printf(
"Error: Region %d has no y coordinate.\n", firstnumber + i);
1153 regionlist[index++] = (REAL) strtod(stringptr, &stringptr);
1155 stringptr = findnextnumber(stringptr);
1156 if (*stringptr ==
'\0') {
1157 printf(
"Error: Region %d has no z coordinate.\n", firstnumber + i);
1160 regionlist[index++] = (REAL) strtod(stringptr, &stringptr);
1162 stringptr = findnextnumber(stringptr);
1163 if (*stringptr ==
'\0') {
1164 printf(
"Error: Region %d has no region attrib.\n", firstnumber + i);
1167 regionlist[index++] = (REAL) strtod(stringptr, &stringptr);
1169 stringptr = findnextnumber(stringptr);
1170 if (*stringptr ==
'\0') {
1171 regionlist[index] = regionlist[index - 1];
1173 regionlist[index] = (REAL) strtod(stringptr, &stringptr);
1177 if (i < numberofregions) {
1187 assert(mesh_dim == 2);
1191 facetlist =
new facet[numberoffacets];
1192 facetmarkerlist = (
int *) NULL;
1193 f = &(facetlist[0]);
1196 stringptr = readnumberline(inputline, polyfile, inpolyfilename);
1198 f->numberofpolygons = (int) strtol (stringptr, &stringptr, 0);
1199 if (f->numberofpolygons > 0) {
1200 f->polygonlist =
new polygon[f->numberofpolygons];
1203 for (j = 0; j < f->numberofpolygons; j++) {
1204 p = &(f->polygonlist[j]);
1207 stringptr = readnumberline(inputline, polyfile, inpolyfilename);
1208 stringptr = findnextnumber(stringptr);
1209 p->numberofvertices = 2;
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);
1216 stringptr = readnumberline(inputline, polyfile, inpolyfilename);
1217 f->numberofholes = (int) strtol (stringptr, &stringptr, 0);
1218 if (f->numberofholes > 0) {
1220 f->holelist =
new REAL[f->numberofholes * 3];
1222 for (j = 0; j < f->numberofholes; j++) {
1224 stringptr = readnumberline(inputline, polyfile, inpolyfilename);
1225 stringptr = findnextnumber(stringptr);
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;
1263 bool tetgenio::load_off(
char* filename)
1267 tetgenio::polygon *p;
1268 char infilename[FILENAMESIZE];
1269 char buffer[INPUTLINESIZE];
1272 int nverts = 0, iverts = 0;
1273 int nfaces = 0, ifaces = 0;
1275 int line_count = 0, i;
1277 strncpy(infilename, filename, 1024 - 1);
1278 infilename[FILENAMESIZE - 1] =
'\0';
1279 if (infilename[0] ==
'\0') {
1280 printf(
"Error: No filename.\n");
1283 if (strcmp(&infilename[strlen(infilename) - 4],
".off") != 0) {
1284 strcat(infilename,
".off");
1287 if (!(fp = fopen(infilename,
"r"))) {
1288 printf(
"File I/O Error: Unable to open file %s\n", infilename);
1291 printf(
"Opening %s.\n", infilename);
1296 while ((bufferp = readline(buffer, fp, &line_count)) != NULL) {
1300 bufferp = strstr(bufferp,
"OFF");
1301 if (bufferp != NULL) {
1303 bufferp = findnextnumber(bufferp);
1304 if (*bufferp ==
'\0') {
1306 bufferp = readline(buffer, fp, &line_count);
1308 if ((sscanf(bufferp,
"%d%d%d", &nverts, &nfaces, &nedges) != 3)
1310 printf(
"Syntax error reading header on line %d in file %s\n",
1311 line_count, infilename);
1317 numberofpoints = nverts;
1318 pointlist =
new REAL[nverts * 3];
1321 numberoffacets = nfaces;
1322 facetlist =
new tetgenio::facet[nfaces];
1325 }
else if (iverts < nverts) {
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);
1335 coord[i] = (REAL) strtod(bufferp, &bufferp);
1336 bufferp = findnextnumber(bufferp);
1339 }
else if (ifaces < nfaces) {
1341 f = &facetlist[ifaces];
1344 f->numberofpolygons = 1;
1345 f->polygonlist =
new tetgenio::polygon[1];
1346 p = &f->polygonlist[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);
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);
1366 p->vertexlist[i] = (int) strtol(bufferp, &bufferp, 0);
1371 printf(
"Found extra text starting at line %d in file %s\n", line_count,
1381 if (iverts != nverts) {
1382 printf(
"Expected %d vertices, but read only %d vertices in file %s\n",
1383 nverts, iverts, infilename);
1388 if (ifaces != nfaces) {
1389 printf(
"Expected %d faces, but read only %d faces in file %s\n",
1390 nfaces, ifaces, infilename);
1416 bool tetgenio::load_ply(
char* filename)
1420 tetgenio::polygon *p;
1421 char infilename[FILENAMESIZE];
1422 char buffer[INPUTLINESIZE];
1423 char *bufferp, *str;
1425 int endheader = 0, format = 0;
1426 int nverts = 0, iverts = 0;
1427 int nfaces = 0, ifaces = 0;
1428 int line_count = 0, i;
1430 strncpy(infilename, filename, FILENAMESIZE - 1);
1431 infilename[FILENAMESIZE - 1] =
'\0';
1432 if (infilename[0] ==
'\0') {
1433 printf(
"Error: No filename.\n");
1436 if (strcmp(&infilename[strlen(infilename) - 4],
".ply") != 0) {
1437 strcat(infilename,
".ply");
1440 if (!(fp = fopen(infilename,
"r"))) {
1441 printf(
"Error: Unable to open file %s\n", infilename);
1444 printf(
"Opening %s.\n", infilename);
1449 while ((bufferp = readline(buffer, fp, &line_count)) != NULL) {
1452 str = strstr(bufferp,
"end_header");
1454 if (!str) str = strstr(bufferp,
"End_header");
1455 if (!str) str = strstr(bufferp,
"End_Header");
1462 if (nverts == 0 || nfaces == 0) {
1464 str = strstr(bufferp,
"element");
1465 if (!str) str = strstr(bufferp,
"Element");
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);
1476 str = strstr(bufferp,
"vertex");
1477 if (!str) str = strstr(bufferp,
"Vertex");
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);
1486 nverts = (int) strtol(bufferp, &bufferp, 0);
1489 numberofpoints = nverts;
1490 pointlist =
new REAL[nverts * 3];
1496 str = strstr(bufferp,
"face");
1497 if (!str) str = strstr(bufferp,
"Face");
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);
1506 nfaces = (int) strtol(bufferp, &bufferp, 0);
1509 numberoffacets = nfaces;
1510 facetlist =
new tetgenio::facet[nfaces];
1518 str = strstr(bufferp,
"format");
1519 if (!str) str = strstr(bufferp,
"Format");
1522 bufferp = findnextfield(str);
1524 str = strstr(bufferp,
"ascii");
1525 if (!str) str = strstr(bufferp,
"ASCII");
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);
1536 }
else if (iverts < nverts) {
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);
1546 coord[i] = (REAL) strtod(bufferp, &bufferp);
1547 bufferp = findnextnumber(bufferp);
1550 }
else if (ifaces < nfaces) {
1552 f = &facetlist[ifaces];
1555 f->numberofpolygons = 1;
1556 f->polygonlist =
new tetgenio::polygon[1];
1557 p = &f->polygonlist[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);
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);
1577 p->vertexlist[i] = (int) strtol(bufferp, &bufferp, 0);
1582 printf(
"Found extra text starting at line %d in file %s\n", line_count,
1592 if (iverts != nverts) {
1593 printf(
"Expected %d vertices, but read only %d vertices in file %s\n",
1594 nverts, iverts, infilename);
1599 if (ifaces != nfaces) {
1600 printf(
"Expected %d faces, but read only %d faces in file %s\n",
1601 nfaces, ifaces, infilename);
1626 bool tetgenio::load_stl(
char* filename)
1629 tetgenmesh::list *plist;
1631 tetgenio::polygon *p;
1632 char infilename[FILENAMESIZE];
1633 char buffer[INPUTLINESIZE];
1634 char *bufferp, *str;
1637 int nverts = 0, iverts = 0;
1639 int line_count = 0, i;
1641 strncpy(infilename, filename, FILENAMESIZE - 1);
1642 infilename[FILENAMESIZE - 1] =
'\0';
1643 if (infilename[0] ==
'\0') {
1644 printf(
"Error: No filename.\n");
1647 if (strcmp(&infilename[strlen(infilename) - 4],
".stl") != 0) {
1648 strcat(infilename,
".stl");
1651 if (!(fp = fopen(infilename,
"r"))) {
1652 printf(
"Error: Unable to open file %s\n", infilename);
1655 printf(
"Opening %s.\n", infilename);
1658 plist =
new tetgenmesh::list(
sizeof(
double) * 3, NULL, 1024);
1660 while ((bufferp = readline(buffer, fp, &line_count)) != NULL) {
1665 bufferp = strstr(bufferp,
"solid");
1666 if (bufferp != NULL) {
1673 bufferp = strstr(bufferp,
"endsolid");
1674 if (bufferp != NULL) {
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",
1691 coord[i] = (REAL) strtod(bufferp, &bufferp);
1699 nverts = plist->len();
1701 if (nverts == 0 || (nverts % 3 != 0)) {
1702 printf(
"Error: Wrong number of vertices in file %s.\n", infilename);
1706 numberofpoints = nverts;
1707 pointlist =
new REAL[nverts * 3];
1708 for (i = 0; i < nverts; i++) {
1709 coord = (
double *) (* plist)[i];
1711 pointlist[iverts] = (REAL) coord[0];
1712 pointlist[iverts + 1] = (REAL) coord[1];
1713 pointlist[iverts + 2] = (REAL) coord[2];
1716 nfaces = (int) (nverts / 3);
1717 numberoffacets = nfaces;
1718 facetlist =
new tetgenio::facet[nfaces];
1722 iverts = firstnumber;
1723 for (i = 0; i < nfaces; i++) {
1727 f->numberofpolygons = 1;
1728 f->polygonlist =
new tetgenio::polygon[1];
1729 p = &f->polygonlist[0];
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;
1757 bool tetgenio::load_medit(
char* filename)
1760 tetgenio::facet *tmpflist, *f;
1761 tetgenio::polygon *p;
1762 char infilename[FILENAMESIZE];
1763 char buffer[INPUTLINESIZE];
1764 char *bufferp, *str;
1774 strncpy(infilename, filename, FILENAMESIZE - 1);
1775 infilename[FILENAMESIZE - 1] =
'\0';
1776 if (infilename[0] ==
'\0') {
1777 printf(
"Error: No filename.\n");
1780 if (strcmp(&infilename[strlen(infilename) - 5],
".mesh") != 0) {
1781 strcat(infilename,
".mesh");
1784 if (!(fp = fopen(infilename,
"r"))) {
1785 printf(
"Error: Unable to open file %s\n", infilename);
1788 printf(
"Opening %s.\n", infilename);
1793 while ((bufferp = readline(buffer, fp, &line_count)) != NULL) {
1794 if (*bufferp ==
'#')
continue;
1795 if (dimension == 0) {
1797 str = strstr(bufferp,
"Dimension");
1798 if (!str) str = strstr(bufferp,
"dimension");
1799 if (!str) str = strstr(bufferp,
"DIMENSION");
1802 bufferp = findnextnumber(str);
1803 if (*bufferp ==
'\0') {
1805 bufferp = readline(buffer, fp, &line_count);
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);
1814 mesh_dim = dimension;
1819 str = strstr(bufferp,
"Vertices");
1820 if (!str) str = strstr(bufferp,
"vertices");
1821 if (!str) str = strstr(bufferp,
"VERTICES");
1824 bufferp = findnextnumber(str);
1825 if (*bufferp ==
'\0') {
1827 bufferp = readline(buffer, fp, &line_count);
1829 nverts = (int) strtol(bufferp, &bufferp, 0);
1832 numberofpoints = nverts;
1833 pointlist =
new REAL[nverts * 3];
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);
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);
1853 if ((j < 2) || (dimension == 3)) {
1854 coord[j] = (REAL) strtod(bufferp, &bufferp);
1856 assert((j == 2) && (dimension == 2));
1859 bufferp = findnextnumber(bufferp);
1868 str = strstr(bufferp,
"Triangles");
1869 if (!str) str = strstr(bufferp,
"triangles");
1870 if (!str) str = strstr(bufferp,
"TRIANGLES");
1874 str = strstr(bufferp,
"Quadrilaterals");
1875 if (!str) str = strstr(bufferp,
"quadrilaterals");
1876 if (!str) str = strstr(bufferp,
"QUADRILATERALS");
1881 if (corners == 3 || corners == 4) {
1883 bufferp = findnextnumber(str);
1884 if (*bufferp ==
'\0') {
1886 bufferp = readline(buffer, fp, &line_count);
1888 nfaces = strtol(bufferp, &bufferp, 0);
1891 if (numberoffacets > 0) {
1893 tmpflist =
new tetgenio::facet[numberoffacets + nfaces];
1894 tmpfmlist =
new int[numberoffacets + nfaces];
1896 for (i = 0; i < numberoffacets; i++) {
1900 tmpfmlist[i] = facetmarkerlist[i];
1903 delete [] facetlist;
1904 delete [] facetmarkerlist;
1906 facetlist = tmpflist;
1907 facetmarkerlist = tmpfmlist;
1910 facetlist =
new tetgenio::facet[nfaces];
1911 facetmarkerlist =
new int[nfaces];
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);
1926 f->numberofpolygons = 1;
1927 f->polygonlist =
new tetgenio::polygon[1];
1928 p = &f->polygonlist[0];
1930 p->numberofvertices = corners;
1932 p->vertexlist =
new int[p->numberofvertices];
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);
1941 p->vertexlist[j] = (int) strtol(bufferp, &bufferp, 0);
1942 if (firstnumber == 1) {
1944 if (p->vertexlist[j] == 0) {
1949 bufferp = findnextnumber(bufferp);
1952 facetmarkerlist[i] = 0;
1953 if (*bufferp !=
'\0') {
1954 facetmarkerlist[i] = (int) strtol(bufferp, &bufferp, 0);
1958 numberoffacets += nfaces;
1981 bool tetgenio::load_plc(
char* filename,
int object)
1983 enum tetgenbehavior::objecttype type;
1985 type = (
enum tetgenbehavior::objecttype)
object;
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);
2000 return load_poly(filename);
2014 bool tetgenio::load_tetmesh(
char* filename)
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;
2026 int markers, corner;
2027 int index, attribindex;
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");
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);
2053 stringptr = readnumberline(inputline, infile, infilename);
2055 stringptr = strstr(inputline,
"rbox");
2056 if (stringptr == NULL) {
2059 stringptr = inputline;
2060 numberofpoints = (int) strtol (stringptr, &stringptr, 0);
2061 stringptr = findnextnumber(stringptr);
2062 if (*stringptr ==
'\0') {
2065 mesh_dim = (int) strtol (stringptr, &stringptr, 0);
2067 stringptr = findnextnumber(stringptr);
2068 if (*stringptr ==
'\0') {
2069 numberofpointattributes = 0;
2071 numberofpointattributes = (int) strtol (stringptr, &stringptr, 0);
2073 stringptr = findnextnumber(stringptr);
2074 if (*stringptr ==
'\0') {
2077 markers = (int) strtol (stringptr, &stringptr, 0);
2081 stringptr = inputline;
2083 mesh_dim = (int) strtol (stringptr, &stringptr, 0);
2085 stringptr = readnumberline(inputline, infile, infilename);
2086 numberofpoints = (int) strtol (stringptr, &stringptr, 0);
2092 if (!load_node_call(infile, markers, infilename)) {
2099 if (mesh_dim == 3) {
2100 infilename = inelefilename;
2101 infile = fopen(infilename,
"r");
2102 if (infile != (FILE *) NULL) {
2103 printf(
"Opening %s.\n", infilename);
2106 stringptr = readnumberline(inputline, infile, infilename);
2107 numberoftetrahedra = (int) strtol (stringptr, &stringptr, 0);
2108 stringptr = findnextnumber(stringptr);
2109 if (*stringptr ==
'\0') {
2110 numberofcorners = 4;
2112 numberofcorners = (int) strtol(stringptr, &stringptr, 0);
2114 stringptr = findnextnumber(stringptr);
2115 if (*stringptr ==
'\0') {
2116 numberoftetrahedronattributes = 0;
2118 numberoftetrahedronattributes = (int) strtol(stringptr, &stringptr, 0);
2120 if (numberofcorners != 4 && numberofcorners != 10) {
2121 printf(
"Error: Wrong number of corners %d (should be 4 or 10).\n",
2127 if (numberoftetrahedra > 0) {
2128 tetrahedronlist =
new int[numberoftetrahedra * numberofcorners];
2129 if (tetrahedronlist == (
int *) NULL) {
2130 printf(
"Error: Out of memory.\n");
2134 if (numberoftetrahedronattributes > 0) {
2135 tetrahedronattributelist =
new REAL[numberoftetrahedra *
2136 numberoftetrahedronattributes];
2137 if (tetrahedronattributelist == (REAL *) NULL) {
2138 printf(
"Error: Out of memory.\n");
2146 for (i = 0; i < numberoftetrahedra; i++) {
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);
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",
2162 tetrahedronlist[index++] = corner;
2165 for (j = 0; j < numberoftetrahedronattributes; j++) {
2166 stringptr = findnextnumber(stringptr);
2167 if (*stringptr ==
'\0') {
2170 attrib = (REAL) strtod(stringptr, &stringptr);
2172 tetrahedronattributelist[attribindex++] = attrib;
2180 if (mesh_dim == 3) {
2181 infilename = infacefilename;
2183 infilename = inelefilename;
2185 infile = fopen(infilename,
"r");
2186 if (infile != (FILE *) NULL) {
2187 printf(
"Opening %s.\n", infilename);
2189 stringptr = readnumberline(inputline, infile, infilename);
2190 numberoftrifaces = (int) strtol (stringptr, &stringptr, 0);
2191 stringptr = findnextnumber(stringptr);
2192 if (mesh_dim == 2) {
2194 stringptr = findnextnumber(stringptr);
2196 if (*stringptr ==
'\0') {
2199 markers = (int) strtol (stringptr, &stringptr, 0);
2201 if (numberoftrifaces > 0) {
2202 trifacelist =
new int[numberoftrifaces * 3];
2203 if (trifacelist == (
int *) NULL) {
2204 printf(
"Error: Out of memory.\n");
2208 trifacemarkerlist =
new int[numberoftrifaces * 3];
2209 if (trifacemarkerlist == (
int *) NULL) {
2210 printf(
"Error: Out of memory.\n");
2217 for (i = 0; i < numberoftrifaces; i++) {
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);
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",
2233 trifacelist[index++] = corner;
2237 stringptr = findnextnumber(stringptr);
2238 if (*stringptr ==
'\0') {
2241 attrib = (REAL) strtod(stringptr, &stringptr);
2243 trifacemarkerlist[i] = (int) attrib;
2250 infilename = inedgefilename;
2251 infile = fopen(infilename,
"r");
2252 if (infile != (FILE *) NULL) {
2253 printf(
"Opening %s.\n", infilename);
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");
2266 for (i = 0; i < numberofedges; i++) {
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);
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",
2282 edgelist[index++] = corner;
2289 infilename = involfilename;
2290 infile = fopen(infilename,
"r");
2291 if (infile != (FILE *) NULL) {
2292 printf(
"Opening %s.\n", infilename);
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);
2301 if (volelements > 0) {
2302 tetrahedronvolumelist =
new REAL[volelements];
2303 if (tetrahedronvolumelist == (REAL *) NULL) {
2304 printf(
"Error: Out of memory.\n");
2309 for (i = 0; i < volelements; i++) {
2310 stringptr = readnumberline(inputline, infile, infilename);
2311 stringptr = findnextnumber(stringptr);
2312 if (*stringptr ==
'\0') {
2315 volume = (REAL) strtod(stringptr, &stringptr);
2317 tetrahedronvolumelist[i] = volume;
2339 bool tetgenio::load_voronoi(
char* filename)
2342 char innodefilename[FILENAMESIZE];
2343 char inedgefilename[FILENAMESIZE];
2344 char inputline[INPUTLINESIZE];
2345 char *stringptr, *infilename;
2348 int firstnode, corner;
2353 strcpy(innodefilename, filename);
2354 strcpy(inedgefilename, filename);
2355 strcat(innodefilename,
".v.node");
2356 strcat(inedgefilename,
".v.edge");
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);
2367 stringptr = readnumberline(inputline, infile, infilename);
2369 stringptr = strstr(inputline,
"rbox");
2370 if (stringptr == NULL) {
2373 stringptr = inputline;
2374 numberofvpoints = (int) strtol (stringptr, &stringptr, 0);
2375 stringptr = findnextnumber(stringptr);
2376 if (*stringptr ==
'\0') {
2379 mesh_dim = (int) strtol (stringptr, &stringptr, 0);
2384 stringptr = inputline;
2386 mesh_dim = (int) strtol (stringptr, &stringptr, 0);
2388 stringptr = readnumberline(inputline, infile, infilename);
2389 numberofvpoints = (int) strtol (stringptr, &stringptr, 0);
2393 vpointlist =
new REAL[numberofvpoints * 3];
2394 if (vpointlist == (REAL *) NULL) {
2395 printf(
"Error: Out of memory.\n");
2400 for (i = 0; i < numberofvpoints; i++) {
2401 stringptr = readnumberline(inputline, infile, infilename);
2404 firstnode = (int) strtol (stringptr, &stringptr, 0);
2405 if ((firstnode == 0) || (firstnode == 1)) {
2406 firstnumber = firstnode;
2409 stringptr = findnextnumber(stringptr);
2411 if (*stringptr ==
'\0') {
2412 printf(
"Error: Point %d has no x coordinate.\n", firstnumber + i);
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);
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);
2428 z = (REAL) strtod(stringptr, &stringptr);
2432 vpointlist[index++] = x;
2433 vpointlist[index++] = y;
2434 vpointlist[index++] = z;
2439 infilename = inedgefilename;
2440 infile = fopen(infilename,
"r");
2441 if (infile != (FILE *) NULL) {
2442 printf(
"Opening %s.\n", infilename);
2444 stringptr = readnumberline(inputline, infile, infilename);
2445 numberofvedges = (int) strtol (stringptr, &stringptr, 0);
2446 if (numberofvedges > 0) {
2447 vedgelist =
new voroedge[numberofvedges];
2451 for (i = 0; i < numberofvedges; i++) {
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);
2462 corner = (int) strtol(stringptr, &stringptr, 0);
2463 j == 0 ? vedge->v1 = corner : vedge->v2 = corner;
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);
2473 vedge->vnormal[j] = (REAL) strtod(stringptr, &stringptr);
2475 if (mesh_dim == 2) {
2476 vedge->vnormal[2] = 0.0;
2479 vedge->vnormal[0] = 0.0;
2480 vedge->vnormal[1] = 0.0;
2481 vedge->vnormal[2] = 0.0;
2498 void tetgenio::save_nodes(
char* filename)
2501 char outnodefilename[FILENAMESIZE];
2502 char outmtrfilename[FILENAMESIZE];
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]);
2515 fprintf(fout,
"%d %.16g %.16g %.16g", i + firstnumber,
2516 pointlist[i * 3], pointlist[i * 3 + 1], pointlist[i * 3 + 2]);
2518 for (j = 0; j < numberofpointattributes; j++) {
2519 fprintf(fout,
" %.16g",
2520 pointattributelist[i * numberofpointattributes + j]);
2522 if (pointmarkerlist != NULL) {
2523 fprintf(fout,
" %d", pointmarkerlist[i]);
2525 fprintf(fout,
"\n");
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]);
2539 fprintf(fout,
"\n");
2553 void tetgenio::save_elements(
char* filename)
2556 char outelefilename[FILENAMESIZE];
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]);
2569 for (j = 0; j < numberoftetrahedronattributes; j++) {
2570 fprintf(fout,
" %g",
2571 tetrahedronattributelist[i * numberoftetrahedronattributes + j]);
2573 fprintf(fout,
"\n");
2587 void tetgenio::save_faces(
char* filename)
2590 char outfacefilename[FILENAMESIZE];
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]);
2604 fprintf(fout,
"\n");
2618 void tetgenio::save_edges(
char* filename)
2621 char outedgefilename[FILENAMESIZE];
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]);
2634 fprintf(fout,
"\n");
2648 void tetgenio::save_neighbors(
char* filename)
2651 char outneighborfilename[FILENAMESIZE];
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]);
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]);
2667 fprintf(fout,
"\n");
2683 void tetgenio::save_poly(
char* filename)
2688 char outpolyfilename[FILENAMESIZE];
2691 sprintf(outpolyfilename,
"%s.poly", filename);
2692 printf(
"Saving poly to %s\n", outpolyfilename);
2693 fout = fopen(outpolyfilename,
"w");
2698 fprintf(fout,
"%d %d %d %d\n", 0, mesh_dim, numberofpointattributes,
2699 pointmarkerlist != NULL ? 1 : 0);
2702 if (mesh_dim == 2) {
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]);
2711 fprintf(fout,
"\n");
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);
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 ");
2728 fprintf(fout,
" %d", p->vertexlist[k]);
2730 fprintf(fout,
"\n");
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]);
2741 fprintf(fout,
"%d\n", numberofholes);
2742 for (i = 0; i < numberofholes; i++) {
2744 fprintf(fout,
"%d %.12g %.12g", i + firstnumber, holelist[i * mesh_dim],
2745 holelist[i * mesh_dim + 1]);
2746 if (mesh_dim == 3) {
2748 fprintf(fout,
" %.12g", holelist[i * mesh_dim + 2]);
2750 fprintf(fout,
"\n");
2754 fprintf(fout,
"%d\n", numberofregions);
2755 for (i = 0; i < numberofregions; i++) {
2756 if (mesh_dim == 2) {
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]);
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]);
2787 char* tetgenio::readline(
char *
string, FILE *infile,
int *linenumber)
2793 result = fgets(
string, INPUTLINESIZE - 1, infile);
2794 if (linenumber) (*linenumber)++;
2795 if (result == (
char *) NULL) {
2796 return (
char *) NULL;
2799 while ((*result ==
' ') || (*result ==
'\t')) result++;
2801 }
while (*result ==
'\0');
2814 char* tetgenio::findnextfield(
char *
string)
2820 while ((*result !=
'\0') && (*result !=
' ') && (*result !=
'\t') &&
2821 (*result !=
',') && (*result !=
';')) {
2826 while ((*result ==
' ') || (*result ==
'\t') || (*result ==
',') ||
2842 char* tetgenio::readnumberline(
char *
string, FILE *infile,
char *infilename)
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);
2858 while ((*result !=
'\0') && (*result !=
'#')
2859 && (*result !=
'.') && (*result !=
'+') && (*result !=
'-')
2860 && ((*result <
'0') || (*result >
'9'))) {
2864 }
while ((*result ==
'#') || (*result ==
'\0'));
2878 char* tetgenio::findnextnumber(
char *
string)
2884 while ((*result !=
'\0') && (*result !=
'#') && (*result !=
' ') &&
2885 (*result !=
'\t') && (*result !=
',')) {
2890 while ((*result !=
'\0') && (*result !=
'#')
2891 && (*result !=
'.') && (*result !=
'+') && (*result !=
'-')
2892 && ((*result <
'0') || (*result >
'9'))) {
2896 if (*result ==
'#') {
2906 static REAL PI = 3.14159265358979323846264338327950288419716939937510582;
2918 tetgenbehavior::tetgenbehavior()
2930 maxdihedral = 165.0;
2936 insertaddpoints = 0;
2973 commandline[0] =
'\0';
2974 infilename[0] =
'\0';
2975 outfilename[0] =
'\0';
2976 addinfilename[0] =
'\0';
2977 bgmeshfilename[0] =
'\0';
2986 void tetgenbehavior::versioninfo()
2988 printf(
"Version 1.4.2 (April 16, 2007).\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");
3002 void tetgenbehavior::syntax()
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.");
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");
3047 void tetgenbehavior::usage()
3050 printf(
"A Quality Tetrahedral Mesh Generator and 3D Delaunay ");
3051 printf(
"Triangulator\n");
3054 printf(
"What Can TetGen Do?\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");
3062 printf(
"Command Line Syntax:\n");
3064 printf(
" Below is the command line syntax of TetGen with a list of ");
3066 printf(
" descriptions. Underscores indicate that numbers may optionally\n");
3067 printf(
" follow certain switches. Do not leave any space between a ");
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");
3078 printf(
"Examples of How to Use TetGen:\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");
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");
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");
3095 printf(
"Please send bugs/comments to Hang Si <si@wias-berlin.de>\n");
3122 bool tetgenbehavior::parse_commandline(
int argc,
char **argv)
3129 char workstring[1024];
3135 commandline[0] =
'\0';
3138 strcpy(commandline, argv[0]);
3139 strcat(commandline,
" ");
3145 for (i = startindex; i < argc; i++) {
3147 strcat(commandline, argv[i]);
3148 strcat(commandline,
" ");
3149 if (startindex == 1) {
3151 if (argv[i][0] !=
'-') {
3152 strncpy(infilename, argv[i], 1024 - 1);
3153 infilename[1024 - 1] =
'\0';
3159 for (j = startindex; argv[i][j] !=
'\0'; j++) {
3160 if (argv[i][j] ==
'p') {
3162 }
else if (argv[i][j] ==
'r') {
3164 }
else if (argv[i][j] ==
'R') {
3166 }
else if (argv[i][j] ==
'q') {
3168 if (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
3169 (argv[i][j + 1] ==
'.')) {
3171 while (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
3172 (argv[i][j + 1] ==
'.')) {
3174 workstring[k] = argv[i][j];
3177 workstring[k] =
'\0';
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);
3190 }
else if (argv[i][j] ==
'm') {
3192 if (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
3193 (argv[i][j + 1] ==
'.')) {
3195 while (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
3196 (argv[i][j + 1] ==
'.')) {
3198 workstring[k] = argv[i][j];
3201 workstring[k] =
'\0';
3203 alpha1 = (REAL) strtod(workstring, (
char **) NULL);
3204 }
else if (metric == 2) {
3205 alpha2 = (REAL) strtod(workstring, (
char **) NULL);
3208 }
else if (argv[i][j] ==
'a') {
3209 if (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
3210 (argv[i][j + 1] ==
'.')) {
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] ==
'+')) {
3217 workstring[k] = argv[i][j];
3220 workstring[k] =
'\0';
3221 maxvolume = (REAL) strtod(workstring, (
char **) NULL);
3225 }
else if (argv[i][j] ==
'A') {
3227 }
else if (argv[i][j] ==
'i') {
3228 insertaddpoints = 1;
3229 }
else if (argv[i][j] ==
'd') {
3231 }
else if (argv[i][j] ==
'z') {
3233 }
else if (argv[i][j] ==
'f') {
3235 }
else if (argv[i][j] ==
'e') {
3237 }
else if (argv[i][j] ==
'n') {
3239 }
else if (argv[i][j] ==
'v') {
3241 }
else if (argv[i][j] ==
'g') {
3243 }
else if (argv[i][j] ==
'G') {
3245 }
else if (argv[i][j] ==
'O') {
3247 }
else if (argv[i][j] ==
'M') {
3249 }
else if (argv[i][j] ==
'Y') {
3251 }
else if (argv[i][j] ==
'J') {
3253 }
else if (argv[i][j] ==
'B') {
3255 }
else if (argv[i][j] ==
'N') {
3257 }
else if (argv[i][j] ==
'E') {
3259 }
else if (argv[i][j] ==
'F') {
3261 }
else if (argv[i][j] ==
'I') {
3263 }
else if (argv[i][j] ==
'o') {
3264 if (argv[i][j + 1] ==
'2') {
3268 }
else if (argv[i][j] ==
'S') {
3269 if (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
3270 (argv[i][j + 1] ==
'.')) {
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] ==
'+')) {
3276 workstring[k] = argv[i][j];
3279 workstring[k] =
'\0';
3280 steiner = (int) strtol(workstring, (
char **) NULL, 0);
3282 }
else if (argv[i][j] ==
's') {
3284 if (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
3285 (argv[i][j + 1] ==
'.')) {
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] ==
'+')) {
3291 workstring[k] = argv[i][j];
3294 workstring[k] =
'\0';
3296 optlevel = (int) strtol(workstring, (
char **) NULL, 0);
3297 }
else if (scount == 2) {
3298 optpasses = (int) strtol(workstring, (
char **) NULL, 0);
3301 }
else if (argv[i][j] ==
'D') {
3303 }
else if (argv[i][j] ==
'T') {
3304 if (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
3305 (argv[i][j + 1] ==
'.')) {
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] ==
'+')) {
3311 workstring[k] = argv[i][j];
3314 workstring[k] =
'\0';
3315 epsilon = (REAL) strtod(workstring, (
char **) NULL);
3317 }
else if (argv[i][j] ==
'C') {
3319 }
else if (argv[i][j] ==
'X') {
3321 }
else if (argv[i][j] ==
'Q') {
3323 }
else if (argv[i][j] ==
'V') {
3328 }
else if ((argv[i][j] ==
'h') || (argv[i][j] ==
'H') ||
3329 (argv[i][j] ==
'?')) {
3333 printf(
"Warning: Unknown switch -%c.\n", argv[i][j]);
3338 if (startindex == 0) {
3340 strcpy(infilename,
"tetgen-tmpfile");
3342 if (infilename[0] ==
'\0') {
3348 if (!strcmp(&infilename[strlen(infilename) - 5],
".node")) {
3349 infilename[strlen(infilename) - 5] =
'\0';
3351 }
else if (!strcmp(&infilename[strlen(infilename) - 5],
".poly")) {
3352 infilename[strlen(infilename) - 5] =
'\0';
3355 }
else if (!strcmp(&infilename[strlen(infilename) - 6],
".smesh")) {
3356 infilename[strlen(infilename) - 6] =
'\0';
3359 }
else if (!strcmp(&infilename[strlen(infilename) - 4],
".off")) {
3360 infilename[strlen(infilename) - 4] =
'\0';
3363 }
else if (!strcmp(&infilename[strlen(infilename) - 4],
".ply")) {
3364 infilename[strlen(infilename) - 4] =
'\0';
3367 }
else if (!strcmp(&infilename[strlen(infilename) - 4],
".stl")) {
3368 infilename[strlen(infilename) - 4] =
'\0';
3371 }
else if (!strcmp(&infilename[strlen(infilename) - 5],
".mesh")) {
3372 infilename[strlen(infilename) - 5] =
'\0';
3375 }
else if (!strcmp(&infilename[strlen(infilename) - 4],
".ele")) {
3376 infilename[strlen(infilename) - 4] =
'\0';
3381 plc = plc || diagnose;
3382 useshelles = plc || refine || coarse || quality;
3383 goodratio = minratio;
3384 goodratio *= goodratio;
3387 if (plc && refine) {
3388 printf(
"Error: Switch -r cannot use together with -p.\n");
3391 if (refine && (plc || noiterationnum)) {
3392 printf(
"Error: Switches %s cannot use together with -r.\n",
3396 if (diagnose && (quality || insertaddpoints || (order == 2) || neighout
3398 printf(
"Error: Switches %s cannot use together with -d.\n",
3399 "-q, -i, -o2, -n, and -C");
3405 if (!refine && !plc) {
3410 if (refine || !plc) {
3414 if (fixedvolume || varvolume) {
3420 goodangle = cos(minangle * PI / 180.0);
3421 goodangle *= goodangle;
3424 strcpy(workstring, infilename);
3426 while (workstring[j] !=
'\0') {
3427 if ((workstring[j] ==
'.') && (workstring[j + 1] !=
'\0')) {
3433 if (increment > 0) {
3436 if ((workstring[j] >=
'0') && (workstring[j] <=
'9')) {
3437 meshnumber = meshnumber * 10 + (int) (workstring[j] -
'0');
3442 }
while (workstring[j] !=
'\0');
3444 if (noiterationnum) {
3445 strcpy(outfilename, infilename);
3446 }
else if (increment == 0) {
3447 strcpy(outfilename, infilename);
3448 strcat(outfilename,
".1");
3450 workstring[increment] =
'%';
3451 workstring[increment + 1] =
'd';
3452 workstring[increment + 2] =
'\0';
3453 sprintf(outfilename, workstring, meshnumber + 1);
3456 strcpy(addinfilename, infilename);
3457 strcat(addinfilename,
".a");
3459 strcpy(bgmeshfilename, infilename);
3460 strcat(bgmeshfilename,
".b");
3483 int tetgenmesh::compare_2_ints(
const void* x,
const void* y) {
3484 if (* (
int *) x < * (
int *) y) {
3486 }
else if (* (
int *) x > * (
int *) y) {
3495 int tetgenmesh::compare_2_longs(
const void* x,
const void* y) {
3496 if (* (
long *) x < * (
long *) y) {
3498 }
else if (* (
long *) x > * (
long *) y) {
3506 int tetgenmesh::compare_2_unsignedlongs(
const void* x,
const void* y) {
3507 if (* (
unsigned long *) x < * (
unsigned long *) y) {
3509 }
else if (* (
unsigned long *) x > * (
unsigned long *) y) {
3531 void tetgenmesh::set_compfunc(
char* str,
int* itbytes, compfunc* pcomp)
3534 if (str[strlen(str) - 1] ==
'*') {
3535 *itbytes =
sizeof(
unsigned long);
3536 *pcomp = &compare_2_unsignedlongs;
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;
3551 printf(
"Error in set_compfunc(): unknown type %s.\n", str);
3566 void tetgenmesh::list::
3567 listinit(
int itbytes, compfunc pcomp,
int mitems,
int exsize)
3570 assert(itbytes > 0 && mitems > 0 && exsize > 0);
3572 itembytes = itbytes;
3575 expandsize = exsize;
3576 base = (
char *) malloc(maxitems * itembytes);
3577 if (base == (
char *) NULL) {
3578 printf(
"Error: Out of memory.\n");
3595 void* tetgenmesh::list::append(
void *appitem)
3598 if (items == maxitems) {
3599 char* newbase = (
char *) realloc(base, (maxitems + expandsize) *
3601 if (newbase == (
char *) NULL) {
3602 printf(
"Error: Out of memory.\n");
3606 maxitems += expandsize;
3608 if (appitem != (
void *) NULL) {
3609 memcpy(base + items * itembytes, appitem, itembytes);
3612 return (
void *) (base + (items - 1) * itembytes);
3626 void* tetgenmesh::list::insert(
int pos,
void* insitem)
3629 return append(insitem);
3632 if (items == maxitems) {
3633 char* newbase = (
char *) realloc(base, (maxitems + expandsize) *
3635 if (newbase == (
char *) NULL) {
3636 printf(
"Error: Out of memory.\n");
3640 maxitems += expandsize;
3643 memmove(base + (pos + 1) * itembytes,
3644 base + pos * itembytes,
3645 (items - pos) * itembytes);
3647 if (insitem != (
void *) NULL) {
3648 memcpy(base + pos * itembytes, insitem, itembytes);
3651 return (
void *) (base + pos * itembytes);
3665 void tetgenmesh::list::del(
int pos,
int order)
3668 if (pos >= 0 && pos < items - 1) {
3671 memmove(base + pos * itembytes,
3672 base + (pos + 1) * itembytes,
3673 (items - pos - 1) * itembytes);
3676 memcpy(base + pos * itembytes,
3677 base + (items - 1) * itembytes,
3696 int tetgenmesh::list::hasitem(
void* checkitem)
3700 for (i = 0; i < items; i++) {
3701 if (comp != (compfunc) NULL) {
3702 if ((* comp)((
void *)(base + i * itembytes), checkitem) == 0) {
3718 void tetgenmesh::list::sort()
3720 qsort((
void *) base, (
size_t) items, (
size_t) itembytes, comp);
3729 tetgenmesh::memorypool::memorypool()
3731 firstblock = nowblock = (
void **) NULL;
3732 nextitem = (
void *) NULL;
3733 deaditemstack = (
void *) NULL;
3734 pathblock = (
void **) NULL;
3735 pathitem = (
void *) NULL;
3736 itemwordtype = POINTER;
3738 itembytes = itemwords = 0;
3740 items = maxitems = 0l;
3741 unallocateditems = 0;
3745 tetgenmesh::memorypool::
3746 memorypool(
int bytecount,
int itemcount,
enum wordtype wtype,
int alignment)
3748 poolinit(bytecount, itemcount, wtype, alignment);
3757 tetgenmesh::memorypool::~memorypool()
3759 while (firstblock != (
void **) NULL) {
3760 nowblock = (
void **) *(firstblock);
3762 firstblock = nowblock;
3782 void tetgenmesh::memorypool::
3783 poolinit(
int bytecount,
int itemcount,
enum wordtype wtype,
int alignment)
3788 itemwordtype = wtype;
3789 wordsize = (itemwordtype == POINTER) ?
sizeof(
void *) :
sizeof(REAL);
3795 if (alignment > wordsize) {
3796 alignbytes = alignment;
3798 alignbytes = wordsize;
3800 if ((
int)
sizeof(
void *) > alignbytes) {
3801 alignbytes = (int)
sizeof(
void *);
3803 itemwords = ((bytecount + alignbytes - 1) / alignbytes)
3804 * (alignbytes / wordsize);
3805 itembytes = itemwords * wordsize;
3806 itemsperblock = itemcount;
3811 firstblock = (
void **) malloc(itemsperblock * itembytes +
sizeof(
void *)
3813 if (firstblock == (
void **) NULL) {
3814 printf(
"Error: Out of memory.\n");
3818 *(firstblock) = (
void *) NULL;
3832 void tetgenmesh::memorypool::restart()
3834 unsigned long alignptr;
3840 nowblock = firstblock;
3842 alignptr = (
unsigned long) (nowblock + 1);
3845 (alignptr + (
unsigned long) alignbytes -
3846 (alignptr % (
unsigned long) alignbytes));
3848 unallocateditems = itemsperblock;
3850 deaditemstack = (
void *) NULL;
3859 void* tetgenmesh::memorypool::alloc()
3863 unsigned long alignptr;
3867 if (deaditemstack != (
void *) NULL) {
3868 newitem = deaditemstack;
3869 deaditemstack = * (
void **) deaditemstack;
3872 if (unallocateditems == 0) {
3874 if (*nowblock == (
void *) NULL) {
3876 newblock = (
void **) malloc(itemsperblock * itembytes +
sizeof(
void *)
3878 if (newblock == (
void **) NULL) {
3879 printf(
"Error: Out of memory.\n");
3882 *nowblock = (
void *) newblock;
3884 *newblock = (
void *) NULL;
3887 nowblock = (
void **) *nowblock;
3890 alignptr = (
unsigned long) (nowblock + 1);
3893 (alignptr + (
unsigned long) alignbytes -
3894 (alignptr % (
unsigned long) alignbytes));
3896 unallocateditems = itemsperblock;
3901 if (itemwordtype == POINTER) {
3902 nextitem = (
void *) ((
void **) nextitem + itemwords);
3904 nextitem = (
void *) ((REAL *) nextitem + itemwords);
3921 void tetgenmesh::memorypool::dealloc(
void *dyingitem)
3924 *((
void **) dyingitem) = deaditemstack;
3925 deaditemstack = dyingitem;
3937 void tetgenmesh::memorypool::traversalinit()
3939 unsigned long alignptr;
3942 pathblock = firstblock;
3944 alignptr = (
unsigned long) (pathblock + 1);
3947 (alignptr + (
unsigned long) alignbytes -
3948 (alignptr % (
unsigned long) alignbytes));
3950 pathitemsleft = itemsperblock;
3965 void* tetgenmesh::memorypool::traverse()
3968 unsigned long alignptr;
3971 if (pathitem == nextitem) {
3972 return (
void *) NULL;
3975 if (pathitemsleft == 0) {
3977 pathblock = (
void **) *pathblock;
3979 alignptr = (
unsigned long) (pathblock + 1);
3982 (alignptr + (
unsigned long) alignbytes -
3983 (alignptr % (
unsigned long) alignbytes));
3985 pathitemsleft = itemsperblock;
3989 if (itemwordtype == POINTER) {
3990 pathitem = (
void *) ((
void **) pathitem + itemwords);
3992 pathitem = (
void *) ((REAL *) pathitem + itemwords);
4007 void tetgenmesh::link::linkinit(
int bytecount, compfunc pcomp,
int itemcount)
4010 assert(bytecount > 0 && itemcount > 0);
4013 linkitembytes = bytecount;
4021 poolinit(bytecount + 2 *
sizeof(
void **), itemcount + 2, POINTER, 0);
4024 head = (
void **) alloc();
4025 tail = (
void **) alloc();
4026 *head = (
void *) tail;
4029 *(tail + 1) = (
void *) head;
4030 nextlinkitem = *head;
4045 void tetgenmesh::link::clear()
4051 head = (
void **) alloc();
4052 tail = (
void **) alloc();
4053 *head = (
void *) tail;
4056 *(tail + 1) = (
void *) head;
4057 nextlinkitem = *head;
4073 bool tetgenmesh::link::move(
int numberofnodes)
4078 nownode = (
void **) nextlinkitem;
4079 if (numberofnodes > 0) {
4082 while ((i < numberofnodes) && *nownode) {
4083 nownode = (
void **) *nownode;
4086 if (*nownode == NULL)
return false;
4087 nextlinkitem = (
void *) nownode;
4088 curpos += numberofnodes;
4089 }
else if (numberofnodes < 0) {
4092 numberofnodes = -numberofnodes;
4093 while ((i < numberofnodes) && *(nownode + 1)) {
4094 nownode = (
void **) *(nownode + 1);
4097 if (*(nownode + 1) == NULL)
return false;
4098 nextlinkitem = (
void *) nownode;
4099 curpos -= numberofnodes;
4116 bool tetgenmesh::link::locate(
int pos)
4118 int headdist, taildist, curdist;
4119 int abscurdist, mindist;
4121 if (pos < 1 || pos > linkitems)
return false;
4124 taildist = linkitems - pos;
4125 curdist = pos - curpos;
4126 abscurdist = curdist >= 0 ? curdist : -curdist;
4128 if (headdist > taildist) {
4129 if (taildist > abscurdist) {
4133 mindist = -taildist;
4138 if (headdist > abscurdist) {
4147 return move(mindist);
4160 void* tetgenmesh::link::add(
void* newitem)
4162 void **newnode = tail;
4163 if (newitem != (
void *) NULL) {
4164 memcpy((
void *)(newnode + 2), newitem, linkitembytes);
4166 tail = (
void **) alloc();
4168 *newnode = (
void*) tail;
4169 *(tail + 1) = (
void*) newnode;
4171 return (
void *)(newnode + 2);
4186 void* tetgenmesh::link::insert(
int pos,
void* insitem)
4189 return add(insitem);
4192 void **nownode = (
void **) nextlinkitem;
4195 void **newnode = (
void **) alloc();
4196 if (insitem != (
void *) NULL) {
4197 memcpy((
void *)(newnode + 2), insitem, linkitembytes);
4200 *(
void **)(*(nownode + 1)) = (
void *) newnode;
4201 *newnode = (
void *) nownode;
4202 *(newnode + 1) = *(nownode + 1);
4203 *(nownode + 1) = (
void *) newnode;
4207 nextlinkitem = (
void *) newnode;
4208 return (
void *)(newnode + 2);
4220 void* tetgenmesh::link::deletenode(
void** deadnode)
4222 void **nextnode = (
void **) *deadnode;
4223 void **prevnode = (
void **) *(deadnode + 1);
4224 *prevnode = (
void *) nextnode;
4225 *(nextnode + 1) = (
void *) prevnode;
4227 dealloc((
void *) deadnode);
4230 nextlinkitem = (
void *) nextnode;
4231 return (
void *)(deadnode + 2);
4244 void* tetgenmesh::link::del(
int pos)
4246 if (!locate(pos) || (linkitems == 0)) {
4247 return (
void *) NULL;
4249 return deletenode((
void **) nextlinkitem);
4262 void* tetgenmesh::link::getitem()
4264 if (nextlinkitem == (
void *) tail)
return NULL;
4265 void **nownode = (
void **) nextlinkitem;
4266 nextlinkitem = *nownode;
4268 return (
void *)(nownode + 2);
4280 void* tetgenmesh::link::getnitem(
int pos)
4282 if (!locate(pos))
return NULL;
4283 return (
void *)((
void **) nextlinkitem + 2);
4296 int tetgenmesh::link::hasitem(
void* checkitem)
4302 pathitem = getitem();
4307 if ((* comp)(pathitem, checkitem) == 0) {
4311 pathitem = getitem();
4330 int tetgenmesh::plus1mod3[3] = {1, 2, 0};
4331 int tetgenmesh::minus1mod3[3] = {2, 0, 1};
4336 int tetgenmesh::ve[6] = { 2, 5, 4, 1, 0, 3 };
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 };
4350 int tetgenmesh::locver2org[4][6] = {
4356 int tetgenmesh::locver2dest[4][6] = {
4362 int tetgenmesh::locver2apex[4][6] = {
4371 int tetgenmesh::loc2oppo[4] = { 3, 2, 0, 1 };
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} }
4388 int tetgenmesh::locver2edge[4][6] = {
4395 int tetgenmesh::edge2locver[6][2] = {
4417 #define Orient(V) ((V) Div2)
4421 #define EdgeRing(V) ((V) Mod2)
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);
4441 inline tetgenmesh::tetrahedron tetgenmesh::encode(triface& t) {
4442 return (tetrahedron) ((
unsigned long) t.tet | (
unsigned long) t.loc);
4447 inline void tetgenmesh::sym(triface& t1, triface& t2) {
4448 tetrahedron ptr = t1.tet[t1.loc];
4452 inline void tetgenmesh::symself(triface& t) {
4453 tetrahedron ptr = t.tet[t.loc];
4459 inline void tetgenmesh::bond(triface& t1, triface& t2) {
4460 t1.tet[t1.loc] = encode(t2);
4461 t2.tet[t2.loc] = encode(t1);
4469 inline void tetgenmesh::dissolve(triface& t) {
4470 t.tet[t.loc] = (tetrahedron) dummytet;
4476 inline tetgenmesh::point tetgenmesh::org(triface& t) {
4477 return (point) t.tet[locver2org[t.loc][t.ver] + 4];
4480 inline tetgenmesh::point tetgenmesh::dest(triface& t) {
4481 return (point) t.tet[locver2dest[t.loc][t.ver] + 4];
4484 inline tetgenmesh::point tetgenmesh::apex(triface& t) {
4485 return (point) t.tet[locver2apex[t.loc][t.ver] + 4];
4488 inline tetgenmesh::point tetgenmesh::oppo(triface& t) {
4489 return (point) t.tet[loc2oppo[t.loc] + 4];
4492 inline void tetgenmesh::setorg(triface& t, point pointptr) {
4493 t.tet[locver2org[t.loc][t.ver] + 4] = (tetrahedron) pointptr;
4496 inline void tetgenmesh::setdest(triface& t, point pointptr) {
4497 t.tet[locver2dest[t.loc][t.ver] + 4] = (tetrahedron) pointptr;
4500 inline void tetgenmesh::setapex(triface& t, point pointptr) {
4501 t.tet[locver2apex[t.loc][t.ver] + 4] = (tetrahedron) pointptr;
4504 inline void tetgenmesh::setoppo(triface& t, point pointptr) {
4505 t.tet[loc2oppo[t.loc] + 4] = (tetrahedron) pointptr;
4515 inline void tetgenmesh::esym(triface& t1, triface& t2) {
4518 t2.ver = t1.ver + (EdgeRing(t1.ver) ? -1 : 1);
4521 inline void tetgenmesh::esymself(triface& t) {
4522 t.ver += (EdgeRing(t.ver) ? -1 : 1);
4527 inline void tetgenmesh::enext(triface& t1, triface& t2) {
4530 t2.ver = ve[t1.ver];
4533 inline void tetgenmesh::enextself(triface& t) {
4539 inline void tetgenmesh::enext2(triface& t1, triface& t2) {
4542 t2.ver = ve[ve[t1.ver]];
4545 inline void tetgenmesh::enext2self(triface& t) {
4546 t.ver = ve[ve[t.ver]];
4553 inline bool tetgenmesh::fnext(triface& t1, triface& t2)
4556 t2.loc = locver2nextf[t1.loc][t1.ver][0];
4560 t2.ver = locver2nextf[t1.loc][t1.ver][1];
4565 if (t2.tet != dummytet) {
4571 for (i = 0; (i < 3) && (org(t2) != torg); i++) {
4577 t2.loc = locver2nextf[tloc][tver][0];
4578 t2.ver = locver2nextf[tloc][tver][1];
4581 return t2.tet != dummytet;
4584 inline bool tetgenmesh::fnextself(triface& t1)
4589 t2.loc = locver2nextf[t1.loc][t1.ver][0];
4593 t2.ver = locver2nextf[t1.loc][t1.ver][1];
4599 if (t2.tet != dummytet) {
4605 for (i = 0; (i < 3) && (org(t2) != torg); i++) {
4608 t1.loc = locver2nextf[t2.loc][t2.ver][0];
4609 t1.ver = locver2nextf[t2.loc][t2.ver][1];
4613 return t2.tet != dummytet;
4619 inline void tetgenmesh::enextfnext(triface& t1, triface& t2) {
4624 inline void tetgenmesh::enextfnextself(triface& t) {
4629 inline void tetgenmesh::enext2fnext(triface& t1, triface& t2) {
4634 inline void tetgenmesh::enext2fnextself(triface& t) {
4643 inline void tetgenmesh::infect(triface& t) {
4644 t.tet[0] = (tetrahedron) ((
unsigned long) t.tet[0] | (
unsigned long) 4l);
4647 inline void tetgenmesh::uninfect(triface& t) {
4648 t.tet[0] = (tetrahedron) ((
unsigned long) t.tet[0] & ~ (
unsigned long) 4l);
4653 inline bool tetgenmesh::infected(triface& t) {
4654 return (((
unsigned long) t.tet[0] & (
unsigned long) 4l) != 0);
4659 inline REAL tetgenmesh::elemattribute(tetrahedron* ptr,
int attnum) {
4660 return ((REAL *) (ptr))[elemattribindex + attnum];
4663 inline void tetgenmesh::
4664 setelemattribute(tetrahedron* ptr,
int attnum, REAL value){
4665 ((REAL *) (ptr))[elemattribindex + attnum] = value;
4670 inline REAL tetgenmesh::volumebound(tetrahedron* ptr) {
4671 return ((REAL *) (ptr))[volumeboundindex];
4674 inline void tetgenmesh::setvolumebound(tetrahedron* ptr, REAL value) {
4675 ((REAL *) (ptr))[volumeboundindex] = value;
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);
4699 inline tetgenmesh::shellface tetgenmesh::sencode(face& s) {
4700 return (shellface) ((
unsigned long) s.sh | (
unsigned long) s.shver);
4706 inline void tetgenmesh::spivot(face& s1, face& s2) {
4707 shellface sptr = s1.sh[Orient(s1.shver)];
4711 inline void tetgenmesh::spivotself(face& s) {
4712 shellface sptr = s.sh[Orient(s.shver)];
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);
4727 inline void tetgenmesh::sbond1(face& s1, face& s2) {
4728 s1.sh[Orient(s1.shver)] = sencode(s2);
4734 inline void tetgenmesh::sdissolve(face& s) {
4735 s.sh[Orient(s.shver)] = (shellface) dummysh;
4741 inline tetgenmesh::point tetgenmesh::sorg(face& s) {
4742 return (point) s.sh[3 + vo[s.shver]];
4745 inline tetgenmesh::point tetgenmesh::sdest(face& s) {
4746 return (point) s.sh[3 + vd[s.shver]];
4749 inline tetgenmesh::point tetgenmesh::sapex(face& s) {
4750 return (point) s.sh[3 + va[s.shver]];
4753 inline void tetgenmesh::setsorg(face& s, point pointptr) {
4754 s.sh[3 + vo[s.shver]] = (shellface) pointptr;
4757 inline void tetgenmesh::setsdest(face& s, point pointptr) {
4758 s.sh[3 + vd[s.shver]] = (shellface) pointptr;
4761 inline void tetgenmesh::setsapex(face& s, point pointptr) {
4762 s.sh[3 + va[s.shver]] = (shellface) pointptr;
4768 inline void tetgenmesh::sesym(face& s1, face& s2) {
4770 s2.shver = s1.shver + (EdgeRing(s1.shver) ? -1 : 1);
4773 inline void tetgenmesh::sesymself(face& s) {
4774 s.shver += (EdgeRing(s.shver) ? -1 : 1);
4777 inline void tetgenmesh::senext(face& s1, face& s2) {
4779 s2.shver = ve[s1.shver];
4782 inline void tetgenmesh::senextself(face& s) {
4783 s.shver = ve[s.shver];
4786 inline void tetgenmesh::senext2(face& s1, face& s2) {
4788 s2.shver = ve[ve[s1.shver]];
4791 inline void tetgenmesh::senext2self(face& s) {
4792 s.shver = ve[ve[s.shver]];
4797 inline void tetgenmesh::sfnext(face& s1, face& s2) {
4798 getnextsface(&s1, &s2);
4801 inline void tetgenmesh::sfnextself(face& s) {
4802 getnextsface(&s, NULL);
4808 inline tetgenmesh::badface* tetgenmesh::shell2badface(face& s) {
4809 return (badface*) s.sh[11];
4812 inline void tetgenmesh::setshell2badface(face& s, badface* value) {
4813 s.sh[11] = (shellface) value;
4818 inline REAL tetgenmesh::areabound(face& s) {
4819 return ((REAL *) (s.sh))[areaboundindex];
4822 inline void tetgenmesh::setareabound(face& s, REAL value) {
4823 ((REAL *) (s.sh))[areaboundindex] = value;
4829 inline int tetgenmesh::shellmark(face& s) {
4830 return ((
int *) (s.sh))[shmarkindex];
4833 inline void tetgenmesh::setshellmark(face& s,
int value) {
4834 ((
int *) (s.sh))[shmarkindex] = value;
4839 inline enum tetgenmesh::shestype tetgenmesh::shelltype(face& s) {
4840 return (
enum shestype) ((
int *) (s.sh))[shmarkindex + 1];
4843 inline void tetgenmesh::setshelltype(face& s,
enum shestype value) {
4844 ((
int *) (s.sh))[shmarkindex + 1] = (
int) value;
4849 inline int tetgenmesh::shellpbcgroup(face& s) {
4850 return ((
int *) (s.sh))[shmarkindex + 2];
4853 inline void tetgenmesh::setshellpbcgroup(face& s,
int value) {
4854 ((
int *) (s.sh))[shmarkindex + 2] = value;
4860 inline void tetgenmesh::sinfect(face& s) {
4861 s.sh[6] = (shellface) ((
unsigned long) s.sh[6] | (
unsigned long) 4l);
4864 inline void tetgenmesh::suninfect(face& s) {
4865 s.sh[6] = (shellface)((
unsigned long) s.sh[6] & ~(
unsigned long) 4l);
4870 inline bool tetgenmesh::sinfected(face& s) {
4871 return (((
unsigned long) s.sh[6] & (
unsigned long) 4l) != 0);
4884 inline void tetgenmesh::tspivot(triface& t, face& s) {
4885 shellface sptr = (shellface) t.tet[8 + t.loc];
4891 inline void tetgenmesh::stpivot(face& s, triface& t) {
4892 tetrahedron ptr = (tetrahedron) s.sh[6 + EdgeRing(s.shver)];
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);
4905 inline void tetgenmesh::tsdissolve(triface& t) {
4906 t.tet[8 + t.loc] = (tetrahedron) dummysh;
4911 inline void tetgenmesh::stdissolve(face& s) {
4912 s.sh[6 + EdgeRing(s.shver)] = (shellface) dummytet;
4925 inline void tetgenmesh::sspivot(face& s, face& edge) {
4926 shellface sptr = (shellface) s.sh[8 + Orient(s.shver)];
4927 sdecode(sptr, edge);
4932 inline void tetgenmesh::ssbond(face& s, face& edge) {
4933 s.sh[8 + Orient(s.shver)] = sencode(edge);
4934 edge.sh[0] = sencode(s);
4939 inline void tetgenmesh::ssdissolve(face& s) {
4940 s.sh[8 + Orient(s.shver)] = (shellface) dummysh;
4951 inline void tetgenmesh::tsspivot1(triface& t, face& seg)
4953 shellface sptr = (shellface) t.tet[8 + locver2edge[t.loc][t.ver]];
4959 inline void tetgenmesh::tssbond1(triface& t, face& seg)
4961 t.tet[8 + locver2edge[t.loc][t.ver]] = (tetrahedron) sencode(seg);
4964 inline void tetgenmesh::tssdissolve1(triface& t)
4966 t.tet[8 + locver2edge[t.loc][t.ver]] = (tetrahedron) dummysh;
4977 inline int tetgenmesh::pointmark(point pt) {
4978 return ((
int *) (pt))[pointmarkindex];
4981 inline void tetgenmesh::setpointmark(point pt,
int value) {
4982 ((
int *) (pt))[pointmarkindex] = value;
4987 inline enum tetgenmesh::verttype tetgenmesh::pointtype(point pt) {
4988 return (
enum verttype) ((
int *) (pt))[pointmarkindex + 1];
4991 inline void tetgenmesh::setpointtype(point pt,
enum verttype value) {
4992 ((
int *) (pt))[pointmarkindex + 1] = (
int) value;
4998 inline tetgenmesh::tetrahedron tetgenmesh::point2tet(point pt) {
4999 return ((tetrahedron *) (pt))[point2simindex];
5002 inline void tetgenmesh::setpoint2tet(point pt, tetrahedron value) {
5003 ((tetrahedron *) (pt))[point2simindex] = value;
5006 inline tetgenmesh::shellface tetgenmesh::point2sh(point pt) {
5007 return (shellface) ((tetrahedron *) (pt))[point2simindex + 1];
5010 inline void tetgenmesh::setpoint2sh(point pt, shellface value) {
5011 ((tetrahedron *) (pt))[point2simindex + 1] = (tetrahedron) value;
5014 inline tetgenmesh::point tetgenmesh::point2ppt(point pt) {
5015 return (point) ((tetrahedron *) (pt))[point2simindex + 2];
5018 inline void tetgenmesh::setpoint2ppt(point pt, point value) {
5019 ((tetrahedron *) (pt))[point2simindex + 2] = (tetrahedron) value;
5022 inline tetgenmesh::tetrahedron tetgenmesh::point2bgmtet(point pt) {
5023 return ((tetrahedron *) (pt))[point2simindex + 3];
5026 inline void tetgenmesh::setpoint2bgmtet(point pt, tetrahedron value) {
5027 ((tetrahedron *) (pt))[point2simindex + 3] = value;
5032 inline tetgenmesh::point tetgenmesh::point2pbcpt(point pt) {
5033 return (point) ((tetrahedron *) (pt))[point2pbcptindex];
5036 inline void tetgenmesh::setpoint2pbcpt(point pt, point value) {
5037 ((tetrahedron *) (pt))[point2pbcptindex] = (tetrahedron) value;
5052 inline void tetgenmesh::adjustedgering(triface& t,
int direction) {
5053 if (EdgeRing(t.ver) != direction) {
5058 inline void tetgenmesh::adjustedgering(face& s,
int direction) {
5059 if (EdgeRing(s.shver) != direction) {
5066 inline bool tetgenmesh::isdead(triface* t) {
5067 if (t->tet == (tetrahedron *) NULL)
return true;
5068 else return t->tet[4] == (tetrahedron) NULL;
5071 inline bool tetgenmesh::isdead(face* s) {
5072 if (s->sh == (shellface *) NULL)
return true;
5073 else return s->sh[3] == (shellface) NULL;
5079 inline bool tetgenmesh::isfacehaspoint(triface* t, point testpoint) {
5080 return ((org(*t) == testpoint) || (dest(*t) == testpoint) ||
5081 (apex(*t) == testpoint));
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);
5093 inline bool tetgenmesh::isfacehasedge(face* s, point tend1, point tend2) {
5094 return (isfacehaspoint(s, tend1) && isfacehaspoint(s, tend2));
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;
5115 void tetgenmesh::getnextsface(face* s1, face* s2)
5117 face neighsh, spinsh;
5120 sspivot(*s1, testseg);
5121 if (testseg.sh != dummysh) {
5123 if (sorg(testseg) == sorg(*s1)) {
5124 spivot(*s1, neighsh);
5130 }
while (spinsh.sh != s1->sh);
5133 spivot(*s1, neighsh);
5135 if (sorg(neighsh) != sorg(*s1)) {
5138 if (s2 != (face *) NULL) {
5161 void tetgenmesh::tsspivot(triface* checkedge, face* checkseg)
5168 spintet = *checkedge;
5169 tapex = apex(*checkedge);
5172 tspivot(spintet, parentsh);
5174 if ((parentsh.sh != dummysh) && (sapex(parentsh) != NULL)) {
5176 findedge(&parentsh, org(*checkedge), dest(*checkedge));
5177 sspivot(parentsh, *checkseg);
5178 if (checkseg->sh != dummysh) {
5180 if (sorg(*checkseg) != org(*checkedge)) {
5181 sesymself(*checkseg);
5186 if (!fnextself(spintet)) {
5189 esym(*checkedge, spintet);
5190 if (!fnextself(spintet)) {
5195 }
while ((apex(spintet) != tapex) && (hitbdry < 2));
5197 checkseg->sh = dummysh;
5211 void tetgenmesh::sstpivot(face* checkseg, triface* retedge)
5216 sdecode(checkseg->sh[0], parentsh);
5218 assert(parentsh.sh != dummysh);
5221 stpivot(parentsh, *retedge);
5222 if (retedge->tet == dummytet) {
5223 sesymself(parentsh);
5224 stpivot(parentsh, *retedge);
5226 assert(retedge->tet != dummytet);
5230 findedge(retedge, sorg(*checkseg), sdest(*checkseg));
5243 bool tetgenmesh::findorg(triface* tface, point dorg)
5245 if (org(*tface) == dorg) {
5248 if (dest(*tface) == dorg) {
5252 if (apex(*tface) == dorg) {
5256 if (oppo(*tface) == dorg) {
5258 adjustedgering(*tface, CCW);
5269 bool tetgenmesh::findorg(face* sface, point dorg)
5271 if (sorg(*sface) == dorg) {
5274 if (sdest(*sface) == dorg) {
5278 if (sapex(*sface) == dorg) {
5279 senext2self(*sface);
5297 void tetgenmesh::findedge(triface* tface, point eorg, point edest)
5301 for (i = 0; i < 3; i++) {
5302 if (org(*tface) == eorg) {
5303 if (dest(*tface) == edest) {
5308 if (org(*tface) == edest) {
5309 if (dest(*tface) == eorg) {
5319 printf(
"Internalerror in findedge(): Unable to find an edge in tet.\n");
5323 void tetgenmesh::findedge(face* sface, point eorg, point edest)
5327 for (i = 0; i < 3; i++) {
5328 if (sorg(*sface) == eorg) {
5329 if (sdest(*sface) == edest) {
5334 if (sorg(*sface) == edest) {
5335 if (sdest(*sface) == eorg) {
5344 printf(
"Internalerror in findedge(): Unable to find an edge in subface.\n");
5358 void tetgenmesh::findface(triface *fface, point forg, point fdest, point fapex)
5361 enum finddirectionresult collinear;
5364 if (!isdead(fface)) {
5366 if (org(*fface) == forg && dest(*fface) == fdest &&
5367 apex(*fface) == fapex)
return;
5370 if (!isdead(&recenttet)) *fface = recenttet;
5373 if (!isdead(fface)) {
5374 if (!findorg(fface, forg)) {
5376 preciselocate(forg, fface, tetrahedrons->items);
5379 if (org(*fface) == forg) {
5380 collinear = finddirection(fface, fdest, tetrahedrons->items);
5381 if (collinear == RIGHTCOLLINEAR) {
5383 }
else if (collinear == LEFTCOLLINEAR) {
5386 }
else if (collinear == TOPCOLLINEAR) {
5393 if ((org(*fface) == forg) && (dest(*fface) == fdest)) {
5398 if (apex(spintet) == fapex) {
5401 if (org(spintet) != org(*fface)) {
5407 if (!fnextself(spintet)) {
5410 esym(*fface, spintet);
5411 if (!fnextself(spintet)) {
5416 }
while (hitbdry < 2 && apex(spintet) != apex(*fface));
5421 if (isdead(fface) || (org(*fface) != forg) || (dest(*fface) != fdest) ||
5422 (apex(*fface) != fapex)) {
5426 if (b->verbose > 1) {
5427 printf(
"Warning in findface(): Perform a brute-force searching.\n");
5429 enum verttype forgty, fdestty, fapexty;
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) {
5441 for (i = 0; i < 4; i++) {
5442 if (pointtype((point) fface->tet[4 + i]) == DEADVERTEX) share ++;
5446 if (pointtype((point) fface->tet[4]) != DEADVERTEX) {
5448 }
else if (pointtype((point) fface->tet[5]) != DEADVERTEX) {
5450 }
else if (pointtype((point) fface->tet[6]) != DEADVERTEX) {
5455 findedge(fface, forg, fdest);
5458 fface->tet = tetrahedrontraverse();
5460 setpointtype(forg, forgty);
5461 setpointtype(fdest, fdestty);
5462 setpointtype(fapex, fapexty);
5463 if (fface->tet == (tetrahedron *) NULL) {
5465 printf(
"Internal error: Fail to find the indicated face.\n");
5480 void tetgenmesh::getonextseg(face* s, face* lseg)
5482 face checksh, checkseg;
5489 senext2self(checksh);
5491 sspivot(checksh, checkseg);
5492 if (checkseg.sh != dummysh)
break;
5494 spivotself(checksh);
5497 assert(checksh.sh != s->sh);
5499 if (sorg(checksh) != forg) {
5502 assert(sorg(checksh) == forg);
5506 if (sorg(checkseg) != forg) sesymself(checkseg);
5520 void tetgenmesh::getseghasorg(face* sseg, point dorg)
5526 checkpt = sorg(nextseg);
5527 while ((checkpt != dorg) && (pointtype(checkpt) == FREESEGVERTEX)) {
5529 senext2self(nextseg);
5530 spivotself(nextseg);
5532 if (sdest(nextseg) != checkpt) sesymself(nextseg);
5533 checkpt = sorg(nextseg);
5535 if (checkpt == dorg) {
5540 checkpt = sdest(nextseg);
5541 while ((checkpt != dorg) && (pointtype(checkpt) == FREESEGVERTEX)) {
5543 senextself(nextseg);
5544 spivotself(nextseg);
5546 if (sorg(nextseg) != checkpt) sesymself(nextseg);
5547 checkpt = sdest(nextseg);
5549 if (checkpt == dorg) {
5550 sesym(nextseg, *sseg);
5554 printf(
"Internalerror in getseghasorg(): Unable to find the subseg.\n");
5564 tetgenmesh::point tetgenmesh::getsubsegfarorg(face* sseg)
5569 checkpt = sorg(*sseg);
5570 senext2(*sseg, prevseg);
5571 spivotself(prevseg);
5573 while (prevseg.sh != dummysh) {
5575 if (sdest(prevseg) != checkpt) sesymself(prevseg);
5576 checkpt = sorg(prevseg);
5577 senext2self(prevseg);
5578 spivotself(prevseg);
5589 tetgenmesh::point tetgenmesh::getsubsegfardest(face* sseg)
5594 checkpt = sdest(*sseg);
5595 senext(*sseg, nextseg);
5596 spivotself(nextseg);
5598 while (nextseg.sh != dummysh) {
5600 if (sorg(nextseg) != checkpt) sesymself(nextseg);
5601 checkpt = sdest(nextseg);
5602 senextself(nextseg);
5603 spivotself(nextseg);
5616 void tetgenmesh::printtet(triface* tface)
5618 triface tmpface, prtface;
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)");
5632 while(facecount < 4) {
5633 tmpface.loc = facecount;
5634 sym(tmpface, prtface);
5635 if(prtface.tet == dummytet) {
5636 printf(
" [%i] Outer space.\n", facecount);
5638 printf(
" [%i] x%lx loc(%i).", facecount,
5639 (
unsigned long)(prtface.tet), prtface.loc);
5640 if (infected(prtface)) {
5641 printf(
" (infected)");
5648 tmppt = org(*tface);
5649 if(tmppt == (point) NULL) {
5650 printf(
" Org [%i] NULL\n", locver2org[tface->loc][tface->ver]);
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));
5656 tmppt = dest(*tface);
5657 if(tmppt == (point) NULL) {
5658 printf(
" Dest[%i] NULL\n", locver2dest[tface->loc][tface->ver]);
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));
5664 tmppt = apex(*tface);
5665 if(tmppt == (point) NULL) {
5666 printf(
" Apex[%i] NULL\n", locver2apex[tface->loc][tface->ver]);
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));
5672 tmppt = oppo(*tface);
5673 if(tmppt == (point) NULL) {
5674 printf(
" Oppo[%i] NULL\n", loc2oppo[tface->loc]);
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));
5681 if (b->useshelles) {
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) {
5708 void tetgenmesh::printsh(face* sface)
5714 if (sapex(*sface) != NULL) {
5715 printf(
"subface x%lx, ver %d, mark %d:",
5716 (
unsigned long)(sface->sh), sface->shver, shellmark(*sface));
5718 printf(
"Subsegment x%lx, ver %d, mark %d:",
5719 (
unsigned long)(sface->sh), sface->shver, shellmark(*sface));
5721 if (sinfected(*sface)) {
5722 printf(
" (infected)");
5724 if (shell2badface(*sface)) {
5725 printf(
" (queued)");
5727 if (sapex(*sface) != NULL) {
5728 if (shelltype(*sface) == SHARP) {
5732 if (shelltype(*sface) == SHARP) {
5737 if (shellpbcgroup(*sface) >= 0) {
5738 printf(
" (pbc %d)", shellpbcgroup(*sface));
5743 sdecode(sface->sh[0], prtsh);
5744 if (prtsh.sh == dummysh) {
5745 printf(
" [0] = No shell\n");
5747 printf(
" [0] = x%lx %d\n", (
unsigned long)(prtsh.sh), prtsh.shver);
5749 sdecode(sface->sh[1], prtsh);
5750 if (prtsh.sh == dummysh) {
5751 printf(
" [1] = No shell\n");
5753 printf(
" [1] = x%lx %d\n", (
unsigned long)(prtsh.sh), prtsh.shver);
5755 sdecode(sface->sh[2], prtsh);
5756 if (prtsh.sh == dummysh) {
5757 printf(
" [2] = No shell\n");
5759 printf(
" [2] = x%lx %d\n", (
unsigned long)(prtsh.sh), prtsh.shver);
5762 printpoint = sorg(*sface);
5763 if (printpoint == (point) NULL)
5764 printf(
" Org [%d] = NULL\n", vo[sface->shver]);
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]);
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));
5777 if (sapex(*sface) != NULL) {
5778 printpoint = sapex(*sface);
5779 if (printpoint == (point) NULL)
5780 printf(
" Apex[%d] = NULL\n", va[sface->shver]);
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));
5786 decode(sface->sh[6], prttet);
5787 if (prttet.tet == dummytet) {
5788 printf(
" [6] = Outer space\n");
5790 printf(
" [6] = x%lx %d\n",
5791 (
unsigned long)(prttet.tet), prttet.loc);
5793 decode(sface->sh[7], prttet);
5794 if (prttet.tet == dummytet) {
5795 printf(
" [7] = Outer space\n");
5797 printf(
" [7] = x%lx %d\n",
5798 (
unsigned long)(prttet.tet), prttet.loc);
5801 sdecode(sface->sh[8], prtsh);
5802 if (prtsh.sh == dummysh) {
5803 printf(
" [8] = No subsegment\n");
5805 printf(
" [8] = x%lx %d\n",
5806 (
unsigned long)(prtsh.sh), prtsh.shver);
5808 sdecode(sface->sh[9], prtsh);
5809 if (prtsh.sh == dummysh) {
5810 printf(
" [9] = No subsegment\n");
5812 printf(
" [9] = x%lx %d\n",
5813 (
unsigned long)(prtsh.sh), prtsh.shver);
5815 sdecode(sface->sh[10], prtsh);
5816 if (prtsh.sh == dummysh) {
5817 printf(
" [10]= No subsegment\n");
5819 printf(
" [10]= x%lx %d\n",
5820 (
unsigned long)(prtsh.sh), prtsh.shver);
5848 void tetgenmesh::makepoint2tetmap()
5853 if (b->verbose > 0) {
5854 printf(
" Constructing mapping from points to tetrahedra.\n");
5858 points->traversalinit();
5859 pointptr = pointtraverse();
5860 while (pointptr != (point) NULL) {
5861 setpoint2tet(pointptr, (tetrahedron) NULL);
5862 pointptr = pointtraverse();
5865 tetrahedrons->traversalinit();
5866 tetloop.tet = tetrahedrontraverse();
5867 while (tetloop.tet != (tetrahedron *) NULL) {
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));
5879 tetloop.tet = tetrahedrontraverse();
5894 void tetgenmesh::makeindex2pointmap(point*& idx2verlist)
5899 if (b->verbose > 0) {
5900 printf(
" Constructing mapping from indices to points.\n");
5903 idx2verlist =
new point[points->items];
5905 points->traversalinit();
5906 pointloop = pointtraverse();
5908 while (pointloop != (point) NULL) {
5909 idx2verlist[idx] = pointloop;
5911 pointloop = pointtraverse();
5933 void tetgenmesh::makesegmentmap(
int*& idx2seglist, shellface**& segsperverlist)
5938 if (b->verbose > 0) {
5939 printf(
" Constructing mapping from points to segments.\n");
5943 idx2seglist =
new int[points->items + 1];
5944 for (i = 0; i < points->items + 1; i++) idx2seglist[i] = 0;
5948 subsegs->traversalinit();
5949 shloop = shellfacetraverse(subsegs);
5950 while (shloop != (shellface *) NULL) {
5952 for (i = 0; i < 2; i++) {
5953 j = pointmark((point) shloop[3 + i]) - in->firstnumber;
5956 shloop = shellfacetraverse(subsegs);
5962 for (i = 0; i < points->items; i++) {
5963 k = idx2seglist[i + 1];
5964 idx2seglist[i + 1] = idx2seglist[i] + j;
5968 segsperverlist =
new shellface*[idx2seglist[i]];
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;
5978 shloop = shellfacetraverse(subsegs);
5981 for (i = points->items - 1; i >= 0; i--) {
5982 idx2seglist[i + 1] = idx2seglist[i];
6006 makesubfacemap(
int*& idx2facelist, shellface**& facesperverlist)
6011 if (b->verbose > 0) {
6012 printf(
" Constructing mapping from points to subfaces.\n");
6016 idx2facelist =
new int[points->items + 1];
6017 for (i = 0; i < points->items + 1; i++) idx2facelist[i] = 0;
6021 subfaces->traversalinit();
6022 shloop = shellfacetraverse(subfaces);
6023 while (shloop != (shellface *) NULL) {
6025 for (i = 0; i < 3; i++) {
6026 j = pointmark((point) shloop[3 + i]) - in->firstnumber;
6029 shloop = shellfacetraverse(subfaces);
6033 j = idx2facelist[0];
6034 idx2facelist[0] = 0;
6035 for (i = 0; i < points->items; i++) {
6036 k = idx2facelist[i + 1];
6037 idx2facelist[i + 1] = idx2facelist[i] + j;
6041 facesperverlist =
new shellface*[idx2facelist[i]];
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;
6051 shloop = shellfacetraverse(subfaces);
6054 for (i = points->items - 1; i >= 0; i--) {
6055 idx2facelist[i + 1] = idx2facelist[i];
6057 idx2facelist[0] = 0;
6079 maketetrahedronmap(
int*& idx2tetlist, tetrahedron**& tetsperverlist)
6081 tetrahedron *tetloop;
6084 if (b->verbose > 0) {
6085 printf(
" Constructing mapping from points to tetrahedra.\n");
6089 idx2tetlist =
new int[points->items + 1];
6090 for (i = 0; i < points->items + 1; i++) idx2tetlist[i] = 0;
6094 tetrahedrons->traversalinit();
6095 tetloop = tetrahedrontraverse();
6096 while (tetloop != (tetrahedron *) NULL) {
6098 for (i = 0; i < 4; i++) {
6099 j = pointmark((point) tetloop[4 + i]) - in->firstnumber;
6102 tetloop = tetrahedrontraverse();
6108 for (i = 0; i < points->items; i++) {
6109 k = idx2tetlist[i + 1];
6110 idx2tetlist[i + 1] = idx2tetlist[i] + j;
6114 tetsperverlist =
new tetrahedron*[idx2tetlist[i]];
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;
6124 tetloop = tetrahedrontraverse();
6127 for (i = points->items - 1; i >= 0; i--) {
6128 idx2tetlist[i + 1] = idx2tetlist[i];
6143 inline REAL tetgenmesh::dot(REAL* v1, REAL* v2)
6145 return v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2];
6150 inline void tetgenmesh::cross(REAL* v1, REAL* v2, REAL* n)
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];
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,
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;
6171 static void m4xm4(REAL m1[4][4], REAL m2[4][4])
6176 for (i = 0; i < 4; i++) {
6177 for (j = 0; j < 4; j++) {
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];
6181 for (j = 0; j < 4; j++)
6187 static void m4xv4(REAL v2[4], REAL m[4][4], REAL v1[4])
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];
6219 bool tetgenmesh::lu_decmp(REAL lu[4][4],
int n,
int* ps, REAL* d,
int N)
6222 REAL pivot, biggest, mult, tempf;
6228 for (i = N; i < n + N; i++) {
6231 for (j = N; j < n + N; j++)
6232 if (biggest < (tempf = fabs(lu[i][j])))
6235 scales[i] = 1.0 / biggest;
6243 for (k = N; k < n + N - 1; k++) {
6246 for (i = k; i < n + N; i++) {
6247 if (biggest < (tempf = fabs(lu[ps[i]][k]) * scales[ps[i]])) {
6252 if (biggest == 0.0) {
6255 if (pivotindex != k) {
6257 ps[k] = ps[pivotindex];
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;
6267 for (j = k + 1; j < n + N; j++)
6268 lu[ps[i]][j] -= mult * lu[ps[k]][j];
6274 return lu[ps[n + N - 1]][n + N - 1] != 0.0;
6292 void tetgenmesh::lu_solve(REAL lu[4][4],
int n,
int* ps, REAL* b,
int N)
6297 for (i = N; i < n + N; i++) X[i] = 0.0;
6300 for (i = N; i < n + N; i++) {
6302 for (j = N; j < i + N; j++)
6303 dot += lu[ps[i]][j] * X[j];
6304 X[i] = b[ps[i]] - dot;
6308 for (i = n + N - 1; i >= N; i--) {
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];
6315 for (i = N; i < n + N; i++) b[i] = X[i];
6346 enum tetgenmesh::interresult tetgenmesh::edge_vert_col_inter(REAL* A, REAL* B,
6354 }
else if (P[i] > A[i]) {
6357 }
else if (P[i] > B[i]) {
6367 }
else if (A[i] > B[i]) {
6370 }
else if (P[i] > B[i]) {
6373 }
else if (P[i] > A[i]) {
6409 enum tetgenmesh::interresult tetgenmesh:: edge_edge_cop_inter(REAL* A, REAL* B,
6410 REAL* P, REAL* Q, REAL* R)
6412 REAL s1, s2, s3, s4;
6417 s1 = orient3d(A, B, R, P);
6418 s2 = orient3d(A, B, R, Q);
6419 if (s1 * s2 > 0.0) {
6423 s3 = orient3d(P, Q, R, A);
6424 s4 = orient3d(P, Q, R, B);
6425 if (s3 * s4 > 0.0) {
6434 enum interresult abp, abq;
6435 enum interresult pqa, pqb;
6439 abp = edge_vert_col_inter(A, B, P);
6440 if (abp == INTERSECT) {
6446 abq = edge_vert_col_inter(A, B, Q);
6447 if (abq == INTERSECT) {
6451 if (abp == SHAREVERTEX && abq == SHAREVERTEX) {
6455 pqa = edge_vert_col_inter(P, Q, A);
6456 if (pqa == INTERSECT) {
6460 pqb = edge_vert_col_inter(P, Q, B);
6461 if (pqb == INTERSECT) {
6465 if (abp == SHAREVERTEX || abq == SHAREVERTEX) {
6475 assert((abp == DISJOINT) && (abp == abq && abq == pqa && pqa == pqb));
6481 assert(abp == SHAREVERTEX || abp == DISJOINT);
6489 abq = edge_vert_col_inter(A, B, Q);
6491 assert(abq == SHAREVERTEX || abq == DISJOINT || abq == INTERSECT);
6503 pqa = edge_vert_col_inter(P, Q, A);
6506 assert(pqa != SHAREVERTEX);
6507 assert(pqa == INTERSECT || pqa == DISJOINT);
6516 pqb = edge_vert_col_inter(P, Q, B);
6519 assert(pqb != SHAREVERTEX);
6520 assert(pqb == INTERSECT || pqb == DISJOINT);
6558 enum tetgenmesh::interresult tetgenmesh::tri_vert_cop_inter(REAL* A, REAL* B,
6559 REAL* C, REAL* P, REAL* R)
6565 assert(R != (REAL *) NULL);
6569 s1 = orient3d(A, B, C, R);
6573 sign = s1 < 0.0 ? 1 : -1;
6576 s1 = orient3d(A, B, R, P) * sign;
6581 s2 = orient3d(B, C, R, P) * sign;
6586 s3 = orient3d(C, A, R, P) * sign;
6644 enum tetgenmesh::interresult tetgenmesh::tri_edge_cop_inter(REAL* A, REAL* B,
6645 REAL* C, REAL* P, REAL* Q, REAL* R)
6647 enum interresult abpq, bcpq, capq;
6648 enum interresult abcp, abcq;
6651 abpq = edge_edge_cop_inter(A, B, P, Q, R);
6652 if (abpq == INTERSECT || abpq == SHAREEDGE) {
6655 bcpq = edge_edge_cop_inter(B, C, P, Q, R);
6656 if (bcpq == INTERSECT || bcpq == SHAREEDGE) {
6659 capq = edge_edge_cop_inter(C, A, P, Q, R);
6660 if (capq == INTERSECT || capq == SHAREEDGE) {
6665 abcp = tri_vert_cop_inter(A, B, C, P, R);
6666 if (abcp == INTERSECT) {
6669 abcq = tri_vert_cop_inter(A, B, C, Q, R);
6670 if (abcq == INTERSECT) {
6676 if (abpq == SHAREVERTEX) {
6679 assert(abcp ^ abcq);
6683 if (bcpq == SHAREVERTEX) {
6686 assert(abcp ^ abcq);
6690 if (capq == SHAREVERTEX) {
6693 assert(abcp ^ abcq);
6716 enum tetgenmesh::interresult tetgenmesh::tri_edge_inter_tail(REAL* A, REAL* B,
6717 REAL* C, REAL* P, REAL* Q, REAL s1, REAL s2)
6722 if (s1 * s2 > 0.0) {
6727 if (s1 * s2 < 0.0) {
6731 sign = s1 < 0.0 ? 1 : -1;
6732 s3 = orient3d(A, B, P, Q) * sign;
6737 s4 = orient3d(B, C, P, Q) * sign;
6742 s5 = orient3d(C, A, P, Q) * sign;
6788 if (s1 != 0.0 || s2 != 0.0) {
6795 return tri_vert_cop_inter(A, B, C, P, Q);
6801 return tri_vert_cop_inter(A, B, C, Q, P);
6808 REAL ax, ay, az, bx, by, bz;
6816 N[0] = ay * bz - by * az;
6817 N[1] = az * bx - bz * ax;
6818 N[2] = ax * by - bx * ay;
6821 assert((fabs(N[0]) + fabs(N[1]) + fabs(N[2])) > 0.0);
6831 if (R[0] == A[0] && R[1] == A[1] && R[2] == A[2]) {
6833 for (i = 0; i < 3; i++) {
6835 assert (R[i] == A[i]);
6840 N[i] += (j * macheps);
6842 N[i] -= (j * macheps);
6846 }
while (R[i] == A[i]);
6850 return tri_edge_cop_inter(A, B, C, P, Q, R);
6863 enum tetgenmesh::interresult tetgenmesh::tri_edge_inter(REAL* A, REAL* B,
6864 REAL* C, REAL* P, REAL* Q)
6869 s1 = orient3d(A, B, C, P);
6870 s2 = orient3d(A, B, C, Q);
6872 return tri_edge_inter_tail(A, B, C, P, Q, s1, s2);
6885 enum tetgenmesh::interresult tetgenmesh::tri_tri_inter(REAL* A, REAL* B,
6886 REAL* C, REAL* O, REAL* P, REAL* Q)
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)) {
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)) {
6907 enum interresult abcop, abcpq, abcqo;
6910 abcop = tri_edge_inter_tail(A, B, C, O, P, s_o, s_p);
6911 if (abcop == INTERSECT) {
6913 }
else if (abcop == SHAREEDGE) {
6916 abcpq = tri_edge_inter_tail(A, B, C, P, Q, s_p, s_q);
6917 if (abcpq == INTERSECT) {
6919 }
else if (abcpq == SHAREEDGE) {
6922 abcqo = tri_edge_inter_tail(A, B, C, Q, O, s_q, s_o);
6923 if (abcqo == INTERSECT) {
6925 }
else if (abcqo == SHAREEDGE) {
6928 if (shareedge == 3) {
6934 assert(shareedge == 0 || shareedge == 1);
6938 enum interresult opqab, opqbc, opqca;
6940 opqab = tri_edge_inter_tail(O, P, Q, A, B, s_a, s_b);
6941 if (opqab == INTERSECT) {
6944 opqbc = tri_edge_inter_tail(O, P, Q, B, C, s_b, s_c);
6945 if (opqbc == INTERSECT) {
6948 opqca = tri_edge_inter_tail(O, P, Q, C, A, s_c, s_a);
6949 if (opqca == INTERSECT) {
6955 if (abcop == SHAREEDGE) {
6957 assert(abcpq == SHAREVERTEX && abcqo == SHAREVERTEX);
6962 if (abcpq == SHAREEDGE) {
6964 assert(abcop == SHAREVERTEX && abcqo == SHAREVERTEX);
6969 if (abcqo == SHAREEDGE) {
6971 assert(abcop == SHAREVERTEX && abcpq == SHAREVERTEX);
6978 if (abcop == SHAREVERTEX) {
6980 if (abcpq == SHAREVERTEX) {
6983 assert(abcqo != SHAREVERTEX);
6988 assert(abcqo == SHAREVERTEX);
6993 if (abcpq == SHAREVERTEX) {
6996 assert(abcqo == SHAREVERTEX);
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)
7023 det = insphere(pa, pb, pc, pd, pe);
7030 REAL sign, det_c, det_d;
7031 int idx[5], perm, tmp;
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;
7043 for (i = 0; i < n - 1; i++) {
7044 for (j = 0; j < n - 1 - i; j++) {
7045 if (idx[j + 1] < idx[j]) {
7047 idx[j] = idx[j + 1];
7057 sign = (perm % 2 == 0) ? 1.0 : -1.0;
7058 det_c = orient3d(p[1], p[2], p[3], p[4]);
7060 return sign * det_c;
7062 det_d = orient3d(p[0], p[2], p[3], p[4]);
7063 return -sign * det_d;
7076 bool tetgenmesh::iscollinear(REAL* A, REAL* B, REAL* C, REAL eps)
7093 Lv = abx * abx + aby * aby + abz * abz;
7095 if (Lv < q)
return true;
7096 Lw = acx * acx + acy * acy + acz * acz;
7098 if (Lw < q)
return true;
7099 dd = abx * acx + aby * acy + abz * acz;
7101 d = (dd * dd) / (Lv * Lw);
7102 if (d > 1.0) d = 1.0;
7120 iscoplanar(REAL* k, REAL* l, REAL* m, REAL* n, REAL vol6, REAL eps)
7125 if (vol6 == 0.0)
return true;
7130 L = sqrt(x * x + y * y + z * z);
7134 L += sqrt(x * x + y * y + z * z);
7138 L += sqrt(x * x + y * y + z * z);
7142 L += sqrt(x * x + y * y + z * z);
7146 L += sqrt(x * x + y * y + z * z);
7150 L += sqrt(x * x + y * y + z * z);
7155 q = fabs(vol6) / (L * L * L);
7172 iscospheric(REAL* k, REAL* l, REAL* m, REAL* n, REAL* o, REAL vol24, REAL eps)
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);
7191 q = fabs(vol24) / (L * L * L * L);
7205 inline REAL tetgenmesh::distance(REAL* p1, REAL* p2)
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]));
7226 REAL tetgenmesh::shortdistance(REAL* p, REAL* e1, REAL* e2)
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];
7238 len = sqrt(dot(v1, v1));
7247 return sqrt(dot(v2, v2) - l_p * l_p);
7256 REAL tetgenmesh::shortdistance(REAL* p, REAL* e1, REAL* e2, REAL* e3)
7260 projpt2face(p, e1, e2, e3, prj);
7261 return distance(p, prj);
7277 REAL tetgenmesh::interiorangle(REAL* o, REAL* p1, REAL* p2, REAL* n)
7279 REAL v1[3], v2[3], np[3];
7280 REAL theta, costheta, lenlen;
7281 REAL ori, len1, len2;
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;
7294 assert(lenlen != 0.0);
7296 costheta = dot(v1, v2) / lenlen;
7297 if (costheta > 1.0) {
7299 }
else if (costheta < -1.0) {
7302 theta = acos(costheta);
7305 np[0] = o[0] + n[0];
7306 np[1] = o[1] + n[1];
7307 np[2] = o[2] + n[2];
7309 ori = orient3d(p1, o, np, p2);
7311 theta = 2 * PI - theta;
7324 void tetgenmesh::projpt2edge(REAL* p, REAL* e1, REAL* e2, REAL* prj)
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];
7336 len = sqrt(dot(v1, v1));
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];
7356 void tetgenmesh::projpt2face(REAL* p, REAL* f1, REAL* f2, REAL* f3, REAL* prj)
7358 REAL fnormal[3], v1[3];
7362 facenormal(f1, f2, f3, fnormal, &len);
7370 v1[0] = p[0] - f1[0];
7371 v1[1] = p[1] - f1[1];
7372 v1[2] = p[2] - f1[2];
7374 dist = dot(fnormal, v1);
7377 prj[0] = p[0] - dist * fnormal[0];
7378 prj[1] = p[1] - dist * fnormal[1];
7379 prj[2] = p[2] - dist * fnormal[2];
7393 void tetgenmesh::facenormal(REAL* pa, REAL* pb, REAL* pc, REAL* n, REAL* nlen)
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];
7405 if (nlen != (REAL *) NULL) {
7406 *nlen = sqrt(dot(n, n));
7426 void tetgenmesh::edgeorthonormal(REAL* e1, REAL* e2, REAL* op, REAL* n)
7428 REAL v1[3], v2[3], fn[3];
7432 v1[0] = e2[0] - e1[0];
7433 v1[1] = e2[1] - e1[1];
7434 v1[2] = e2[2] - e1[2];
7436 v2[0] = op[0] - e1[0];
7437 v2[1] = op[1] - e1[1];
7438 v2[2] = op[2] - e1[2];
7444 len = sqrt(dot(n, n));
7461 REAL tetgenmesh::facedihedral(REAL* pa, REAL* pb, REAL* pc1, REAL* pc2)
7468 facenormal(pa, pb, pc1, n1, &n1len);
7469 facenormal(pa, pb, pc2, n2, &n2len);
7470 costheta = dot(n1, n2) / (n1len * n2len);
7472 if (costheta > 1.0) {
7474 }
else if (costheta < -1.0) {
7477 theta = acos(costheta);
7478 ori = orient3d(pa, pb, pc1, pc2);
7480 theta = 2 * PI - theta;
7497 void tetgenmesh::tetalldihedral(point pa, point pb, point pc, point pd,
7498 REAL* cosdd, REAL* cosmaxd, REAL* cosmind)
7500 REAL N[4][3], cosd, len;
7507 tetallnormal(pa, pb, pc, pd, N, NULL);
7509 for (i = 0; i < 4; i++) {
7510 len = sqrt(dot(N[i], N[i]));
7512 for (j = 0; j < 3; j++) N[i][j] /= len;
7516 for (i = 0; i < 6; i++) {
7518 case 0: f1 = 2; f2 = 3;
break;
7519 case 1: f1 = 0; f2 = 3;
break;
7520 case 2: f1 = 1; f2 = 3;
break;
7521 case 3: f1 = 1; f2 = 2;
break;
7522 case 4: f1 = 2; f2 = 0;
break;
7523 case 5: f1 = 0; f2 = 1;
break;
7525 cosd = -dot(N[f1], N[f2]);
7526 if (cosdd) cosdd[i] = cosd;
7528 if (cosmaxd) *cosmaxd = cosd;
7529 if (cosmind) *cosmind = cosd;
7531 if (cosmaxd) *cosmaxd = cosd < *cosmaxd ? cosd : *cosmaxd;
7532 if (cosmind) *cosmind = cosd > *cosmind ? cosd : *cosmind;
7546 void tetgenmesh::tetallnormal(point pa, point pb, point pc, point pd,
7547 REAL N[4][3], REAL* volume)
7549 REAL A[4][4], rhs[4], D;
7554 for (i = 0; i < 3; i++) A[0][i] = pa[i] - pd[i];
7555 for (i = 0; i < 3; i++) A[1][i] = pb[i] - pd[i];
7556 for (i = 0; i < 3; i++) A[2][i] = pc[i] - pd[i];
7558 lu_decmp(A, 3, indx, &D, 0);
7559 if (volume != NULL) {
7561 *volume = fabs((A[indx[0]][0] * A[indx[1]][1] * A[indx[2]][2])) / 6.0;
7563 for (j = 0; j < 3; j++) {
7564 for (i = 0; i < 3; i++) rhs[i] = 0.0;
7566 lu_solve(A, 3, indx, rhs, 0);
7567 for (i = 0; i < 3; i++) N[j][i] = rhs[i];
7570 for (i = 0; i < 3; i++) N[3][i] = - N[0][i] - N[1][i] - N[2][i];