47 char workstring[1024];
60 for (i = startindex; i < argc; i++) {
64 if (startindex == 1) {
66 if (argv[i][0] !=
'-') {
73 for (j = startindex; argv[i][j] !=
'\0'; j++) {
74 if (argv[i][j] ==
'p') {
76 if (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
77 (argv[i][j + 1] ==
'.')) {
79 while (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
80 (argv[i][j + 1] ==
'.')) {
82 workstring[k] = argv[i][j];
88 if ((argv[i][j + 1] ==
'/') || (argv[i][j + 1] ==
',')) {
90 if (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
91 (argv[i][j + 1] ==
'.')) {
93 while (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
94 (argv[i][j + 1] ==
'.') || (argv[i][j + 1] ==
'e') ||
95 (argv[i][j + 1] ==
'-') || (argv[i][j + 1] ==
'+')) {
97 workstring[k] = argv[i][j];
100 workstring[k] =
'\0';
104 if ((argv[i][j + 1] ==
'/') || (argv[i][j + 1] ==
',')) {
106 if (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
107 (argv[i][j + 1] ==
'.')) {
109 while (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
110 (argv[i][j + 1] ==
'.')) {
112 workstring[k] = argv[i][j];
115 workstring[k] =
'\0';
119 }
else if (argv[i][j] ==
's') {
121 }
else if (argv[i][j] ==
'Y') {
123 if ((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) {
127 if ((argv[i][j + 1] ==
'/') || (argv[i][j + 1] ==
',')) {
129 if ((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) {
134 if ((argv[i][j + 1] ==
'/') || (argv[i][j + 1] ==
',')) {
136 if ((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) {
141 }
else if (argv[i][j] ==
'r') {
143 }
else if (argv[i][j] ==
'q') {
145 if (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
146 (argv[i][j + 1] ==
'.')) {
148 while (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
149 (argv[i][j + 1] ==
'.')) {
151 workstring[k] = argv[i][j];
154 workstring[k] =
'\0';
157 if ((argv[i][j + 1] ==
'/') || (argv[i][j + 1] ==
',')) {
159 if (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
160 (argv[i][j + 1] ==
'.')) {
162 while (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
163 (argv[i][j + 1] ==
'.')) {
165 workstring[k] = argv[i][j];
168 workstring[k] =
'\0';
172 }
else if (argv[i][j] ==
'R') {
174 if ((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) {
178 if ((argv[i][j + 1] ==
'/') || (argv[i][j + 1] ==
',')) {
180 if (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
181 (argv[i][j + 1] ==
'.')) {
183 while (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
184 (argv[i][j + 1] ==
'.')) {
186 workstring[k] = argv[i][j];
189 workstring[k] =
'\0';
193 }
else if (argv[i][j] ==
'w') {
195 if ((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) {
199 }
else if (argv[i][j] ==
'b') {
202 if (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
203 (argv[i][j + 1] ==
'.')) {
205 while (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
206 (argv[i][j + 1] ==
'.')) {
208 workstring[k] = argv[i][j];
211 workstring[k] =
'\0';
214 if ((argv[i][j + 1] ==
'/') || (argv[i][j + 1] ==
',')) {
216 if (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
217 (argv[i][j + 1] ==
'.')) {
219 while (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
220 (argv[i][j + 1] ==
'.')) {
222 workstring[k] = argv[i][j];
225 workstring[k] =
'\0';
229 if ((argv[i][j + 1] ==
'/') || (argv[i][j + 1] ==
',')) {
231 if (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
232 (argv[i][j + 1] ==
'.') || (argv[i][j + 1] ==
'-')) {
234 while (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
235 (argv[i][j + 1] ==
'.') || (argv[i][j + 1] ==
'-')) {
237 workstring[k] = argv[i][j];
240 workstring[k] =
'\0';
244 if ((argv[i][j + 1] ==
'/') || (argv[i][j + 1] ==
',')) {
246 if (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
247 (argv[i][j + 1] ==
'.') || (argv[i][j + 1] ==
'-')) {
249 while (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
250 (argv[i][j + 1] ==
'.') || (argv[i][j + 1] ==
'-')) {
252 workstring[k] = argv[i][j];
255 workstring[k] =
'\0';
266 }
else if (argv[i][j] ==
'l') {
268 }
else if (argv[i][j] ==
'L') {
270 }
else if (argv[i][j] ==
'm') {
272 }
else if (argv[i][j] ==
'a') {
273 if (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
274 (argv[i][j + 1] ==
'.')) {
277 while (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
278 (argv[i][j + 1] ==
'.') || (argv[i][j + 1] ==
'e') ||
279 (argv[i][j + 1] ==
'-') || (argv[i][j + 1] ==
'+')) {
281 workstring[k] = argv[i][j];
284 workstring[k] =
'\0';
289 }
else if (argv[i][j] ==
'A') {
291 }
else if (argv[i][j] ==
'D') {
293 if ((argv[i][j + 1] >=
'1') && (argv[i][j + 1] <=
'3')) {
294 reflevel = (argv[i][j + 1] -
'1') + 1;
297 }
else if (argv[i][j] ==
'i') {
299 }
else if (argv[i][j] ==
'd') {
301 }
else if (argv[i][j] ==
'c') {
303 }
else if (argv[i][j] ==
'M') {
306 if ((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'1')) {
310 if ((argv[i][j + 1] ==
'/') || (argv[i][j + 1] ==
',')) {
312 if ((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'1')) {
317 }
else if (argv[i][j] ==
'X') {
318 if (argv[i][j + 1] ==
'1') {
324 }
else if (argv[i][j] ==
'z') {
325 if (argv[i][j + 1] ==
'1') {
331 }
else if (argv[i][j] ==
'f') {
333 }
else if (argv[i][j] ==
'e') {
335 }
else if (argv[i][j] ==
'n') {
337 }
else if (argv[i][j] ==
'v') {
339 }
else if (argv[i][j] ==
'g') {
341 }
else if (argv[i][j] ==
'k') {
343 }
else if (argv[i][j] ==
'J') {
345 }
else if (argv[i][j] ==
'B') {
347 }
else if (argv[i][j] ==
'N') {
349 }
else if (argv[i][j] ==
'E') {
351 }
else if (argv[i][j] ==
'F') {
353 }
else if (argv[i][j] ==
'I') {
355 }
else if (argv[i][j] ==
'S') {
356 if (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
357 (argv[i][j + 1] ==
'.')) {
359 while (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
360 (argv[i][j + 1] ==
'.') || (argv[i][j + 1] ==
'e') ||
361 (argv[i][j + 1] ==
'-') || (argv[i][j + 1] ==
'+')) {
363 workstring[k] = argv[i][j];
366 workstring[k] =
'\0';
367 steinerleft = (int) strtol(workstring, (
char **) NULL, 0);
369 }
else if (argv[i][j] ==
'o') {
370 if (argv[i][j + 1] ==
'2') {
374 if ((argv[i][j + 1] ==
'/') || (argv[i][j + 1] ==
',')) {
376 if (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
377 (argv[i][j + 1] ==
'.')) {
379 while (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
380 (argv[i][j + 1] ==
'.')) {
382 workstring[k] = argv[i][j];
385 workstring[k] =
'\0';
389 }
else if (argv[i][j] ==
'O') {
390 if ((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) {
394 if ((argv[i][j + 1] ==
'/') || (argv[i][j + 1] ==
',')) {
396 if ((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'7')) {
401 }
else if (argv[i][j] ==
'T') {
402 if (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
403 (argv[i][j + 1] ==
'.')) {
405 while (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
406 (argv[i][j + 1] ==
'.') || (argv[i][j + 1] ==
'e') ||
407 (argv[i][j + 1] ==
'-') || (argv[i][j + 1] ==
'+')) {
409 workstring[k] = argv[i][j];
412 workstring[k] =
'\0';
413 epsilon = (
REAL) strtod(workstring, (
char **) NULL);
415 }
else if (argv[i][j] ==
'C') {
417 }
else if (argv[i][j] ==
'Q') {
419 }
else if (argv[i][j] ==
'V') {
421 }
else if (argv[i][j] ==
'x') {
422 if (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
423 (argv[i][j + 1] ==
'.')) {
425 while (((argv[i][j + 1] >=
'0') && (argv[i][j + 1] <=
'9')) ||
426 (argv[i][j + 1] ==
'.') || (argv[i][j + 1] ==
'e') ||
427 (argv[i][j + 1] ==
'-') || (argv[i][j + 1] ==
'+')) {
429 workstring[k] = argv[i][j];
432 workstring[k] =
'\0';
441 }
else if ((argv[i][j] ==
'h') || (argv[i][j] ==
'H') ||
442 (argv[i][j] ==
'?')) {
444 printf(
"Warning: Unknown switch -%c.\n", argv[i][j]);
449 if (startindex == 0) {
516 printf(
"Error: Switches -w cannot use together with -p or -r.\n");
567 while (workstring[j] !=
'\0') {
568 if ((workstring[j] ==
'.') && (workstring[j + 1] !=
'\0')) {
577 if ((workstring[j] >=
'0') && (workstring[j] <=
'9')) {
578 meshnumber = meshnumber * 10 + (int) (workstring[j] -
'0');
583 }
while (workstring[j] !=
'\0');
587 }
else if (increment == 0) {
591 workstring[increment] =
'%';
592 workstring[increment + 1] =
'd';
593 workstring[increment + 2] =
'\0';
634 int tetgenmesh::esymtbl[12] = {9, 6, 11, 4, 3, 7, 1, 5, 10, 0, 8, 2};
640 int tetgenmesh:: orgpivot[12] = {7, 7, 5, 5, 6, 4, 4, 6, 5, 6, 7, 4};
641 int tetgenmesh::destpivot[12] = {6, 4, 4, 6, 5, 6, 7, 4, 7, 7, 5, 5};
642 int tetgenmesh::apexpivot[12] = {5, 6, 7, 4, 7, 7, 5, 5, 6, 4, 4, 6};
643 int tetgenmesh::oppopivot[12] = {4, 5, 6, 7, 4, 5, 6, 7, 4, 5, 6, 7};
648 int tetgenmesh::ver2edge[12] = {0, 1, 2, 3, 3, 5, 1, 5, 4, 0, 4, 2};
653 int tetgenmesh::epivot[12] = {4, 5, 2, 11, 4, 5, 2, 11, 4, 5, 2, 11};
677 int soffset, toffset;
682 for (i = 0; i < 12; i++) {
683 for (j = 0; j < 12; j++) {
684 bondtbl[i][j] = (j & 3) + (((i & 12) + (j & 12)) % 12);
690 for (i = 0; i < 12; i++) {
691 for (j = 0; j < 12; j++) {
692 fsymtbl[i][j] = (j + 12 - (i & 12)) % 12;
697 for (i = 0; i < 12; i++) {
701 for (i = 0; i < 12; i++) {
702 for (j = 0; j < 12; j++) {
707 for (i = 0; i < 12; i++) {
712 for (i = 0; i < 12; i++) {
717 for (i = 0; i < 12; i++) {
723 for (i = 0; i < 12; i++) {
724 for (j = 0; j < 6; j++) {
726 soffset = (6 - ((i & 12) >> 1)) % 6;
727 toffset = (12 - ((j & 6) << 1)) % 12;
729 soffset = (i & 12) >> 1;
730 toffset = (j & 6) << 1;
732 tsbondtbl[i][j] = (j & 1) + (((j & 6) + soffset) % 6);
733 stbondtbl[i][j] = (i & 3) + (((i & 12) + toffset) % 12);
739 for (i = 0; i < 12; i++) {
740 for (j = 0; j < 6; j++) {
742 soffset = (i & 12) >> 1;
743 toffset = (j & 6) << 1;
745 soffset = (6 - ((i & 12) >> 1)) % 6;
746 toffset = (12 - ((j & 6) << 1)) % 12;
748 tspivottbl[i][j] = (j & 1) + (((j & 6) + soffset) % 6);
749 stpivottbl[i][j] = (i & 3) + (((i & 12) + toffset) % 12);
782 objectbytes = sizeofobject > 1 ? sizeofobject : 1;
784 log2objectsperblock = log2objperblk;
786 objectsperblock = ((int) 1) << log2objectsperblock;
787 objectsperblockmark = objectsperblock - 1;
792 toparray = (
char **) NULL;
807 poolinit(sizeofobject, log2objperblk);
815 if (toparray != (
char **) NULL) {
817 for (i = 0; i < toparraylen; i++) {
819 if (toparray[i] != (
char *) NULL) {
821 free((
void *) toparray[i]);
825 free((
void *) toparray);
829 toparray = (
char **) NULL;
856 topindex = objectindex >> log2objectsperblock;
858 if (toparray == (
char **) NULL) {
861 newsize = topindex + 128;
862 toparray = (
char **) malloc((
size_t) (newsize *
sizeof(
char *)));
863 toparraylen = newsize;
864 for (i = 0; i < newsize; i++) {
865 toparray[i] = (
char *) NULL;
868 totalmemory = newsize * (uintptr_t)
sizeof(
char *);
869 }
else if (topindex >= toparraylen) {
871 newsize = 3 * toparraylen;
872 if (topindex >= newsize) {
873 newsize = topindex + 128;
877 newarray = (
char **) malloc((
size_t) (newsize *
sizeof(
char *)));
878 for (i = 0; i < toparraylen; i++) {
879 newarray[i] = toparray[i];
881 for (i = toparraylen; i < newsize; i++) {
882 newarray[i] = (
char *) NULL;
886 totalmemory += (newsize - toparraylen) *
sizeof(
char *);
888 toparraylen = newsize;
892 block = toparray[topindex];
893 if (block == (
char *) NULL) {
895 block = (
char *) malloc((
size_t) (objectsperblock * objectbytes));
896 toparray[topindex] = block;
898 totalmemory += objectsperblock * objectbytes;
918 if (toparray == (
char **) NULL) {
919 return (
void *) NULL;
923 topindex = objectindex >> log2objectsperblock;
925 if (topindex >= toparraylen) {
926 return (
void *) NULL;
930 block = toparray[topindex];
931 if (block == (
char *) NULL) {
932 return (
void *) NULL;
938 return (
void *)(block + (objectindex & (objectsperblock - 1)) * objectbytes);
952 int newindex = objects;
953 *newptr = (
void *) (getblock(objects) +
954 (objects & (objectsperblock - 1)) * objectbytes);
969 firstblock = nowblock = (
void **) NULL;
970 nextitem = (
void *) NULL;
971 deaditemstack = (
void *) NULL;
972 pathblock = (
void **) NULL;
973 pathitem = (
void *) NULL;
975 itembytes = itemwords = 0;
977 items = maxitems = 0l;
978 unallocateditems = 0;
985 poolinit(bytecount, itemcount, wsize, alignment);
996 while (firstblock != (
void **) NULL) {
997 nowblock = (
void **) *(firstblock);
999 firstblock = nowblock;
1027 if (alignment > wordsize) {
1028 alignbytes = alignment;
1030 alignbytes = wordsize;
1032 if ((
int)
sizeof(
void *) > alignbytes) {
1033 alignbytes = (int)
sizeof(
void *);
1035 itemwords = ((bytecount + alignbytes - 1) / alignbytes)
1036 * (alignbytes / wordsize);
1037 itembytes = itemwords * wordsize;
1038 itemsperblock = itemcount;
1043 firstblock = (
void **) malloc(itemsperblock * itembytes +
sizeof(
void *)
1045 if (firstblock == (
void **) NULL) {
1049 *(firstblock) = (
void *) NULL;
1071 nowblock = firstblock;
1073 alignptr = (uintptr_t) (nowblock + 1);
1076 (alignptr + (uintptr_t) alignbytes -
1077 (alignptr % (uintptr_t) alignbytes));
1079 unallocateditems = itemsperblock;
1081 deaditemstack = (
void *) NULL;
1098 if (deaditemstack != (
void *) NULL) {
1099 newitem = deaditemstack;
1100 deaditemstack = * (
void **) deaditemstack;
1103 if (unallocateditems == 0) {
1105 if (*nowblock == (
void *) NULL) {
1107 newblock = (
void **) malloc(itemsperblock * itembytes +
sizeof(
void *)
1109 if (newblock == (
void **) NULL) {
1112 *nowblock = (
void *) newblock;
1114 *newblock = (
void *) NULL;
1117 nowblock = (
void **) *nowblock;
1120 alignptr = (uintptr_t) (nowblock + 1);
1123 (alignptr + (uintptr_t) alignbytes -
1124 (alignptr % (uintptr_t) alignbytes));
1126 unallocateditems = itemsperblock;
1131 nextitem = (
void *) ((uintptr_t) nextitem + itembytes);
1150 *((
void **) dyingitem) = deaditemstack;
1151 deaditemstack = dyingitem;
1168 pathblock = firstblock;
1170 alignptr = (uintptr_t) (pathblock + 1);
1173 (alignptr + (uintptr_t) alignbytes -
1174 (alignptr % (uintptr_t) alignbytes));
1176 pathitemsleft = itemsperblock;
1197 if (pathitem == nextitem) {
1198 return (
void *) NULL;
1201 if (pathitemsleft == 0) {
1203 pathblock = (
void **) *pathblock;
1205 alignptr = (uintptr_t) (pathblock + 1);
1208 (alignptr + (uintptr_t) alignbytes -
1209 (alignptr % (uintptr_t) alignbytes));
1211 pathitemsleft = itemsperblock;
1215 pathitem = (
void *) ((uintptr_t) pathitem + itembytes);
1236 printf(
" Constructing mapping from indices to points.\n");
1243 idx =
in->firstnumber;
1244 while (pointloop != (
point) NULL) {
1245 idx2verlist[idx++] = pointloop;
1265 face*& facperverlist)
1271 printf(
" Making a map from points to subfaces.\n");
1276 for (i = 0; i <
points->
items + 1; i++) idx2faclist[i] = 0;
1288 if (shloop.
sh[5] != NULL) {
1299 k = idx2faclist[i + 1];
1300 idx2faclist[i + 1] = idx2faclist[i] + j;
1305 facperverlist =
new face[idx2faclist[i]];
1313 facperverlist[idx2faclist[j]] = shloop;
1316 if (shloop.
sh[5] != NULL) {
1319 facperverlist[idx2faclist[j]] = shloop;
1323 facperverlist[idx2faclist[j]] = shloop;
1328 facperverlist[idx2faclist[j]] = shloop;
1336 idx2faclist[i + 1] = idx2faclist[i];
1354 if (dyingtetrahedron[8] != NULL) {
1357 if (dyingtetrahedron[9] != NULL) {
1379 }
while ((newtetrahedron[4] == (
tetrahedron) NULL) ||
1381 return newtetrahedron;
1393 }
while (newtetrahedron[4] == (
tetrahedron) NULL);
1394 return newtetrahedron;
1409 pool->
dealloc((
void *) dyingsh);
1425 if (newshellface == (
shellface *) NULL) {
1428 }
while (newshellface[3] == (
shellface) NULL);
1429 return newshellface;
1459 if (newpoint == (
point) NULL) {
1460 return (
point) NULL;
1477 newtet->
tet[0] = NULL;
1478 newtet->
tet[1] = NULL;
1479 newtet->
tet[2] = NULL;
1480 newtet->
tet[3] = NULL;
1482 newtet->
tet[4] = NULL;
1483 newtet->
tet[5] = NULL;
1484 newtet->
tet[6] = NULL;
1485 newtet->
tet[7] = NULL;
1487 newtet->
tet[8] = NULL;
1488 newtet->
tet[9] = NULL;
1514 newface->
sh[0] = NULL;
1515 newface->
sh[1] = NULL;
1516 newface->
sh[2] = NULL;
1518 newface->
sh[3] = NULL;
1519 newface->
sh[4] = NULL;
1520 newface->
sh[5] = NULL;
1522 newface->
sh[6] = NULL;
1523 newface->
sh[7] = NULL;
1524 newface->
sh[8] = NULL;
1526 newface->
sh[9] = NULL;
1527 newface->
sh[10] = NULL;
1557 (*pnewpoint)[3 + i] = 0.0;
1593 int pointsize = 0, elesize = 0, shsize = 0;
1597 printf(
" Initializing memorypools.\n");
1612 if (
addin != NULL) {
1626 if (
in->segmentconstraintlist ||
in->facetconstraintlist) {
1654 bgm->
in->numberofpointmtrs :
in->numberofpointmtrs;
1788 printf(
" Size of a tetrahedron: %d (%d) bytes.\n", elesize,
1809 shmarkindex = (shsize +
sizeof(int) - 1) /
sizeof(int);
1821 printf(
" Size of a shellface: %d (%d) bytes.\n", shsize,
1894 sign =
insphere(pa, pb, pc, pd, pe);
1900 point pt[5], swappt;
1919 for (i = 0; i < n; i++) {
1921 swappt = pt[i]; pt[i] = pt[i+1]; pt[i+1] = swappt;
1926 }
while (count > 0);
1928 oriA =
orient3d(pt[1], pt[2], pt[3], pt[4]);
1931 if ((swaps % 2) != 0) oriA = -oriA;
1935 oriB = -
orient3d(pt[0], pt[2], pt[3], pt[4]);
1940 if ((swaps % 2) != 0) oriB = -oriB;
1967 sign =
orient4d(pa, pb, pc, pd, pe,
1968 aheight, bheight, cheight, dheight, eheight);
1974 point pt[5], swappt;
1993 for (i = 0; i < n; i++) {
1995 swappt = pt[i]; pt[i] = pt[i+1]; pt[i+1] = swappt;
2000 }
while (count > 0);
2002 oriA =
orient3d(pt[1], pt[2], pt[3], pt[4]);
2005 if ((swaps % 2) != 0) oriA = -oriA;
2009 oriB = -
orient3d(pt[0], pt[2], pt[3], pt[4]);
2014 if ((swaps % 2) != 0) oriB = -oriB;
2040 #define SETVECTOR3(V, a0, a1, a2) (V)[0] = (a0); (V)[1] = (a1); (V)[2] = (a2)
2042 #define SWAP2(a0, a1, tmp) (tmp) = (a0); (a0) = (a1); (a1) = (tmp)
2045 point R,
int level,
int *types,
int *pos)
2051 REAL s1, s2, s3, s4;
2060 len = sqrt(
dot(n, n));
2070 R[0] = A[0] + len * n[0];
2071 R[1] = A[1] + len * n[1];
2072 R[2] = A[2] + len * n[2];
2297 s1 =
orient3d(U[0], U[2], R, V[1]);
2298 s2 =
orient3d(U[1], U[2], R, V[0]);
2337 s3 =
orient3d(U[0], U[2], R, V[0]);
2338 s4 =
orient3d(U[1], U[2], R, V[1]);
2444 }
else if (z1 == 2) {
2548 }
else if (z1 == 3) {
2656 REAL sP,
REAL sQ,
int level,
int *types,
int *pos)
2725 return tri_edge_2d(A, B, C, P, Q, R, level, types, pos);
2728 s1 =
orient3d(U[0], U[1], V[0], V[1]);
2733 s2 =
orient3d(U[1], U[2], V[0], V[1]);
2738 s3 =
orient3d(U[2], U[0], V[0], V[1]);
2854 point R,
int level,
int *types,
int *pos)
2862 return tri_edge_tail(A, B, C, P, Q, R, sP, sQ, level, types, pos);
2878 int types[2], pos[4];
2881 ni =
tri_edge_tail(A, B, C, P, Q, NULL, s_p, s_q, 1, types, pos);
2891 }
else if (ni == 4) {
2920 if ((s_o * s_p > 0.0) && (s_o * s_q > 0.0)) {
2928 if ((s_a * s_b > 0.0) && (s_a * s_c > 0.0)) {
2933 int abcop, abcpq, abcqo;
2954 if (shareedge == 3) {
2960 int opqab, opqbc, opqca;
3036 for (i = N; i < n + N; i++) {
3039 for (j = N; j < n + N; j++)
3040 if (biggest < (tempf = fabs(lu[i][j])))
3043 scales[i] = 1.0 / biggest;
3051 for (k = N; k < n + N - 1; k++) {
3054 for (i = k; i < n + N; i++) {
3055 if (biggest < (tempf = fabs(lu[ps[i]][k]) * scales[ps[i]])) {
3060 if (biggest == 0.0) {
3063 if (pivotindex != k) {
3065 ps[k] = ps[pivotindex];
3071 pivot = lu[ps[k]][k];
3072 for (i = k + 1; i < n + N; i++) {
3073 lu[ps[i]][k] =
mult = lu[ps[i]][k] / pivot;
3075 for (j = k + 1; j < n + N; j++)
3076 lu[ps[i]][j] -=
mult * lu[ps[k]][j];
3082 return lu[ps[n + N - 1]][n + N - 1] != 0.0;
3105 for (i = N; i < n + N; i++) X[i] = 0.0;
3108 for (i = N; i < n + N; i++) {
3110 for (j = N; j < i + N; j++)
3111 dot += lu[ps[i]][j] * X[j];
3112 X[i] =
b[ps[i]] -
dot;
3116 for (i = n + N - 1; i >= N; i--) {
3118 for (j = i + 1; j < n + N; j++)
3119 dot += lu[ps[i]][j] * X[j];
3120 X[i] = (X[i] -
dot) / lu[ps[i]][i];
3123 for (i = N; i < n + N; i++)
b[i] = X[i];
3140 REAL area2[2], n1[3], n2[3],
c[3];
3145 area2[0] =
dot(n1, n1);
3147 area2[1] =
dot(n2, n2);
3149 if (area2[0] > area2[1]) {
3165 if (fabs(sign) / r < b->
epsilon) {
3192 REAL v1[3], v2[3], v3[3], *pv1, *pv2;
3195 v1[0] = pb[0] - pa[0];
3196 v1[1] = pb[1] - pa[1];
3197 v1[2] = pb[2] - pa[2];
3198 v2[0] = pa[0] - pc[0];
3199 v2[1] = pa[1] - pc[1];
3200 v2[2] = pa[2] - pc[2];
3205 v3[0] = pc[0] - pb[0];
3206 v3[1] = pc[1] - pb[1];
3207 v3[2] = pc[2] - pb[2];
3227 *lav = (sqrt(L1) + sqrt(L2) + sqrt(L3)) / 3.0;
3260 v1[0] = e2[0] - e1[0];
3261 v1[1] = e2[1] - e1[1];
3262 v1[2] = e2[2] - e1[2];
3263 v2[0] = p[0] - e1[0];
3264 v2[1] = p[1] - e1[1];
3265 v2[2] = p[2] - e1[2];
3267 len = sqrt(
dot(v1, v1));
3274 return sqrt(
dot(v2, v2) - l_p * l_p);
3288 A[0][0] = pb[0] - pa[0];
3289 A[0][1] = pb[1] - pa[1];
3290 A[0][2] = pb[2] - pa[2];
3291 A[1][0] = pc[0] - pa[0];
3292 A[1][1] = pc[1] - pa[1];
3293 A[1][2] = pc[2] - pa[2];
3295 cross(A[0], A[1], A[2]);
3297 return 0.5 * sqrt(
dot(A[2], A[2]));
3306 adx = pa[0] - pd[0];
3307 bdx = pb[0] - pd[0];
3308 cdx = pc[0] - pd[0];
3309 ady = pa[1] - pd[1];
3310 bdy = pb[1] - pd[1];
3311 cdy = pc[1] - pd[1];
3312 adz = pa[2] - pd[2];
3313 bdz = pb[2] - pd[2];
3314 cdz = pc[2] - pd[2];
3316 return adx * (bdy * cdz - bdz * cdy)
3317 + bdx * (cdy * adz - cdz * ady)
3318 + cdx * (ady * bdz - adz * bdy);
3336 REAL v1[3], v2[3], np[3];
3337 REAL theta, costheta, lenlen;
3338 REAL ori, len1, len2;
3341 v1[0] = p1[0] - o[0];
3342 v1[1] = p1[1] - o[1];
3343 v1[2] = p1[2] - o[2];
3344 v2[0] = p2[0] - o[0];
3345 v2[1] = p2[1] - o[1];
3346 v2[2] = p2[2] - o[2];
3347 len1 = sqrt(
dot(v1, v1));
3348 len2 = sqrt(
dot(v2, v2));
3349 lenlen = len1 * len2;
3351 costheta =
dot(v1, v2) / lenlen;
3352 if (costheta > 1.0) {
3354 }
else if (costheta < -1.0) {
3357 theta = acos(costheta);
3360 np[0] = o[0] + n[0];
3361 np[1] = o[1] + n[1];
3362 np[2] = o[2] + n[2];
3366 theta = 2 *
PI - theta;
3384 v1[0] = e2[0] - e1[0];
3385 v1[1] = e2[1] - e1[1];
3386 v1[2] = e2[2] - e1[2];
3387 v2[0] = p[0] - e1[0];
3388 v2[1] = p[1] - e1[1];
3389 v2[2] = p[2] - e1[2];
3391 len = sqrt(
dot(v1, v1));
3397 prj[0] = e1[0] + l_p * v1[0];
3398 prj[1] = e1[1] + l_p * v1[1];
3399 prj[2] = e1[2] + l_p * v1[2];
3410 REAL fnormal[3], v1[3];
3415 len = sqrt(fnormal[0]*fnormal[0] + fnormal[1]*fnormal[1] +
3416 fnormal[2]*fnormal[2]);
3421 v1[0] = p[0] - f1[0];
3422 v1[1] = p[1] - f1[1];
3423 v1[2] = p[2] - f1[2];
3425 dist =
dot(fnormal, v1);
3428 prj[0] = p[0] - dist * fnormal[0];
3429 prj[1] = p[1] - dist * fnormal[1];
3430 prj[2] = p[2] - dist * fnormal[2];
3447 REAL N[4][3], vol, cosd, len;
3448 int f1 = 0, f2 = 0, i, j;
3457 for (i = 0; i < 4; i++) {
3458 len = sqrt(
dot(N[i], N[i]));
3460 for (j = 0; j < 3; j++) N[i][j] /= len;
3475 for (i = 0; i < 4; i++) {
3476 len = sqrt(
dot(N[i], N[i]));
3478 for (j = 0; j < 3; j++) N[i][j] /= len;
3488 if (cosdd != NULL) {
3489 for (i = 0; i < 6; i++) {
3494 if (cosmaxd != NULL) {
3497 if (cosmind != NULL) {
3505 for (i = 0; i < 6; i++) {
3507 case 0: f1 = 0; f2 = 1;
break;
3508 case 1: f1 = 1; f2 = 2;
break;
3509 case 2: f1 = 2; f2 = 3;
break;
3510 case 3: f1 = 0; f2 = 3;
break;
3511 case 4: f1 = 2; f2 = 0;
break;
3512 case 5: f1 = 1; f2 = 3;
break;
3514 cosd = -
dot(N[f1], N[f2]);
3515 if (cosd < -1.0) cosd = -1.0;
3516 if (cosd > 1.0) cosd = 1.0;
3517 if (cosdd) cosdd[i] = cosd;
3518 if (cosmaxd || cosmind) {
3520 if (cosmaxd) *cosmaxd = cosd;
3521 if (cosmind) *cosmind = cosd;
3523 if (cosmaxd) *cosmaxd = cosd < *cosmaxd ? cosd : *cosmaxd;
3524 if (cosmind) *cosmind = cosd > *cosmind ? cosd : *cosmind;
3545 REAL A[4][4], rhs[4],
D;
3550 for (i = 0; i < 3; i++) A[0][i] = pa[i] - pd[i];
3551 for (i = 0; i < 3; i++) A[1][i] = pb[i] - pd[i];
3552 for (i = 0; i < 3; i++) A[2][i] = pc[i] - pd[i];
3556 if (volume != NULL) {
3558 *volume = fabs((A[indx[0]][0] * A[indx[1]][1] * A[indx[2]][2])) / 6.0;
3560 for (j = 0; j < 3; j++) {
3561 for (i = 0; i < 3; i++) rhs[i] = 0.0;
3564 for (i = 0; i < 3; i++) N[j][i] = rhs[i];
3567 for (i = 0; i < 3; i++) N[3][i] = - N[0][i] - N[1][i] - N[2][i];
3570 if (volume != NULL) {
3587 REAL V[6][3], edgelength[6], longlen;
3588 REAL vda[3], vdb[3], vdc[3];
3589 REAL N[4][3], A[4][4], rhs[4],
D;
3590 REAL H[4], volume, minheightinv;
3595 for (i = 0; i < 3; i++) V[0][i] = pa[i] - pd[i];
3596 for (i = 0; i < 3; i++) V[1][i] = pb[i] - pd[i];
3597 for (i = 0; i < 3; i++) V[2][i] = pc[i] - pd[i];
3598 for (i = 0; i < 3; i++) V[3][i] = pb[i] - pa[i];
3599 for (i = 0; i < 3; i++) V[4][i] = pc[i] - pb[i];
3600 for (i = 0; i < 3; i++) V[5][i] = pa[i] - pc[i];
3603 for (i = 0; i < 6; i++) edgelength[i] =
dot(V[i], V[i]);
3606 longlen = edgelength[0];
3607 for (i = 1; i < 6; i++) {
3608 longlen = edgelength[i] > longlen ? edgelength[i] : longlen;
3612 for (i = 0; i < 3; i++) A[0][i] = vda[i] = pa[i] - pd[i];
3613 for (i = 0; i < 3; i++) A[1][i] = vdb[i] = pb[i] - pd[i];
3614 for (i = 0; i < 3; i++) A[2][i] = vdc[i] = pc[i] - pd[i];
3618 volume = (A[indx[0]][0] * A[indx[1]][1] * A[indx[2]][2]) / 6.0;
3620 if (volume == 0.0)
return 1.0e+200;
3623 for (j = 0; j < 3; j++) {
3624 for (i = 0; i < 3; i++) rhs[i] = 0.0;
3627 for (i = 0; i < 3; i++) N[j][i] = rhs[i];
3630 for (i = 0; i < 3; i++) N[3][i] = - N[0][i] - N[1][i] - N[2][i];
3632 for (i = 0; i < 4; i++) {
3634 H[i] = sqrt(
dot(N[i], N[i]));
3642 minheightinv = H[0];
3643 for (i = 1; i < 4; i++) {
3644 if (H[i] > minheightinv) minheightinv = H[i];
3647 return sqrt(longlen) * minheightinv;
3668 REAL A[4][4], rhs[4],
D;
3672 A[0][0] = pb[0] - pa[0];
3673 A[0][1] = pb[1] - pa[1];
3674 A[0][2] = pb[2] - pa[2];
3675 A[1][0] = pc[0] - pa[0];
3676 A[1][1] = pc[1] - pa[1];
3677 A[1][2] = pc[2] - pa[2];
3679 A[2][0] = pd[0] - pa[0];
3680 A[2][1] = pd[1] - pa[1];
3681 A[2][2] = pd[2] - pa[2];
3683 cross(A[0], A[1], A[2]);
3687 rhs[0] = 0.5 *
dot(A[0], A[0]);
3688 rhs[1] = 0.5 *
dot(A[1], A[1]);
3690 rhs[2] = 0.5 *
dot(A[2], A[2]);
3698 if (radius != (
REAL *) NULL) *radius = 0.0;
3702 if (cent != (
REAL *) NULL) {
3703 cent[0] = pa[0] + rhs[0];
3704 cent[1] = pa[1] + rhs[1];
3705 cent[2] = pa[2] + rhs[2];
3707 if (radius != (
REAL *) NULL) {
3708 *radius = sqrt(rhs[0] * rhs[0] + rhs[1] * rhs[1] + rhs[2] * rhs[2]);
3727 REAL A[4][4], rhs[4],
D;
3731 A[0][0] = 1.0; A[0][1] = pa[0]; A[0][2] = pa[1]; A[0][3] = pa[2];
3732 A[1][0] = 1.0; A[1][1] = pb[0]; A[1][2] = pb[1]; A[1][3] = pb[2];
3733 A[2][0] = 1.0; A[2][1] = pc[0]; A[2][2] = pc[1]; A[2][3] = pc[2];
3734 A[3][0] = 1.0; A[3][1] = pd[0]; A[3][2] = pd[1]; A[3][3] = pd[2];
3737 rhs[0] = 0.5 * aheight;
3738 rhs[1] = 0.5 * bheight;
3739 rhs[2] = 0.5 * cheight;
3740 rhs[3] = 0.5 * dheight;
3745 if (radius != (
REAL *) NULL) *radius = 0.0;
3750 if (orthocent != (
REAL *) NULL) {
3751 orthocent[0] = rhs[1];
3752 orthocent[1] = rhs[2];
3753 orthocent[2] = rhs[3];
3755 if (radius != (
REAL *) NULL) {
3761 *radius = sqrt(rhs[1] * rhs[1] + rhs[2] * rhs[2] + rhs[3] * rhs[3]
3793 REAL n[3], det, det1;
3798 det = n[0] * (e2[0] - e1[0]) + n[1] * (e2[1] - e1[1])
3799 + n[2] * (e2[2] - e1[2]);
3802 det1 = n[0] * (pa[0] - e1[0]) + n[1] * (pa[1] - e1[1])
3803 + n[2] * (pa[2] - e1[2]);
3805 ip[0] = e1[0] + *u * (e2[0] - e1[0]);
3806 ip[1] = e1[1] + *u * (e2[1] - e1[1]);
3807 ip[2] = e1[2] + *u * (e2[2] - e1[2]);
3828 REAL vab[3], vcd[3], vca[3];
3829 REAL vab_vab, vcd_vcd, vab_vcd;
3830 REAL vca_vab, vca_vcd;
3834 for (i = 0; i < 3; i++) {
3835 vab[i] = B[i] - A[i];
3836 vcd[i] =
D[i] - C[i];
3837 vca[i] = A[i] - C[i];
3840 vab_vab =
dot(vab, vab);
3841 vcd_vcd =
dot(vcd, vcd);
3842 vab_vcd =
dot(vab, vcd);
3844 det = vab_vab * vcd_vcd - vab_vcd * vab_vcd;
3846 eps = det / (fabs(vab_vab * vcd_vcd) + fabs(vab_vcd * vab_vcd));
3851 vca_vab =
dot(vca, vab);
3852 vca_vcd =
dot(vca, vcd);
3854 *tp = (vcd_vcd * (- vca_vab) + vab_vcd * vca_vcd) / det;
3855 *tq = (vab_vcd * (- vca_vab) + vab_vab * vca_vcd) / det;
3857 for (i = 0; i < 3; i++) P[i] = A[i] + (*tp) * vab[i];
3858 for (i = 0; i < 3; i++) Q[i] = C[i] + (*tq) * vcd[i];
3886 REAL *p4, *p5, *p6, *p7;
3887 REAL w4, w5, w6, w7;
3902 vol[0] =
orient4d(p5, p6, p4, p3, p7, w5, w6, w4, 0, w7);
3903 vol[1] =
orient4d(p3, p6, p2, p0, p1, 0, w6, 0, 0, 0);
3904 vol[2] =
orient4d(p4, p6, p3, p0, p1, w4, w6, 0, 0, 0);
3905 vol[3] =
orient4d(p6, p5, p4, p3, p1, w6, w5, w4, 0, 0);
3907 return fabs(vol[0]) + fabs(vol[1]) + fabs(vol[2]) + fabs(vol[3]);
3919 point *ppt, pa, pb, pc;
3920 REAL v1[3], v2[3], n[3];
3921 REAL lab, len, A, area;
3931 for (i = 1; i < facpoints->
objects; i++) {
3933 x = (*ppt)[0] - pa[0];
3934 y = (*ppt)[1] - pa[1];
3935 z = (*ppt)[2] - pa[2];
3936 len = x * x + y * y +
z *
z;
3945 printf(
"Warning: All points of a facet are coincident with %d.\n",
3952 v1[0] = pb[0] - pa[0];
3953 v1[1] = pb[1] - pa[1];
3954 v1[2] = pb[2] - pa[2];
3956 for (i = 1; i < facpoints->
objects; i++) {
3958 v2[0] = (*ppt)[0] - pa[0];
3959 v2[1] = (*ppt)[1] - pa[1];
3960 v2[2] = (*ppt)[2] - pa[2];
3971 printf(
"Warning: All points of a facet are collinaer with [%d, %d].\n",
3979 len = sqrt(
dot(n, n));
4009 REAL len, len1, len2;
4013 len1 = sqrt(
dot(n1, n1));
4015 len2 = sqrt(
dot(n2, n2));
4044 point pa, pb, pc, pd;
4052 printf(
"Found two %s self-intersecting facets (dihedral angle %12.5E).\n",
4053 dihedang > 0 ?
"nearly" :
"exactly",dihedang);
4054 printf(
" 1st: [%d, %d, %d] #%d\n",
4056 printf(
" 2nd: [%d, %d, %d] #%d\n",
4059 printf(
"The dihedral angle between them is %g degree.\n",
4060 dihedang /
PI * 180.0);
4063 printf(
"Hint: You may use Mesh.AngleToleranceFacetOverlap to decrease the dihedral angle tolerance %g (degree)",
4069 printf(
"Found two overlapping facets.\n");
4071 printf(
"Found two duplicated facets.\n");
4073 printf(
" 1st: [%d, %d, %d] #%d\n",
4075 printf(
" 2nd: [%d, %d, %d] #%d\n",
4112 point forg = NULL, fdest = NULL, fapex = NULL;
4113 int etype = 0, geomtag = 0, facemark = 0;
4115 if (iedge != NULL) {
4116 if (iedge->
sh[5] != NULL) {
4120 fapex =
sapex(*iedge);
4128 spivot(*iedge, parentsh);
4129 if (parentsh.
sh != NULL) {
4141 if (colseg.
sh != iedge->
sh) {
4143 spivot(colseg, parentsh);
4144 printf(
"PLC Error: Two segments are overlapping.\n");
4145 printf(
" Segment 1: [%d, %d] #%d (%d)\n",
pointmark(
sorg(colseg)),
4148 printf(
" Segment 2: [%d, %d] #%d (%d)\n",
pointmark(forg),
4165 }
else if (etype == 2) {
4166 printf(
"PLC Error: A segment lies in a facet.\n");
4169 printf(
" Facet: [%d,%d,%d] #%d\n",
pointmark(forg),
4188 if (colface.
sh != iedge->
sh) {
4189 printf(
"PLC Error: Two facets are overlapping.\n");
4190 printf(
" Facet 1: [%d,%d,%d] #%d\n",
pointmark(forg),
4192 printf(
" Facet 2: [%d,%d,%d] #%d\n",
pointmark(
sorg(colface)),
4218 printf(
"PLC Error: A vertex lies in a segment.\n");
4219 printf(
" Vertex: [%d] (%g,%g,%g).\n",
pointmark(pp),pp[0],pp[1],pp[2]);
4220 printf(
" Segment: [%d, %d] #%d (%d)\n",
pointmark(forg),
4236 }
else if (etype == 2) {
4237 printf(
"PLC Error: A vertex lies in a facet.\n");
4238 printf(
" Vertex: [%d] (%g,%g,%g).\n",
pointmark(pp),pp[0],pp[1],pp[2]);
4257 face parentseg, parentsh;
4259 spivot(parentseg, parentsh);
4260 if (parentseg.
sh != NULL) {
4264 printf(
"PLC Error: Two segments intersect at point (%g,%g,%g).\n",
4265 pp[0], pp[1], pp[2]);
4266 printf(
" Segment 1: [%d, %d], #%d (%d)\n",
pointmark(forg),
4268 printf(
" Segment 2: [%d, %d], #%d (%d)\n",
pointmark(p1),
4285 }
else if (etype == 2) {
4286 printf(
"PLC Error: A segment and a facet intersect at point");
4287 printf(
" (%g,%g,%g).\n", pp[0], pp[1], pp[2]);
4288 printf(
" Segment: [%d, %d], #%d (%d)\n",
pointmark(p1),
4291 printf(
" Facet: [%d,%d,%d] #%d\n",
pointmark(forg),
4314 if (parentsh.
sh != NULL) {
4319 printf(
"PLC Error: A segment and a facet intersect at point");
4320 printf(
" (%g,%g,%g).\n", pp[0], pp[1], pp[2]);
4321 printf(
" Segment : [%d, %d], #%d (%d)\n",
pointmark(forg),
4323 printf(
" Facet : [%d, %d, %d] #%d.\n",
pointmark(p1),
4339 }
else if (etype == 2) {
4340 printf(
"PLC Error: Two facets intersect at point (%g,%g,%g).\n",
4341 pp[0], pp[1], pp[2]);
4342 printf(
" Facet 1: [%d, %d, %d] #%d.\n",
pointmark(forg),
4344 printf(
" Facet 2: [%d, %d, %d] #%d.\n",
pointmark(p1),
4378 spivot(checkseg, parentsh);
4382 REAL P[3], Q[3], tp = 0, tq = 0;
4385 printf(
"PLC Error: Two segments intersect at point (%g,%g,%g).\n",
4387 printf(
" Segment 1: [%d, %d] #%d (%d)\n",
pointmark(forg),
4389 printf(
" Segment 2: [%d, %d] #%d (%d)\n",
pointmark(p1),
4406 }
else if (etype == 2) {
4407 printf(
"PLC Error: A segment and a facet intersect at point");
4408 printf(
" (%g,%g,%g).\n", P[0], P[1], P[2]);
4409 printf(
" Segment: [%d, %d] #%d (%d)\n",
pointmark(p1),
4412 printf(
" Facet: [%d, %d, %d] #%d.\n",
pointmark(forg),
4441 printf(
"PLC Error: A segment and a facet intersect at point");
4442 printf(
" (%g,%g,%g).\n", ip[0], ip[1], ip[2]);
4443 printf(
" Segment: [%d, %d] #%d (%d)\n",
pointmark(forg),
4445 printf(
" Facet: [%d, %d, %d] #%d.\n",
pointmark(p1),
4461 }
else if (etype == 2) {
4462 printf(
"PLC Error: Two facets intersect at point (%g,%g,%g).\n",
4463 ip[0], ip[1], ip[2]);
4464 printf(
" Facet 1: [%d, %d, %d] #%d.\n",
pointmark(forg),
4466 printf(
" Facet 2: [%d, %d, %d] #%d.\n",
pointmark(p1),
4507 triface* iedge,
int intflag,
int* types,
int* poss)
4510 point e1 = NULL, e2 = NULL, e3 = NULL;
4511 int etype = 0, geomtag = 0, facemark = 0;
4542 printf(
"PLC Error: A segment and a facet intersect at point");
4543 printf(
" (%g,%g,%g).\n", ip[0], ip[1], ip[2]);
4546 printf(
" Facet: [%d,%d,%d] #%d\n",
pointmark(p1),
4563 printf(
"PLC Error: Two facets intersect at point");
4564 printf(
" (%g,%g,%g).\n", ip[0], ip[1], ip[2]);
4567 printf(
" Facet 2: [%d,%d,%d] #%d\n",
pointmark(p1),
4586 point crosspt = NULL;
4589 }
else if (poss[0] == 1) {
4591 }
else if (poss[0] == 2) {
4598 printf(
"PLC Error: A vertex and a segment intersect at (%g,%g,%g)\n",
4599 crosspt[0], crosspt[1], crosspt[2]);
4600 printf(
" Vertex: #%d\n",
pointmark(crosspt));
4618 printf(
"PLC Error: A vertex and a facet intersect at (%g,%g,%g)\n",
4619 crosspt[0], crosspt[1], crosspt[2]);
4620 printf(
" Vertex: #%d\n",
pointmark(crosspt));
4621 printf(
" Facet: [%d,%d,%d] #%d\n",
pointmark(p1),
4642 }
else if ((types[0] == (
int)
TOUCHFACE) ||
4645 point touchpt = NULL;
4647 touchpt =
org(*iedge);
4648 }
else if (poss[1] == 1) {
4649 touchpt =
dest(*iedge);
4654 printf(
"PLC Error: A vertex and a facet intersect at (%g,%g,%g)\n",
4655 touchpt[0], touchpt[1], touchpt[2]);
4656 printf(
" Vertex: #%d\n",
pointmark(touchpt));
4657 printf(
" Facet: [%d,%d,%d] #%d\n",
pointmark(p1),
4677 }
else if (types[0] == (
int)
SHAREVERT) {
4682 }
else if (intflag == 4) {
4684 printf(
"PLC Error: Two facets are overlapping.\n");
4685 printf(
" Facet 1: [%d,%d,%d] #%d\n",
pointmark(e1),
4687 printf(
" Facet 2: [%d,%d,%d] #%d\n",
pointmark(p1),
4746 triface topcastets[3], botcastets[3];
4748 point pa, pb, pc, pd, pe;
4749 REAL attrib, volume;
4757 newface = fliptets[0];
4758 fliptets[0] = fliptets[1];
4759 fliptets[1] = newface;
4777 pa =
org(fliptets[0]);
4778 pb =
dest(fliptets[0]);
4779 pc =
apex(fliptets[0]);
4780 pd =
oppo(fliptets[0]);
4781 pe =
oppo(fliptets[1]);
4786 for (i = 0; i < 3; i++) {
4787 fnext(fliptets[0], topcastets[i]);
4790 for (i = 0; i < 3; i++) {
4791 fnext(fliptets[1], botcastets[i]);
4796 fliptets[0].
ver = 11;
4797 fliptets[1].
ver = 11;
4803 if (fliptets[0].tet[8] != NULL) {
4805 fliptets[0].
tet[8] = NULL;
4807 if (fliptets[1].tet[8] != NULL) {
4809 fliptets[1].
tet[8] = NULL;
4814 if (fliptets[0].tet[9] != NULL) {
4816 fliptets[0].
tet[9] = NULL;
4818 if (fliptets[1].tet[9] != NULL) {
4820 fliptets[1].
tet[9] = NULL;
4854 for (i = 0; i < 3; i++) {
4867 if (
fc->remove_ndelaunay_edge) {
4868 REAL volneg[2], volpos[3], vol_diff;
4890 vol_diff = volpos[0] + volpos[1] + volpos[2] - volneg[0] - volneg[1];
4891 fc->tetprism_vol_sum += vol_diff;
4895 for (i = 0; i < 3; i++) {
4896 esym(fliptets[i], newface);
4897 bond(newface, fliptets[(i + 1) % 3]);
4900 for (i = 0; i < 3; i++) {
4902 bond(newface, topcastets[i]);
4905 for (i = 0; i < 3; i++) {
4907 bond(newface, botcastets[i]);
4915 for (i = 0; i < 3; i++) {
4921 if (
fc->chkencflag & 1) {
4927 for (i = 0; i < 3; i++) {
4928 eprev(topcastets[i], casface);
4931 enext(fliptets[i], newface);
4934 esym(fliptets[(i + 2) % 3], newface);
4938 if (
fc->chkencflag & 1) {
4944 for (i = 0; i < 3; i++) {
4945 enext(botcastets[i], casface);
4948 eprev(fliptets[i], newface);
4951 esym(fliptets[(i + 2) % 3], newface);
4955 if (
fc->chkencflag & 1) {
4965 for (i = 0; i < 3; i++) {
4967 tspivot(topcastets[i], checksh);
4970 tsbond(newface, checksh);
4971 if (
fc->chkencflag & 2) {
4976 for (i = 0; i < 3; i++) {
4978 tspivot(botcastets[i], checksh);
4981 tsbond(newface, checksh);
4982 if (
fc->chkencflag & 2) {
4989 if (
fc->chkencflag & 4) {
4991 for (i = 0; i < 3; i++) {
5004 if (dummyflag != 0) {
5006 if (dummyflag == -1) {
5008 for (i = 0; i < 3; i++) {
5012 newface = fliptets[1];
5013 fliptets[1] = fliptets[2];
5014 fliptets[2] = newface;
5017 if (dummyflag == 1) {
5019 newface = fliptets[0];
5020 fliptets[0] = fliptets[2];
5021 fliptets[2] = fliptets[1];
5022 fliptets[1] = newface;
5025 newface = fliptets[0];
5026 fliptets[0] = fliptets[1];
5027 fliptets[1] = fliptets[2];
5028 fliptets[2] = newface;
5034 if (
fc->enqflag > 0) {
5036 for (i = 0; i < 3; i++) {
5040 if (
fc->enqflag > 1) {
5041 for (i = 0; i < 3; i++) {
5085 triface topcastets[3], botcastets[3];
5089 point pa, pb, pc, pd, pe;
5090 REAL attrib, volume;
5092 int spivot = -1, scount = 0;
5100 for (i = 0; i < 3; i++) {
5104 newface = fliptets[1];
5105 fliptets[1] = fliptets[2];
5106 fliptets[2] = newface;
5112 newface = fliptets[0];
5113 fliptets[0] = fliptets[1];
5114 fliptets[1] = fliptets[2];
5115 fliptets[2] = newface;
5118 newface = fliptets[0];
5119 fliptets[0] = fliptets[2];
5120 fliptets[2] = fliptets[1];
5121 fliptets[1] = newface;
5128 pa =
apex(fliptets[0]);
5129 pb =
apex(fliptets[1]);
5130 pc =
apex(fliptets[2]);
5131 pd =
dest(fliptets[0]);
5132 pe =
org(fliptets[0]);
5137 for (i = 0; i < 3; i++) {
5139 fsym(casface, topcastets[i]);
5141 for (i = 0; i < 3; i++) {
5143 fsym(casface, botcastets[i]);
5148 for (i = 0; i < 3; i++) {
5149 tspivot(fliptets[i], flipshs[i]);
5150 if (flipshs[i].sh != NULL) {
5161 fliptets[0].
ver = 11;
5162 fliptets[1].
ver = 11;
5167 if (fliptets[0].tet[8] != NULL) {
5169 fliptets[0].
tet[8] = NULL;
5171 if (fliptets[1].tet[8] != NULL) {
5173 fliptets[1].
tet[8] = NULL;
5178 if (fliptets[0].tet[9] != NULL) {
5180 fliptets[0].
tet[9] = NULL;
5182 if (fliptets[1].tet[9] != NULL) {
5184 fliptets[1].
tet[9] = NULL;
5197 for (j = 0; j < 2; j++) {
5239 if (
fc->remove_ndelaunay_edge) {
5240 REAL volneg[3], volpos[2], vol_diff;
5262 vol_diff = volpos[0] + volpos[1] - volneg[0] - volneg[1] - volneg[2];
5263 fc->tetprism_vol_sum += vol_diff;
5267 bond(fliptets[0], fliptets[1]);
5269 for (i = 0; i < 3; i++) {
5270 esym(fliptets[0], newface);
5271 bond(newface, topcastets[i]);
5275 for (i = 0; i < 3; i++) {
5276 esym(fliptets[1], newface);
5277 bond(newface, botcastets[i]);
5283 for (i = 0; i < 3; i++) {
5290 if (
fc->chkencflag & 1) {
5298 for (i = 0; i < 3; i++) {
5299 esym(fliptets[0], newface);
5301 enext(topcastets[i], casface);
5306 if (
fc->chkencflag & 1) {
5313 for (i = 0; i < 3; i++) {
5314 esym(fliptets[1], newface);
5316 eprev(botcastets[i], casface);
5321 if (
fc->chkencflag & 1) {
5332 for (i = 0; i < 3; i++) {
5334 tspivot(topcastets[i], checksh);
5335 esym(fliptets[0], newface);
5337 tsbond(newface, checksh);
5338 if (
fc->chkencflag & 2) {
5345 for (i = 0; i < 3; i++) {
5347 tspivot(botcastets[i], checksh);
5348 esym(fliptets[1], newface);
5350 tsbond(newface, checksh);
5351 if (
fc->chkencflag & 2) {
5361 flipfaces[0] = flipshs[(
spivot + 1) % 3];
5362 flipfaces[1] = flipshs[(
spivot + 2) % 3];
5364 flip22(flipfaces, 0,
fc->chkencflag);
5368 topcastets[0] = fliptets[0];
5369 botcastets[0] = fliptets[1];
5370 for (i = 0; i < ((
spivot + 1) % 3); i++) {
5378 tspivot(topcastets[0], checksh);
5379 if (checksh.
sh == NULL) {
5380 tsbond(topcastets[0], flipfaces[0]);
5383 tsbond(topcastets[0], flipfaces[0]);
5392 tspivot(botcastets[0], checksh);
5393 if (checksh.
sh == NULL) {
5394 tsbond(botcastets[0], flipfaces[1]);
5397 tsbond(botcastets[0], flipfaces[1]);
5405 if (
fc->chkencflag & 4) {
5407 for (i = 0; i < 2; i++) {
5419 if (dummyflag != 0) {
5421 if (dummyflag == -1) {
5423 newface = fliptets[0];
5424 fliptets[0] = fliptets[1];
5425 fliptets[1] = newface;
5428 if (dummyflag == 1) {
5439 if (
fc->enqflag > 0) {
5446 if (
fc->enqflag > 1) {
5453 esym(fliptets[0], newface);
5455 esym(fliptets[1], newface);
5487 triface topcastets[3], botcastet;
5490 point pa, pb, pc, pd, pp;
5492 int spivot = -1, scount = 0;
5496 pa =
org(fliptets[3]);
5497 pb =
dest(fliptets[3]);
5498 pc =
apex(fliptets[3]);
5499 pd =
dest(fliptets[0]);
5500 pp =
org(fliptets[0]);
5505 for (i = 0; i < 3; i++) {
5506 enext(fliptets[i], topcastets[i]);
5510 fsym(fliptets[3], botcastet);
5515 for (i = 0; i < 3; i++) {
5516 fnext(fliptets[3], newface);
5518 if (flipshs[i].sh != NULL) {
5531 for (i = 0; i < 3; i++) {
5532 esym(neightet, newface);
5544 fliptets[0].
ver = 11;
5549 if (fliptets[0].tet[8] != NULL) {
5551 fliptets[0].
tet[8] = NULL;
5556 if (fliptets[0].tet[9] != NULL) {
5558 fliptets[0].
tet[9] = NULL;
5562 for (i = 1; i < 4; i++) {
5601 if (dummyflag > 0) {
5604 }
else if (dummyflag < 0) {
5613 if (
fc->remove_ndelaunay_edge) {
5614 REAL volneg[4], volpos[1], vol_diff;
5615 if (dummyflag > 0) {
5638 }
else if (dummyflag < 0) {
5651 vol_diff = volpos[0] - volneg[0] - volneg[1] - volneg[2] - volneg[3];
5652 fc->tetprism_vol_sum += vol_diff;
5656 for (i = 0; i < 3; i++) {
5657 esym(fliptets[0], newface);
5658 bond(newface, topcastets[i]);
5661 bond(fliptets[0], botcastet);
5666 for (i = 0; i < 3; i++) {
5667 eprev(topcastets[i], newface);
5670 esym(fliptets[0], newface);
5674 if (
fc->chkencflag & 1) {
5680 for (i = 0; i < 3; i++) {
5685 if (
fc->chkencflag & 1) {
5696 for (i = 0; i < 3; i++) {
5698 tspivot(topcastets[i], checksh);
5699 esym(fliptets[0], newface);
5701 tsbond(newface, checksh);
5702 if (
fc->chkencflag & 2) {
5711 tsbond(fliptets[0], checksh);
5712 if (
fc->chkencflag & 2) {
5726 for (i = 0; i < 3; i++) {
5731 for (i = 0; i < 3; i++) {
5738 sesym(flipshs[3], checksh);
5739 tsbond(newface, checksh);
5742 tsbond(fliptets[0], flipshs[3]);
5743 fsym(fliptets[0], newface);
5744 sesym(flipshs[3], checksh);
5745 tsbond(newface, checksh);
5750 if (
fc->chkencflag & 4) {
5760 if (
fc->enqflag > 0) {
5763 if (
fc->enqflag > 1) {
5764 for (i = 0; i < 3; i++) {
5765 esym(fliptets[0], newface);
5807 triface fliptets[3], spintet, flipedge;
5809 point pa, pb, pc, pd, pe, pf;
5811 int hullflag, hulledgeflag;
5812 int reducflag, rejflag;
5813 int reflexlinkedgecount;
5819 pa =
org(abtets[0]);
5820 pb =
dest(abtets[0]);
5824 reflexlinkedgecount = 0;
5826 for (i = 0; i < n; i++) {
5839 pc =
apex(abtets[i]);
5840 pd =
apex(abtets[(i + 1) % n]);
5841 pe =
apex(abtets[(i - 1 + n) % n]);
5852 if (hullflag == 0) {
5862 }
else if (ori == 0) {
5868 pf =
apex(abtets[(i + 2) % n]);
5875 reflexlinkedgecount++;
5886 pf =
apex(abtets[(i + 2) % n]);
5903 if (
getedge(pe, pd, &spintet)) {
5913 if (
fc->checkflipeligibility) {
5920 fliptets[0] = abtets[i];
5921 fsym(fliptets[0], fliptets[1]);
5943 abtets[(i - 1 + n) % n] = fliptets[0];
5944 for (j = i; j < n - 1; j++) {
5945 abtets[j] = abtets[j + 1];
5953 abtets[n - 1].
ver = 0;
5956 abtets[n - 1].
ver |= (1 << 4);
5958 abtets[n - 1].
ver |= (i << 6);
5960 if (
fc->collectnewtets) {
5963 for (j = 1; j < 3; j++) {
5965 *parytet = fliptets[j];
5970 nn =
flipnm(abtets, n - 1, level, abedgepivot,
fc);
5977 if (
fc->unflip || (ori == 0)) {
5984 fliptets[0] = abtets[(i-1 + (n-1)) % (n-1)];
5986 fnext(fliptets[0], fliptets[1]);
5987 fnext(fliptets[1], fliptets[2]);
5991 for (j = 0; j < 2; j++) {
5995 for (j = n - 2; j>= i; j--) {
5996 abtets[j + 1] = abtets[j];
6001 esym(fliptets[1], abtets[(i - 1 + n) % n]);
6002 abtets[i] = fliptets[0];
6004 if (
fc->collectnewtets) {
6023 if (reflexlinkedgecount > 0) {
6028 for (i = 0; i < n; i++) {
6034 pc =
apex(abtets[i]);
6038 pd =
apex(abtets[(i + 1) % n]);
6039 pe =
apex(abtets[(i - 1 + n) % n]);
6051 enext(abtets[i], flipedge);
6059 eprev(abtets[i], flipedge);
6064 if (!edgepivot)
continue;
6070 if (
fc->collectencsegflag) {
6071 face checkseg, *paryseg;
6077 *paryseg = checkseg;
6094 if (spintet.
tet == flipedge.
tet)
break;
6117 tmpabtets[j] = spintet;
6122 if (spintet.
tet == flipedge.
tet)
break;
6126 nn =
flipnm(tmpabtets, n1, level + 1, edgepivot,
fc);
6131 if (edgepivot == 1) {
6133 spintet = tmpabtets[0];
6139 spintet = tmpabtets[1];
6146 abtets[(i - 1 + n) % n] = spintet;
6147 for (j = i; j < n - 1; j++) {
6148 abtets[j] = abtets[j + 1];
6153 abtets[n - 1].
ver = 0;
6155 abtets[n - 1].
ver |= edgepivot;
6157 abtets[n - 1].
ver |= (1 << 5);
6159 abtets[n - 1].
ver |= (i << 6);
6161 abtets[n - 1].
ver |= (n1 << 19);
6166 tmpabtets[0].
ver = (1 << 5);
6169 nn =
flipnm(abtets, n - 1, level, abedgepivot,
fc);
6182 if (edgepivot == 1) {
6184 tmpabtets[0] = abtets[((i-1)+(n-1))%(n-1)];
6188 fsym(tmpabtets[0], tmpabtets[1]);
6191 tmpabtets[1] = abtets[((i-1)+(n-1))%(n-1)];
6195 fsym(tmpabtets[1], tmpabtets[0]);
6202 for (j = n - 2; j >= i; j--) {
6203 abtets[j + 1] = abtets[j];
6205 if (edgepivot == 1) {
6209 fliptets[0] = tmpabtets[1];
6212 fliptets[1] = tmpabtets[0];
6219 fliptets[0] = tmpabtets[1];
6222 fliptets[1] = tmpabtets[0];
6226 for (j = 0; j < 2; j++) {
6230 abtets[(i - 1 + n) % n] = fliptets[0];
6231 abtets[i] = fliptets[1];
6234 delete [] tmpabtets;
6252 for (j = 0; j < nn; j++) {
6256 delete [] tmpabtets;
6268 pc =
apex(abtets[1]);
6269 pd =
apex(abtets[2]);
6270 pe =
apex(abtets[0]);
6273 pc =
apex(abtets[2]);
6274 pd =
apex(abtets[0]);
6275 pe =
apex(abtets[1]);
6278 pc =
apex(abtets[0]);
6279 pd =
apex(abtets[1]);
6280 pe =
apex(abtets[2]);
6288 if (hullflag == 0) {
6311 if (reducflag == 1) {
6314 point searchpt = NULL, chkpt;
6315 REAL bigvol = 0.0, ori1, ori2;
6320 fliptets[0] = abtets[hullflag % 3];
6322 spintet = fliptets[0];
6325 chkpt =
oppo(spintet);
6326 if (chkpt == pb)
break;
6335 if (searchpt != NULL) {
6337 ori1 =
orient3d(pd, pc, searchpt, pa);
6338 ori2 =
orient3d(pd, pc, searchpt, pb);
6339 if (ori1 * ori2 >= 0.0) {
6342 ori1 =
orient3d(pa, pb, searchpt, pc);
6343 ori2 =
orient3d(pa, pb, searchpt, pd);
6344 if (ori1 * ori2 >= 0.0) {
6364 for (j = 0; j < 3; j++) {
6376 }
else if (nn == 2) {
6379 eorgoppo(abtets[(edgepivot + 1) % 3], spintet);
6388 }
else if (nn == 3) {
6393 if (!rejflag &&
fc->checkflipeligibility) {
6404 if (
fc->remove_ndelaunay_edge) {
6408 if ((
fc->tetprism_vol_sum >= 0.0) ||
6409 (fabs(
fc->tetprism_vol_sum) <
fc->bak_tetprism_vol)) {
6413 for (j = 0; j < 3; j++) {
6420 if (
fc->collectnewtets) {
6424 for (j = 0; j < 2; j++) {
6426 *parytet = abtets[j];
6432 if (abedgepivot == 1) {
6433 *parytet = abtets[1];
6435 *parytet = abtets[0];
6481 triface fliptets[3], flipface;
6496 if (
fc->collectnewtets) {
6498 if (abedgepivot == 0) {
6512 for (i = nn; i < n; i++) {
6516 fliptype = ((abtets[i].
ver >> 4) & 3);
6517 if (fliptype == 1) {
6519 t = (abtets[i].
ver >> 6);
6522 printf(
" Recover a 2-to-3 flip at f[%d].\n", t);
6526 fliptets[0] = abtets[((t - 1) + i) % i];
6530 fnext(fliptets[0], fliptets[1]);
6531 fnext(fliptets[1], fliptets[2]);
6537 for (j = i - 1; j >= t; j--) {
6538 abtets[j + 1] = abtets[j];
6543 esym(fliptets[1], abtets[((t-1) + (i+1)) % (i+1)]);
6544 abtets[t] = fliptets[0];
6545 if (
fc->collectnewtets) {
6550 }
else if (fliptype == 2) {
6551 tmpabtets = (
triface *) (abtets[i].tet);
6552 n1 = ((abtets[i].
ver >> 19) & 8191);
6553 edgepivot = (abtets[i].
ver & 3);
6554 t = ((abtets[i].
ver >> 6) & 8191);
6557 printf(
" Recover a %d-to-m flip at e[%d] of f[%d].\n", n1,
6564 if (edgepivot == 1) {
6566 tmpabtets[0] = abtets[(t - 1 + i) % i];
6570 fsym(tmpabtets[0], tmpabtets[1]);
6573 tmpabtets[1] = abtets[(t - 1 + i) % i];
6577 fsym(tmpabtets[1], tmpabtets[0]);
6584 for (j = i - 1; j >= t; j--) {
6585 abtets[j + 1] = abtets[j];
6587 if (edgepivot == 1) {
6591 fliptets[0] = tmpabtets[1];
6594 fliptets[1] = tmpabtets[0];
6601 fliptets[0] = tmpabtets[1];
6604 fliptets[1] = tmpabtets[0];
6609 abtets[((t-1) + (i+1)) % (i+1)] = fliptets[0];
6610 abtets[t] = fliptets[1];
6617 printf(
" Release %d spaces at f[%d].\n", n1, i);
6619 delete [] tmpabtets;
6645 triface *cavetet, spintet, neightet, neineitet, *parytet;
6646 triface oldtet, newtet, newneitet;
6647 face checksh, neighsh, *parysh;
6648 face checkseg, *paryseg;
6649 point *pts, pa, pb, pc, *parypt;
6652 REAL attrib, volume;
6658 printf(
" Insert point %d\n",
pointmark(insertpt));
6662 if (searchtet->
tet != NULL) {
6667 if (searchtet->
tet == NULL) {
6676 loc =
locate(insertpt, searchtet);
6679 ivf->
iloc = (int) loc;
6685 sign =
orient4d_s(pts[4], pts[5], pts[6], pts[7], insertpt,
6686 pts[4][3], pts[5][3], pts[6][3], pts[7][3],
6705 for (i = 0; i < 4; i++) {
6709 *parytet = neightet;
6713 *parytet = *searchtet;
6717 for (i = 0; i < 4; i++) {
6721 *parytet = neightet;
6725 *parytet = *searchtet;
6726 }
else if (loc ==
ONFACE) {
6729 j = (searchtet->
ver & 3);
6730 for (i = 1; i < 4; i++) {
6731 decode(searchtet->
tet[(j + i) % 4], neightet);
6734 *parytet = neightet;
6737 j = (spintet.
ver & 3);
6738 for (i = 1; i < 4; i++) {
6739 decode(spintet.
tet[(j + i) % 4], neightet);
6742 *parytet = neightet;
6749 *parytet = *searchtet;
6752 if ((splitsh != NULL) && (splitsh->
sh != NULL)) {
6759 }
else if (loc ==
ONEDGE) {
6762 spintet = *searchtet;
6768 *parytet = neightet;
6773 *parytet = neightet;
6778 if (spintet.
tet == searchtet->
tet)
break;
6783 if ((splitseg != NULL) && (splitseg->
sh != NULL)) {
6785 splitseg->
shver = 0;
6786 spivot(*splitseg, *splitsh);
6788 if (splitsh != NULL) {
6789 if (splitsh->
sh != NULL) {
6791 pa =
sorg(*splitsh);
6795 if (
sorg(neighsh) != pa) {
6808 if (neighsh.
sh == splitsh->
sh)
break;
6809 if (neighsh.
sh == NULL)
break;
6814 }
else if (loc ==
INSTAR) {
6821 for (j = 0; j < 4; j++) {
6827 *parytet = neightet;
6857 sign =
orient4d_s(pts[4], pts[5], pts[6], pts[7], insertpt,
6858 pts[4][3], pts[5][3], pts[6][3], pts[7][3],
6861 sign =
insphere_s(pts[4], pts[5], pts[6], pts[7], insertpt);
6863 enqflag = (sign < 0.0);
6867 ori =
orient3d(pts[4], pts[5], pts[6], insertpt);
6872 }
else if (ori == 0.0) {
6882 sign =
orient4d_s(pts[4],pts[5],pts[6],pts[7], insertpt,
6883 pts[4][3], pts[5][3], pts[6][3],
6884 pts[7][3], insertpt[3]);
6886 sign =
insphere_s(pts[4],pts[5],pts[6],pts[7], insertpt);
6888 enqflag = (sign < 0.0);
6908 sign =
orient4d_s(pts[4],pts[5],pts[6],pts[7], insertpt,
6909 pts[4][3], pts[5][3], pts[6][3],
6910 pts[7][3], insertpt[3]);
6912 sign =
insphere_s(pts[4],pts[5],pts[6],pts[7], insertpt);
6914 enqflag = (sign < 0.0);
6928 k = (cavetet->
ver & 3);
6929 for (j = 1; j < 4; j++) {
6930 decode(cavetet->
tet[(j + k) % 4], neightet);
6932 *parytet = neightet;
6936 *parytet = *cavetet;
6941 *parytet = *cavetet;
6955 for (j = 0; j < 6; j++) {
6961 *paryseg = checkseg;
6981 for (j = 0; j < 4; j++) {
7016 for (j = 0; j < 3; j++) {
7018 spivot(checksh, neighsh);
7029 sapex(neighsh), insertpt);
7086 *parytet = neightet;
7088 for (j = 0; j < 3; j++) {
7089 esym(neightet, neineitet);
7092 *parytet = neineitet;
7113 if (spintet.
tet == neightet.
tet)
break;
7119 pb =
dest(neightet);
7126 ori =
orient3d(pa, pb, pc, insertpt);
7129 esym(spintet, neineitet);
7130 pc =
apex(neineitet);
7132 ori =
orient3d(pb, pa, pc, insertpt);
7142 if (spintet.
tet == neightet.
tet)
break;
7150 printf(
" Cut tet (%d, %d, %d, %d)\n",
7159 *parytet = neightet;
7161 for (j = 0; j < 3; j++) {
7162 esym(neightet, neineitet);
7165 *parytet = neineitet;
7178 fsym(*cavetet, neightet);
7188 pb =
dest(*cavetet);
7189 pc =
apex(*cavetet);
7192 if ((fabs(volume) / attrib) <
b->
epsilon) {
7195 ori =
orient3d(pa, pb, pc, insertpt);
7201 enqflag = (ori > 0);
7215 *parytet = *cavetet;
7221 for (j = 0; j < 3; j++) {
7222 esym(neightet, neineitet);
7225 *parytet = neineitet;
7244 fsym(*cavetet, neightet);
7248 *parytet = *cavetet;
7261 *parytet = *cavetet;
7303 if (cutshcount > 0) {
7307 if ((splitsh != NULL) && (splitsh->
sh != NULL)) {
7311 }
else if (loc ==
ONEDGE) {
7312 if ((splitseg != NULL) && (splitseg->
sh != NULL)) {
7316 if ((splitsh != NULL) && (splitsh->
sh != NULL)) {
7318 pa =
sorg(*splitsh);
7322 if (
sorg(neighsh) != pa) {
7330 if (neighsh.
sh == splitsh->
sh)
break;
7331 if (neighsh.
sh == NULL)
break;
7366 pts = (
point *) &(searchtet->
tet[4]);
7367 for (i = 0; i < 3; i++) {
7372 pts = (
point *) &(searchtet->
tet[4]);
7373 for (i = 0; i < 4; i++) {
7377 }
else if (loc ==
ONFACE) {
7378 pts = (
point *) &(searchtet->
tet[4]);
7379 for (i = 0; i < 3; i++) {
7387 fsym(*searchtet, spintet);
7390 *parypt =
oppo(spintet);
7392 }
else if (loc ==
ONEDGE) {
7393 spintet = *searchtet;
7395 *parypt =
org(spintet);
7397 *parypt =
dest(spintet);
7401 *parypt =
apex(spintet);
7404 if (spintet.
tet == searchtet->
tet)
break;
7408 int rejptflag = (ivf->
rejflag & 4);
7440 printf(
"Warning: Two points, %d and %d, are very close.\n",
7442 printf(
" Creating a very short edge (len = %g) (< %g).\n",
7444 printf(
" You may try a smaller tolerance (-T) (current is %g)\n",
7446 printf(
" to avoid this warning.\n");
7451 ivf->
iloc = (int) loc;
7458 ivf->
iloc = (int) loc;
7471 pts = (
point *) &(cavetet->
tet[4]);
7472 for (j = 0; j < 4; j++) {
7496 if (len < ivf->smlen) {
7545 neineitet = spintet;
7550 if (spintet.
tet == neightet.
tet)
break;
7558 *paryseg = checkseg;
7567 checkseg = *paryseg;
7570 *paryseg = checkseg;
7592 for (j = 0; j < 2; j++) {
7608 }
else if (k == 1) {
7635 neightet = *cavetet;
7638 fsym(neightet, oldtet);
7667 bond(newtet, neightet);
7685 fsym(oldtet, neightet);
7686 fsym(neightet, newtet);
7689 for (j = 0; j < 3; j++) {
7690 esym(newtet, neightet);
7691 if (neightet.
tet[neightet.
ver & 3] == NULL) {
7698 fsym(spintet, newneitet);
7700 bond(neightet, newneitet);
7703 *parytet = neightet;
7721 fsym(neightet, spintet);
7723 tsbond(spintet, *parysh);
7739 if (spintet.
tet == neightet.
tet)
break;
7745 if (((splitsh != NULL) && (splitsh->
sh != NULL)) ||
7746 ((splitseg != NULL) && (splitseg->
sh != NULL))) {
7757 spivot(*parysh, checksh);
7759 if (checksh.
sh[3] != NULL) {
7770 fsym(spintet, neightet);
7775 if (
apex(spintet) == insertpt)
break;
7778 if (
sorg(checksh) !=
org(spintet)) {
7782 tsbond(spintet, checksh);
7785 tsbond(spintet, checksh);
7794 spivot(*parysh, checksh);
7796 if (checksh.
sh[3] != NULL) {
7819 if (splitseg != NULL) {
7824 checkseg = *paryseg;
7827 spivot(checkseg, checksh);
7828 if (checksh.
sh != NULL) {
7841 if (spintet.
tet == neightet.
tet)
break;
7848 if (splitseg != NULL) {
7852 checkseg = *paryseg;
7858 *paryseg = checkseg;
7864 checkseg = *paryseg;
7871 *paryseg = checkseg;
7889 printf(
" Point #%d is non-regular after the insertion of #%d.\n",
7908 if (splitseg != NULL) {
7930 spivot(*parysh, checksh);
7932 if (checksh.
sh[3] != NULL) {
7957 if (((splitsh != NULL) && (splitsh->
sh != NULL)) ||
7958 ((splitseg != NULL) && (splitseg->
sh != NULL))) {
7967 if (neightet.
tet != NULL) {
7968 if (neightet.
tet[4] != NULL) {
7978 if ((splitseg != NULL) && (splitseg->
sh != NULL)) {
8019 if (((splitsh != NULL) && (splitsh->
sh != NULL)) ||
8020 ((splitseg != NULL) && (splitseg->
sh != NULL))) {
8058 if ((splitseg != NULL) && (splitseg->
sh != NULL)) {
8097 int gc[8], N, mask, travel_bit;
8102 N = (n == 2) ? 4 : 8;
8103 mask = (n == 2) ? 3 : 7;
8106 for (i = 0; i < N; i++) {
8107 gc[i] = i ^ (i >> 1);
8110 for (e = 0; e < N; e++) {
8111 for (d = 0; d < n; d++) {
8116 for (i = 0; i < N; i++) {
8118 k = gc[i] * (travel_bit * 2);
8119 g = ((k | (k / N)) & mask);
8128 for (i = 1; i < N; i++) {
8130 v = (v ^ (v - 1)) >> 1;
8131 for (
c = 0; v;
c++) {
8156 axis = (gc0 ^ gc1) >> 1;
8160 split = 0.5 * (bxmin + bxmax);
8161 }
else if (axis == 1) {
8162 split = 0.5 * (bymin + bymax);
8164 split = 0.5 * (bzmin + bzmax);
8169 d = ((gc0 & (1<<axis)) == 0) ? 1 : -1;
8180 for (; i < arraysize; i++) {
8181 if (vertexarray[i][axis] >= split)
break;
8183 for (; j >= 0; j--) {
8184 if (vertexarray[j][axis] < split)
break;
8187 if (i == (j + 1))
break;
8189 swapvert = vertexarray[i];
8190 vertexarray[i] = vertexarray[j];
8191 vertexarray[j] = swapvert;
8196 for (; i < arraysize; i++) {
8197 if (vertexarray[i][axis] <= split)
break;
8199 for (; j >= 0; j--) {
8200 if (vertexarray[j][axis] > split)
break;
8203 if (i == (j + 1))
break;
8205 swapvert = vertexarray[i];
8206 vertexarray[i] = vertexarray[j];
8207 vertexarray[j] = swapvert;
8219 REAL x1, x2, y1, y2, z1, z2;
8220 int p[9], w, e_w, d_w, k, ei, di;
8221 int n = 3, mask = 7;
8228 bxmin, bxmax, bymin, bymax, bzmin, bzmax);
8230 bxmin, bxmax, bymin, bymax, bzmin, bzmax);
8232 bxmin, bxmax, bymin, bymax, bzmin, bzmax);
8235 bxmin, bxmax, bymin, bymax, bzmin, bzmax) + p[2];
8238 bxmin, bxmax, bymin, bymax, bzmin, bzmax) + p[4];
8241 bxmin, bxmax, bymin, bymax, bzmin, bzmax) + p[4];
8244 bxmin, bxmax, bymin, bymax, bzmin, bzmax) + p[6];
8255 for (w = 0; w < 8; w++) {
8267 k = 2 * ((w - 1) / 2);
8271 e_w = ((k << (d+1)) & mask) | ((k >> (n-d-1)) & mask);
8280 di = (d + d_w + 1) % n;
8283 x1 = 0.5 * (bxmin + bxmax);
8287 x2 = 0.5 * (bxmin + bxmax);
8290 y1 = 0.5 * (bymin + bymax);
8294 y2 = 0.5 * (bymin + bymax);
8297 z1 = 0.5 * (bzmin + bzmax);
8301 z2 = 0.5 * (bzmin + bzmax);
8304 x1, x2, y1, y2, z1, z2, depth+1);
8316 int threshold,
REAL ratio,
int *depth)
8321 if (arraysize >= threshold) {
8323 middle = arraysize * ratio;
8327 hilbert_sort3(&(vertexarray[middle]), arraysize - middle, 0, 0,
8339 unsigned long newrandom;
8341 if (choices >= 714025l) {
8342 newrandom = (
randomseed * 1366l + 150889l) % 714025l;
8343 randomseed = (newrandom * 1366l + 150889l) % 714025l;
8344 newrandom = newrandom * (choices / 714025l) +
randomseed;
8345 if (newrandom >= choices) {
8346 return newrandom - choices;
8373 long sampleblocks, samplesperblock, samplenum;
8374 long tetblocks, i, j;
8375 REAL searchdist, dist;
8378 printf(
" Random sampling tetrahedra for searching point %d.\n",
8383 if (searchtet->
tet == NULL) {
8392 torg =
org(*searchtet);
8393 searchdist = (searchpt[0] - torg[0]) * (searchpt[0] - torg[0]) +
8394 (searchpt[1] - torg[1]) * (searchpt[1] - torg[1]) +
8395 (searchpt[2] - torg[2]) * (searchpt[2] - torg[2]);
8402 dist = (searchpt[0] - torg[0]) * (searchpt[0] - torg[0]) +
8403 (searchpt[1] - torg[1]) * (searchpt[1] - torg[1]) +
8404 (searchpt[2] - torg[2]) * (searchpt[2] - torg[2]);
8405 if (dist < searchdist) {
8425 samplesperblock = 1 + (
samples / tetblocks);
8426 sampleblocks =
samples / samplesperblock;
8428 for (i = 0; i < sampleblocks; i++) {
8429 alignptr = (uintptr_t) (sampleblock + 1);
8433 for (j = 0; j < samplesperblock; j++) {
8434 if (i == tetblocks - 1) {
8443 torg = (
point) tetptr[4];
8444 if (torg != (
point) NULL) {
8445 dist = (searchpt[0] - torg[0]) * (searchpt[0] - torg[0]) +
8446 (searchpt[1] - torg[1]) * (searchpt[1] - torg[1]) +
8447 (searchpt[2] - torg[2]) * (searchpt[2] - torg[2]);
8448 if (dist < searchdist) {
8449 searchtet->
tet = tetptr;
8450 searchtet->
ver = 11;
8455 if (i != tetblocks - 1) j--;
8458 sampleblock = (
void **) *sampleblock;
8487 point torg, tdest, tapex, toppo;
8488 enum {ORGMOVE, DESTMOVE, APEXMOVE} nextmove;
8489 REAL ori, oriorg, oridest, oriapex;
8494 if (searchtet->tet == NULL) {
8507 for (searchtet->ver = 0; searchtet->ver < 4; searchtet->ver++) {
8508 torg =
org(*searchtet);
8509 tdest =
dest(*searchtet);
8510 tapex =
apex(*searchtet);
8511 ori =
orient3d(torg, tdest, tapex, searchpt);
8512 if (ori < 0.0)
break;
8514 if (searchtet->ver == 4) {
8521 toppo =
oppo(*searchtet);
8524 if (toppo == searchpt) {
8533 oriorg =
orient3d(tdest, tapex, toppo, searchpt);
8534 oridest =
orient3d(tapex, torg, toppo, searchpt);
8535 oriapex =
orient3d(torg, tdest, toppo, searchpt);
8546 }
else if (s == 1) {
8547 nextmove = DESTMOVE;
8549 nextmove = APEXMOVE;
8557 nextmove = DESTMOVE;
8567 nextmove = APEXMOVE;
8580 nextmove = DESTMOVE;
8582 nextmove = APEXMOVE;
8586 nextmove = DESTMOVE;
8591 nextmove = APEXMOVE;
8640 if (nextmove == ORGMOVE) {
8642 }
else if (nextmove == DESTMOVE) {
8662 torg =
org(*searchtet);
8663 tdest =
dest(*searchtet);
8664 tapex =
apex(*searchtet);
8684 newflipface->
tt = *flipface;
8688 fstack = newflipface;
8714 triface fliptets[5], *parytet;
8715 point *pts, *parypt, pe;
8729 while (popface != NULL) {
8731 for (i = 4; i < 8; i++) {
8732 if ((pts[i] != newpt) && (pts[i] !=
dummypoint)) {
8749 fliptets[0] = popface->
tt;
8764 fsym(fliptets[0], fliptets[1]);
8765 pts = (
point *) fliptets[1].tet;
8766 ori =
orient3d(pts[4], pts[5], pts[6], newpt);
8771 enext(fliptets[1], fliptets[2]);
8772 eprev(fliptets[1], fliptets[3]);
8775 if (
oppo(fliptets[2]) == newpt) {
8776 if (
oppo(fliptets[3]) == newpt) {
8780 esym(fliptets[2], fliptets[0]);
8781 fnext(fliptets[0], fliptets[1]);
8782 fnext(fliptets[1], fliptets[2]);
8792 if (
oppo(fliptets[3]) == newpt) {
8793 fnext(fliptets[3], fliptets[0]);
8794 fnext(fliptets[0], fliptets[1]);
8795 fnext(fliptets[1], fliptets[2]);
8804 pe =
oppo(fliptets[1]);
8828 fsym(fliptets[0], fliptets[1]);
8840 pts = (
point *) fliptets[1].tet;
8842 sign =
orient4d_s(pts[4], pts[5], pts[6], pts[7], newpt,
8843 pts[4][3], pts[5][3], pts[6][3], pts[7][3],
8846 sign =
insphere_s(pts[4], pts[5], pts[6], pts[7], newpt);
8856 for (i = 0; i < 3; i++) {
8858 if (ori <= 0)
break;
8872 for (i = 0; i < 3; i++) {
8873 fnext(fliptets[i], fliptets[i+1]);
8875 if (fliptets[3].tet == fliptets[0].tet) {
8881 fnext(fliptets[3], fliptets[4]);
8882 if (fliptets[4].tet == fliptets[0].tet) {
8895 fnext(fliptets[3], fliptets[1]);
8896 fnext(fliptets[1], fliptets[2]);
8915 *parytet = fliptets[1];
8951 triface firsttet, tetopa, tetopb, tetopc, tetopd;
8955 printf(
" Create init tet (%d, %d, %d, %d)\n",
pointmark(pa),
8974 bond(firsttet, tetopd);
8975 esym(firsttet, worktet);
8976 bond(worktet, tetopc);
8978 bond(worktet, tetopa);
8980 bond(worktet, tetopb);
8983 esym(tetopc, worktet);
8984 esym(tetopd, worktet1);
8985 bond(worktet, worktet1);
8986 esym(tetopa, worktet);
8988 bond(worktet, worktet1);
8989 esym(tetopb, worktet);
8991 bond(worktet, worktet1);
8994 bond(worktet, worktet1);
8997 bond(worktet, worktet1);
9000 bond(worktet, worktet1);
9036 point *permutarray, swapvertex;
9037 REAL v1[3], v2[3], n[3];
9038 REAL bboxsize, bboxsize2, bboxsize3, ori;
9044 printf(
"Delaunizing vertices...\n");
9048 permutarray =
new point[
in->numberofpoints];
9053 printf(
" Using the input order.\n");
9055 for (i = 0; i <
in->numberofpoints; i++) {
9060 printf(
" Permuting vertices.\n");
9062 srand(
in->numberofpoints);
9063 for (i = 0; i <
in->numberofpoints; i++) {
9064 randindex = rand() % (i + 1);
9065 permutarray[i] = permutarray[randindex];
9070 printf(
" Sorting vertices.\n");
9082 bboxsize2 = bboxsize * bboxsize;
9083 bboxsize3 = bboxsize2 * bboxsize;
9089 if (i ==
in->numberofpoints - 1) {
9090 printf(
"Exception: All vertices are (nearly) identical (Tol = %g).\n",
9097 swapvertex = permutarray[i];
9098 permutarray[i] = permutarray[1];
9099 permutarray[1] = swapvertex;
9106 for (j = 0; j < 3; j++) {
9107 v1[j] = permutarray[1][j] - permutarray[0][j];
9108 v2[j] = permutarray[i][j] - permutarray[0][j];
9111 while ((sqrt(
norm2(n[0], n[1], n[2])) / bboxsize2) <
9114 if (i ==
in->numberofpoints - 1) {
9115 printf(
"Exception: All vertices are (nearly) collinear (Tol = %g).\n",
9119 for (j = 0; j < 3; j++) {
9120 v2[j] = permutarray[i][j] - permutarray[0][j];
9126 swapvertex = permutarray[i];
9127 permutarray[i] = permutarray[2];
9128 permutarray[2] = swapvertex;
9133 ori =
orient3dfast(permutarray[0], permutarray[1], permutarray[2],
9137 if (i ==
in->numberofpoints) {
9138 printf(
"Exception: All vertices are coplanar (Tol = %g).\n",
9142 ori =
orient3dfast(permutarray[0], permutarray[1], permutarray[2],
9147 swapvertex = permutarray[i];
9148 permutarray[i] = permutarray[3];
9149 permutarray[3] = swapvertex;
9156 swapvertex = permutarray[0];
9157 permutarray[0] = permutarray[1];
9158 permutarray[1] = swapvertex;
9166 printf(
" Incrementally inserting vertices.\n");
9182 for (i = 4; i <
in->numberofpoints; i++) {
9191 searchtet.
tet = NULL;
9195 if (
insertpoint(permutarray[i], &searchtet, NULL, NULL, &ivf)) {
9203 swapvertex =
org(searchtet);
9206 printf(
"Warning: Point #%d is coincident with #%d. Ignored!\n",
9214 swapvertex =
org(searchtet);
9216 printf(
"Warning: Point %d is replaced by point %d.\n",
9218 printf(
" Avoid creating a very short edge (len = %g) (< %g).\n",
9220 printf(
" You may try a smaller tolerance (-T) (current is %g)\n",
9222 printf(
" or use the option -M0/1 to avoid such replacement.\n");
9231 printf(
" Point #%d is non-regular, skipped.\n",
9242 delete [] permutarray;
9264 newflipface->
ss = *flipedge;
9265 newflipface->
forg =
sorg(*flipedge);
9283 face bdedges[4], outfaces[4], infaces[4];
9286 point pa, pb, pc, pd;
9289 pa =
sorg(flipfaces[0]);
9290 pb =
sdest(flipfaces[0]);
9291 pc =
sapex(flipfaces[0]);
9292 pd =
sapex(flipfaces[1]);
9294 if (
sorg(flipfaces[1]) != pb) {
9301 senext(flipfaces[0], bdedges[0]);
9302 senext2(flipfaces[0], bdedges[1]);
9303 senext(flipfaces[1], bdedges[2]);
9304 senext2(flipfaces[1], bdedges[3]);
9307 for (i = 0; i < 4; i++) {
9308 spivot(bdedges[i], outfaces[i]);
9309 infaces[i] = outfaces[i];
9310 sspivot(bdedges[i], bdsegs[i]);
9311 if (outfaces[i].sh != NULL) {
9313 spivot(infaces[i], checkface);
9314 while (checkface.
sh != bdedges[i].
sh) {
9315 infaces[i] = checkface;
9316 spivot(infaces[i], checkface);
9346 for (i = 0; i < 4; i++) {
9347 if (outfaces[(3 + i) % 4].sh != NULL) {
9349 if (bdsegs[(3 + i) % 4].sh != NULL) {
9350 bdsegs[(3 + i) % 4].shver = 0;
9351 if (
sorg(bdedges[i]) !=
sorg(bdsegs[(3 + i) % 4])) {
9355 sbond1(bdedges[i], outfaces[(3 + i) % 4]);
9356 sbond1(infaces[(3 + i) % 4], bdedges[i]);
9360 if (bdsegs[(3 + i) % 4].sh != NULL) {
9361 ssbond(bdedges[i], bdsegs[(3 + i) % 4]);
9362 if (chkencflag & 1) {
9371 if (chkencflag & 2) {
9373 for (i = 0; i < 2; i++) {
9382 for (i = 0; i < 4; i++) {
9405 face bdedges[3], outfaces[3], infaces[3];
9411 pa =
sdest(flipfaces[0]);
9412 pb =
sdest(flipfaces[1]);
9413 pc =
sdest(flipfaces[2]);
9418 for (i = 0; i < 3; i++) {
9419 senext(flipfaces[i], bdedges[i]);
9420 spivot(bdedges[i], outfaces[i]);
9421 infaces[i] = outfaces[i];
9422 sspivot(bdedges[i], bdsegs[i]);
9423 if (outfaces[i].sh != NULL) {
9425 spivot(infaces[i], checkface);
9426 while (checkface.
sh != bdedges[i].
sh) {
9427 infaces[i] = checkface;
9428 spivot(infaces[i], checkface);
9458 bdedges[0] = flipfaces[3];
9459 senext(flipfaces[3], bdedges[1]);
9460 senext2(flipfaces[3], bdedges[2]);
9463 for (i = 0; i < 3; i++) {
9464 if (outfaces[i].sh != NULL) {
9466 if (bdsegs[i].sh != NULL) {
9467 bdsegs[i].
shver = 0;
9468 if (
sorg(bdedges[i]) !=
sorg(bdsegs[i])) {
9472 sbond1(bdedges[i], outfaces[i]);
9473 sbond1(infaces[i], bdedges[i]);
9475 if (bdsegs[i].sh != NULL) {
9476 ssbond(bdedges[i], bdsegs[i]);
9484 for (i = 0; i < 3; i++) {
9500 point pa, pb, pc, pd;
9512 flipfaces[0] = popface->
ss;
9514 pb = popface->
fdest;
9519 if (flipfaces[0].sh[3] == NULL)
continue;
9521 if ((
sorg(flipfaces[0]) != pa) || (
sdest(flipfaces[0]) != pb))
continue;
9526 spivot(flipfaces[0], flipfaces[1]);
9527 if (flipfaces[1].sh == NULL)
continue;
9528 pc =
sapex(flipfaces[0]);
9529 pd =
sapex(flipfaces[1]);
9541 printf(
" Performed %ld flips.\n", flipcount);
9579 int iloc,
int bowywat,
int rflag)
9581 face cavesh, neighsh, *parysh;
9582 face newsh, casout, casin;
9590 printf(
" Insert facet point %d.\n",
pointmark(insertpt));
9597 if ((splitseg != NULL) && (splitseg->
sh != NULL)) {
9599 spivot(*splitseg, *searchsh);
9605 if (searchsh->
sh == NULL) {
9609 loc =
slocate(insertpt, searchsh, 1, 1, rflag);
9619 *parysh = *searchsh;
9620 }
else if (loc ==
ONEDGE) {
9621 if ((splitseg != NULL) && (splitseg->
sh != NULL)) {
9622 splitseg->
shver = 0;
9623 pa =
sorg(*splitseg);
9625 pa =
sorg(*searchsh);
9627 if (searchsh->
sh != NULL) {
9629 neighsh = *searchsh;
9643 if (neighsh.
sh == searchsh->
sh)
break;
9644 if (neighsh.
sh == NULL)
break;
9657 neighsh = *searchsh;
9661 if (casout.
sh == NULL) {
9665 *searchsh = neighsh;
9676 pa =
sorg(*searchsh);
9677 pb =
sdest(*searchsh);
9691 sbond1(newsh, *searchsh);
9692 sbond1(*searchsh, newsh);
9694 if (casin.
sh != NULL) {
9706 neighsh = *searchsh;
9710 if (casout.
sh == NULL) {
9711 *searchsh = neighsh;
9718 pa =
sorg(*searchsh);
9719 pb =
sdest(*searchsh);
9722 if (ori >= 0)
break;
9724 }
else if (loc ==
INSTAR) {
9732 for (j = 0; j < 3; j++) {
9735 if (neighsh.
sh != NULL) {
9747 sapex(neighsh), insertpt);
9769 if ((
sorg(cavesh) == insertpt) || (
sdest(cavesh) == insertpt)) {
9798 pb =
sdest(*parysh);
9819 if (casout.
sh != NULL) {
9821 if (checkseg.
sh != NULL) {
9824 if (
sorg(newsh) !=
sorg(checkseg)) {
9829 while (neighsh.
sh != parysh->
sh) {
9837 if (checkseg.
sh != NULL) {
9845 if (newsh.
sh != NULL) {
9862 if (neighsh.
sh == NULL) {
9864 pb =
sdest(*parysh);
9869 if (neighsh.
sh == NULL)
break;
9873 if (neighsh.
sh != NULL) {
9877 sbond(newsh, neighsh);
9883 if (neighsh.
sh == NULL) {
9890 if (neighsh.
sh == NULL)
break;
9894 if (neighsh.
sh != NULL) {
9898 sbond(newsh, neighsh);
9903 if ((loc ==
ONEDGE) || ((splitseg != NULL) && (splitseg->
sh != NULL))
9910 face aseg, bseg, aoutseg, boutseg;
9917 if (
sapex(cavesh) == insertpt) {
9923 spivot(*parysh, neighsh);
9925 if (
sorg(neighsh) !=
sorg(cavesh)) {
9933 for (j = 0; j < 2; j++) {
9945 for (j = 0; j < 2; j++) {
9959 if ((splitseg != NULL) && (splitseg->
sh != NULL)) {
9965 pa =
sorg(*splitseg);
9966 pb =
sdest(*splitseg);
9988 if (boutseg.
sh != NULL) {
9990 sbond(boutseg, aoutseg);
9993 senext(*splitseg, aoutseg);
9995 if (aoutseg.
sh != NULL) {
9997 sbond(boutseg, aoutseg);
10002 sbond(aoutseg, boutseg);
10009 spivot(*parysh, neighsh);
10011 if (
sorg(neighsh) != pa) {
10040 if (
sapex(*parysh) == insertpt) {
10046 if ((splitseg != NULL) && (splitseg->
sh != NULL)) {
10082 face flipfaces[4], spinsh, *parysh;
10083 point pa, pb, pc, pd;
10087 if (parentseg != NULL) {
10090 face startsh, neighsh, nextsh;
10091 face abseg, prevseg, checkseg;
10092 face adjseg1, adjseg2;
10094 senext2(*parentseg, prevseg);
10098 pa =
sorg(prevseg);
10099 pb =
sdest(*parentseg);
10101 printf(
" Remove vertex %d from segment [%d, %d].\n",
10116 if (adjseg1.
sh != NULL) {
10120 sbond(adjseg1, adjseg2);
10123 senext(*parentseg, adjseg1);
10125 if (adjseg1.
sh != NULL) {
10129 sbond(adjseg1, adjseg2);
10137 spivot(*parentseg, *parentsh);
10138 if (parentsh->
sh != NULL) {
10139 spinsh = *parentsh;
10146 if (spinsh.
sh == NULL) {
10149 if (spinsh.
sh == parentsh->
sh)
break;
10158 if (
sorg(startsh) != delpt) {
10166 if (checkseg.
sh != NULL) {
10174 if (neighsh.
sh != startsh.
sh) {
10188 sbond(nextsh, startsh);
10190 sbond(nextsh, neighsh);
10198 spivot(neighsh, startsh);
10230 *parentseg = abseg;
10234 printf(
" Remove vertex %d from surface.\n",
pointmark(delpt));
10240 *parysh = *parentsh;
10251 if (
sorg(*parentsh) != delpt) {
10258 spinsh = *parentsh;
10264 if (spinsh.
sh == parentsh->
sh)
break;
10270 for (i = 0; i < 3; i++) {
10272 flipfaces[i] = *parysh;
10274 flip31(flipfaces, lawson);
10275 for (i = 0; i < 3; i++) {
10281 *parysh = flipfaces[3];
10289 flipfaces[0] = *parysh;
10290 spivot(flipfaces[0], flipfaces[1]);
10291 if (
sorg(flipfaces[0]) !=
sdest(flipfaces[1]))
10295 pa =
sorg(flipfaces[0]);
10296 pb =
sdest(flipfaces[0]);
10297 pc =
sapex(flipfaces[0]);
10298 pd =
sapex(flipfaces[1]);
10303 if (ori1 * ori2 < 0) {
10305 flip22(flipfaces, lawson, 0);
10308 senext2(flipfaces[1], *parentsh);
10311 *parysh = flipfaces[0];
10320 flipfaces[0] = *parysh;
10321 spivot(flipfaces[0], flipfaces[1]);
10322 if (
sorg(flipfaces[0]) !=
sdest(flipfaces[1])) {
10325 flip22(flipfaces, lawson, 0);
10326 senext2(flipfaces[1], *parentsh);
10329 *parysh = flipfaces[0];
10380 face* searchsh,
int aflag,
int cflag,
int rflag)
10385 enum {MOVE_BC, MOVE_CA} nextmove;
10386 REAL ori, ori_bc, ori_ca;
10389 pa =
sorg(*searchsh);
10390 pb =
sdest(*searchsh);
10391 pc =
sapex(*searchsh);
10402 }
else if (ori == 0.0) {
10408 for (i = 0; i < 3; i++) {
10409 pa =
sorg(*searchsh);
10410 pb =
sdest(*searchsh);
10412 if (ori > 0)
break;
10419 pc =
sapex(*searchsh);
10421 if (pc == searchpt) {
10435 nextmove = MOVE_CA;
10437 nextmove = MOVE_BC;
10441 nextmove = MOVE_BC;
10446 nextmove = MOVE_CA;
10472 if (nextmove == MOVE_BC) {
10483 spivot(*searchsh, neighsh);
10484 if (neighsh.sh == NULL) {
10488 if (
sorg(neighsh) !=
sdest(*searchsh)) {
10493 *searchsh = neighsh;
10494 pa =
sorg(*searchsh);
10495 pb =
sdest(*searchsh);
10496 pc =
sapex(*searchsh);
10498 if (pc == searchpt) {
10510 REAL n[3], area_abc, area_abp, area_bcp, area_cap;
10512 pa =
sorg(*searchsh);
10513 pb =
sdest(*searchsh);
10514 pc =
sapex(*searchsh);
10517 area_abc = sqrt(
dot(n, n));
10520 area_bcp = sqrt(
dot(n, n));
10521 if ((area_bcp / area_abc) <
b->
epsilon) {
10526 area_cap = sqrt(
dot(n, n));
10527 if ((area_cap / area_abc) <
b->
epsilon) {
10533 area_abp = sqrt(
dot(n, n));
10534 if ((area_abp / area_abc) <
b->
epsilon) {
10541 if (area_abp == 0) {
10542 if (area_bcp == 0) {
10546 if (area_cap == 0) {
10552 }
else if (area_bcp == 0) {
10553 if (area_cap == 0) {
10560 }
else if (area_cap == 0) {
10588 point endpt,
int insertsegflag,
int reporterrorflag,
int chkencflag)
10590 face flipshs[2], neighsh;
10591 point startpt, pa, pb, pc, pd;
10593 enum {MOVE_AB, MOVE_CA} nextmove;
10594 REAL ori_ab, ori_ca, len;
10597 startpt =
sorg(*searchsh);
10598 nextmove = MOVE_AB;
10601 printf(
" Scout segment (%d, %d).\n",
pointmark(startpt),
10609 pb =
sdest(*searchsh);
10615 pc =
sapex(*searchsh);
10640 nextmove = MOVE_CA;
10642 nextmove = MOVE_AB;
10645 nextmove = MOVE_AB;
10649 nextmove = MOVE_CA;
10677 if (nextmove == MOVE_AB) {
10684 spivot(*searchsh, neighsh);
10685 if (neighsh.sh != NULL) {
10687 senext(neighsh, *searchsh);
10695 *searchsh = neighsh;
10701 *searchsh = neighsh;
10708 *searchsh = neighsh;
10713 if (neighsh.sh != NULL) {
10715 *searchsh = neighsh;
10725 spivot(*searchsh, neighsh);
10727 senext(neighsh, *searchsh);
10733 if (insertsegflag) {
10740 ssbond(*searchsh, newseg);
10741 spivot(*searchsh, neighsh);
10742 if (neighsh.sh != NULL) {
10743 ssbond(neighsh, newseg);
10751 if (reporterrorflag) {
10753 printf(
"PLC Error: A vertex lies in a segment in facet #%d.\n",
10755 printf(
" Vertex: [%d] (%g,%g,%g).\n",
pointmark(pp),pp[0],pp[1],pp[2]);
10763 senext(*searchsh, flipshs[0]);
10765 if (reporterrorflag) {
10766 REAL P[3], Q[3], tp = 0, tq = 0;
10767 linelineint(startpt, endpt, pb, pc, P, Q, &tp, &tq);
10768 printf(
"PLC Error: Two segments intersect at point (%g,%g,%g),",
10770 printf(
" in facet #%d.\n",
shellmark(*searchsh));
10777 spivot(flipshs[0], flipshs[1]);
10781 pa =
sapex(flipshs[1]);
10782 pb =
sapex(flipshs[0]);
10783 pc =
sorg(flipshs[0]);
10784 pd =
sdest(flipshs[0]);
10791 }
else if (ori_ca <= 0) {
10795 *searchsh = flipshs[0];
10798 return sscoutsegment(searchsh, endpt, insertsegflag, reporterrorflag,
10812 face *parysh, searchsh, neighsh;
10822 searchsh = *parysh;
10823 searchsh.
shver = 0;
10824 for (j = 0; j < 3; j++) {
10825 spivot(searchsh, neighsh);
10827 if (neighsh.
sh != NULL) {
10840 *parysh = searchsh;
10849 for (i = 0; i < 3 * holes; i += 3) {
10851 loc =
slocate(&(holelist[i]), &searchsh, 1, 1, 0);
10855 *parysh = searchsh;
10862 searchsh = *parysh;
10863 searchsh.
shver = 0;
10864 for (j = 0; j < 3; j++) {
10865 spivot(searchsh, neighsh);
10866 if (neighsh.
sh != NULL) {
10907 badface *facelink = NULL, *newlinkitem, *f1, *f2;
10908 face *facperverlist, sface;
10909 face subsegloop, testseg;
10911 REAL ori1, ori2, ori3;
10913 REAL cosang, ang, ang_tol;
10918 printf(
" Unifying segments.\n");
10922 if (ang_tol < 0.0) ang_tol = 0.0;
10928 subsegloop.
shver = 0;
10932 torg =
sorg(subsegloop);
10933 tdest =
sdest(subsegloop);
10940 for (k = idx2faclist[idx]; k < idx2faclist[idx + 1]; k++) {
10941 sface = facperverlist[k];
10943 if (sface.
sh[3] == NULL)
continue;
10945 if (
sdest(sface) != tdest) {
10949 if (
sdest(sface) != tdest)
continue;
10966 }
else if (ori3 < 0) {
10973 }
else if (ori2 < 0) {
10988 }
else if (ori1 < 0) {
10992 }
else if (ori2 < 0) {
10998 }
else if (ori3 < 0) {
11021 }
else if (ori2 < 0) {
11029 if (
dot(n1, n2) > 0) {
11032 printf(
"%12.5E %12.5E %12.5E\n",torg[0],torg[1],torg[2]);
11033 printf(
"%12.5E %12.5E %12.5E\n",tdest[0],tdest[1],tdest[2]);
11034 printf(
"%12.5E %12.5E %12.5E\n",
sapex(f1->ss)[0],
sapex(f1->ss)[1],
sapex(f1->ss)[2]);
11037 printf(
"%12.5E %12.5E n1(%12.5E,%12.5E,%12.5E) n2(%12.5E,%12.5E,%12.5E)\n",
11038 ori1,ori2, n1[0], n1[1], n1[2], n2[0], n2[1], n2[2]);
11047 if (sface.
sh[3] != NULL) {
11050 newlinkitem->ss = sface;
11051 newlinkitem->nextitem = f1->nextitem;
11052 f1->nextitem = newlinkitem;
11062 if (
dot(n1, n2) > 0) {
11068 if (sface.
sh[3] != NULL) {
11071 newlinkitem->ss = sface;
11072 newlinkitem->nextitem = NULL;
11073 f1->nextitem = newlinkitem;
11078 newlinkitem->ss = sface;
11079 newlinkitem->nextitem = NULL;
11080 facelink = newlinkitem;
11091 if ((testseg.
sh != subsegloop.
sh) && (testseg.
sh[3] != NULL)) {
11095 ssbond(f1->ss, subsegloop);
11107 cosang =
dot(n1, n2) / (sqrt(
dot(n1, n1)) * sqrt(
dot(n2, n2)));
11109 if (cosang > 1.0) cosang = 1.0;
11110 else if (cosang < -1.0) cosang = -1.0;
11111 ang = acos(cosang);
11112 if (ang < ang_tol) {
11129 if (
b->
quality && (
in->segmentconstraintlist != (
REAL *) NULL)) {
11132 for (k = 0; k <
in->numberofsegmentconstraints; k++) {
11133 e1 = (int)
in->segmentconstraintlist[k * 3];
11134 e2 = (
int)
in->segmentconstraintlist[k * 3 + 1];
11137 len =
in->segmentconstraintlist[k * 3 + 2];
11147 delete [] idx2faclist;
11148 delete [] facperverlist;
11165 face* shperverlist;
11167 face searchsh, neighsh;
11168 face segloop, checkseg, newseg;
11169 point checkpt, pa = NULL, pb = NULL;
11178 printf(
"Inserting edges ...\n");
11185 for (i = 0; i <
in->numberofedges; i++) {
11186 endpts = &(
in->edgelist[(i << 1)]);
11187 if (endpts[0] == endpts[1]) {
11189 printf(
"Warning: Edge #%d is degenerated. Skipped.\n", i);
11195 edgemarker =
in->edgemarkerlist ?
in->edgemarkerlist[i] : -2;
11199 searchsh.
sh = NULL;
11200 idx = endpts[0] -
in->firstnumber;
11201 for (j = idx2shlist[idx]; j < idx2shlist[idx + 1]; j++) {
11202 checkpt =
sdest(shperverlist[j]);
11204 searchsh = shperverlist[j];
11207 checkpt =
sapex(shperverlist[j]);
11209 senext2(shperverlist[j], searchsh);
11216 if (searchsh.
sh != NULL) {
11219 if (checkseg.
sh != NULL) {
11224 pa =
sorg(searchsh);
11225 pb =
sdest(searchsh);
11228 ssbond(searchsh, newseg);
11229 spivot(searchsh, neighsh);
11230 if (neighsh.
sh != NULL) {
11231 ssbond(neighsh, newseg);
11237 pa = idx2verlist[endpts[0]];
11238 pb = idx2verlist[endpts[1]];
11241 printf(
"Warning: Edge #%d is degenerated. Skipped.\n", i);
11250 while (segloop.
sh != NULL) {
11251 ppt = (
point *) &(segloop.
sh[3]);
11252 if (((ppt[0] == pa) && (ppt[1] == pb)) ||
11253 ((ppt[0] == pb) && (ppt[1] == pa))) {
11260 if (newseg.
sh == NULL) {
11268 if (
b->
quality && (
in->segmentconstraintlist != (
REAL *) NULL)) {
11269 for (i = 0; i <
in->numberofsegmentconstraints; i++) {
11270 e1 = (int)
in->segmentconstraintlist[i * 3];
11271 e2 = (
int)
in->segmentconstraintlist[i * 3 + 1];
11274 len =
in->segmentconstraintlist[i * 3 + 2];
11282 delete [] shperverlist;
11283 delete [] idx2shlist;
11294 face parentsh, neighsh, neineish;
11296 point pa, pb, pc, pd;
11298 REAL cosang, cosang_tol;
11303 REAL *paryang = NULL;
11317 spivot(segloop, parentsh);
11318 if (parentsh.
sh != NULL) {
11319 spivot(parentsh, neighsh);
11320 if (neighsh.
sh != NULL) {
11321 spivot(neighsh, neineish);
11322 if (neineish.
sh == parentsh.
sh) {
11326 pa =
sorg(segloop);
11327 pb =
sdest(segloop);
11328 pc =
sapex(parentsh);
11329 pd =
sapex(neighsh);
11333 cosang =
dot(n1, n2) / (sqrt(
dot(n1, n1)) * sqrt(
dot(n2, n2)));
11334 if (cosang < cosang_tol) {
11343 dihedangarray->
newindex((
void **) &paryang);
11360 REAL cosang1, cosang2;
11366 for (i = 0; i < 3; i++) {
11368 senext(shloop, neighsh);
11373 pb =
sdest(shloop);
11374 pc =
sapex(shloop);
11375 for (j = 0; j < 3; j++) n1[j] = pa[j] - pb[j];
11376 for (j = 0; j < 3; j++) n2[j] = pc[j] - pb[j];
11377 cosang =
dot(n1, n2) / (sqrt(
dot(n1, n1)) * sqrt(
dot(n2, n2)));
11378 if (cosang > cosang_tol) {
11383 if (seg1.
sh[6] != NULL) {
11384 paryang = (
REAL *) (seg1.
sh[6]);
11385 cosang1 = *paryang;
11389 if (seg2.
sh[6] != NULL) {
11390 paryang = (
REAL *) (seg2.
sh[6]);
11391 cosang2 = *paryang;
11395 if (cosang1 < cosang_sep_tol) {
11396 if (cosang2 < cosang_sep_tol) {
11397 if (cosang1 < cosang2) {
11406 if (cosang2 < cosang_sep_tol) {
11410 if (segloop.
sh != NULL) {
11413 spivot(segloop, parentsh);
11414 spivot(parentsh, neighsh);
11430 delete dihedangarray;
11460 point pa, pb, pc, pd;
11461 enum {HMOVE, RMOVE, LMOVE} nextmove;
11462 REAL hori, rori, lori;
11467 pa =
org(*searchtet);
11470 decode(searchtet->tet[3], *searchtet);
11472 if ((
point) searchtet->tet[4] == pa) {
11473 searchtet->ver = 11;
11474 }
else if ((
point) searchtet->tet[5] == pa) {
11475 searchtet->ver = 3;
11476 }
else if ((
point) searchtet->tet[6] == pa) {
11477 searchtet->ver = 7;
11479 searchtet->ver = 0;
11483 pb =
dest(*searchtet);
11490 pc =
apex(*searchtet);
11498 size_t loopCount = 0;
11501 pd =
oppo(*searchtet);
11524 hori =
orient3d(pa, pb, pc, endpt);
11525 rori =
orient3d(pb, pa, pd, endpt);
11526 lori =
orient3d(pa, pc, pd, endpt);
11537 }
else if (s == 1) {
11619 if (nextmove == RMOVE) {
11621 }
else if (nextmove == LMOVE) {
11629 pb =
dest(*searchtet);
11630 pc =
apex(*searchtet);
11634 const size_t surelyHangingLoopCountThreshold = 1000000;
11635 if (loopCount > surelyHangingLoopCountThreshold)
11665 int level,
int edgepivot,
11670 int types[2], poss[4];
11675 if (
fc->seg[0] != NULL) {
11677 if (fliptype == 1) {
11682 for (i = 0; i < 3 && !rejflag; i++) {
11686 NULL, 1, types, poss);
11687 if (intflag == 2) {
11694 if (poss[0] == 0) {
11700 }
else if (intflag == 4) {
11704 if (poss[0] == 0) {
11713 }
else if (fliptype == 2) {
11719 if (intflag == 2) {
11726 }
else if (intflag == 4) {
11738 if ((
fc->fac[0] != NULL) && !rejflag) {
11740 if (fliptype == 1) {
11744 NULL, 1, types, poss);
11745 if (intflag == 2) {
11753 }
else if (intflag == 4) {
11756 for (i = 0; i < 2 && !rejflag; i++) {
11768 if ((
fc->remvert != NULL) && !rejflag) {
11771 if (fliptype == 1) {
11773 if ((pd ==
fc->remvert) || (pe ==
fc->remvert)) {
11779 if (
fc->remove_large_angle && !rejflag) {
11781 REAL cosmaxd = 0, diff;
11782 if (fliptype == 1) {
11790 diff = cosmaxd -
fc->cosdihed_in;
11791 if (fabs(diff/
fc->cosdihed_in) <
b->
epsilon) diff = 0.0;
11796 if (cosmaxd < fc->cosdihed_out) {
11797 fc->cosdihed_out = cosmaxd;
11801 diff = cosmaxd -
fc->cosdihed_in;
11802 if (fabs(diff/
fc->cosdihed_in) <
b->
epsilon) diff = 0.0;
11807 if (cosmaxd < fc->cosdihed_out) {
11808 fc->cosdihed_out = cosmaxd;
11813 }
else if (fliptype == 2) {
11821 diff = cosmaxd -
fc->cosdihed_in;
11822 if (fabs(diff/
fc->cosdihed_in) <
b->
epsilon) diff = 0.0;
11827 if (cosmaxd < fc->cosdihed_out) {
11828 fc->cosdihed_out = cosmaxd;
11832 diff = cosmaxd -
fc->cosdihed_in;
11833 if (fabs(diff/
fc->cosdihed_in) <
b->
epsilon) diff = 0.0;
11838 if (cosmaxd < fc->cosdihed_out) {
11839 fc->cosdihed_out = cosmaxd;
11845 if (edgepivot == 1) {
11850 diff = cosmaxd -
fc->cosdihed_in;
11851 if (fabs(diff/
fc->cosdihed_in) <
b->
epsilon) diff = 0.0;
11856 if (cosmaxd < fc->cosdihed_out) {
11857 fc->cosdihed_out = cosmaxd;
11866 diff = cosmaxd -
fc->cosdihed_in;
11867 if (fabs(diff/
fc->cosdihed_in) <
b->
epsilon) diff = 0.0;
11872 if (cosmaxd < fc->cosdihed_out) {
11873 fc->cosdihed_out = cosmaxd;
11908 if (
fc->collectencsegflag) {
11909 face checkseg, *paryseg;
11915 *paryseg = checkseg;
11924 spintet = *flipedge;
11928 if (spintet.
tet == flipedge->
tet)
break;
11943 spintet = *flipedge;
11946 abtets[i] = spintet;
11950 if (spintet.
tet == flipedge->
tet)
break;
11960 for (i = 0; i < nn; i++) {
11964 *flipedge = abtets[0];
11969 int bakunflip =
fc->unflip;
11972 fc->unflip = bakunflip;
11993 triface fliptets[3], flipedge;
11994 point pa, pb, pc, pd, pe;
11998 fliptets[0] = *flipface;
11999 fsym(*flipface, fliptets[1]);
12000 pa =
org(fliptets[0]);
12001 pb =
dest(fliptets[0]);
12002 pc =
apex(fliptets[0]);
12003 pd =
oppo(fliptets[0]);
12004 pe =
oppo(fliptets[1]);
12015 eprev(*flipface, flipedge);
12018 enext(*flipface, flipedge);
12021 flipedge = *flipface;
12060 triface* searchtet,
int fullsearch)
12065 fc.seg[0] = startpt;
12067 fc.checkflipeligibility = 1;
12076 if (
dest(*searchtet) == endpt) {
12126 point pa, pb, pc, pd;
12129 int types[2], poss[4], pos = 0;
12148 neightet = *searchtet;
12149 j = (neightet.
ver & 3);
12150 for (i = j + 1; i < j + 4; i++) {
12151 neightet.
ver = (i % 4);
12152 pa =
org(neightet);
12153 pb =
dest(neightet);
12154 pc =
apex(neightet);
12155 pd =
oppo(neightet);
12156 if (
tri_edge_test(pa,pb,pc,startpt,endpt, pd, 1, types, poss)) {
12172 for (i = 0; i < 2; i++) {
12178 pa =
org(neightet);
12179 pb =
dest(neightet);
12180 pc =
apex(neightet);
12181 pd =
oppo(neightet);
12182 if (
tri_edge_test(pa,pb,pc,startpt,endpt,pd,1, types, poss)) {
12203 for (i = 0; i < pos; i++) {
12209 pd =
org(neightet);
12219 *searchtet = neightet;
12222 bakface.
forg =
org(*searchtet);
12268 if ((searchtet->
tet == NULL) ||
12269 (
org(*searchtet) != bakface.
forg) ||
12277 if (
dest(*searchtet) == bakface.
fdest) {
12278 spintet = *searchtet;
12282 *searchtet = spintet;
12286 if (spintet.
tet == searchtet->
tet) {
12287 searchtet->
tet = NULL;
12291 if (searchtet->
tet != NULL) {
12292 if (
oppo(*searchtet) != bakface.
foppo) {
12294 if (
oppo(*searchtet) != bakface.
foppo) {
12296 searchtet->
tet = NULL;
12302 searchtet->
tet = NULL;
12305 searchtet->
tet = NULL;
12307 if (searchtet->
tet == NULL) {
12348 point pc, pd, steinerpt;
12351 REAL vcd[3], sampt[3], smtpt[3];
12352 REAL maxminvol = 0.0, minvol = 0.0, ori;
12353 int success, maxidx = 0;
12357 pc =
apex(abtets[0]);
12358 pd =
oppo(abtets[n-1]);
12365 for (i = 0; i < n; i++) {
12368 *parytet = worktet;
12371 *parytet = worktet;
12378 for (i = 0; i < 3; i++) vcd[i] = pd[i] - pc[i];
12381 for (it = 1; it < N; it++) {
12382 for (i = 0; i < 3; i++) {
12383 sampt[i] = pc[i] + (stepi * (double) it) * vcd[i];
12391 if (minvol > ori) minvol = ori;
12395 maxminvol = minvol;
12398 if (maxminvol < minvol) {
12399 maxminvol = minvol;
12405 if (maxminvol <= 0) {
12410 for (i = 0; i < 3; i++) {
12411 smtpt[i] = pc[i] + (stepi * (double) maxidx) * vcd[i];
12419 *parytet = faketet1;
12423 *parytet = faketet2;
12460 for (i = 0; i < 3; i++) steinerpt[i] = smtpt[i];
12463 for (i = 0; i < n; i++) {
12466 *parytet = abtets[i];
12468 worktet = abtets[0];
12474 locate(steinerpt, &(abtets[0]));
12475 worktet = abtets[0];
12480 if (
insertpoint(steinerpt, &worktet, NULL, NULL, &ivf)) {
12501 face *paryseg, candseg;
12502 point startpt, endpt, pc, pd;
12505 REAL P[3], Q[3], tp, tq;
12506 REAL len, smlen = 0, split = 0, split_q = 0;
12510 startpt =
sorg(*misseg);
12511 endpt =
sdest(*misseg);
12513 fc.seg[0] = startpt;
12515 fc.checkflipeligibility = 1;
12516 fc.collectencsegflag = 1;
12539 pc =
sorg(*paryseg);
12540 pd =
sdest(*paryseg);
12542 if (
linelineint(startpt, endpt, pc, pd, P, Q, &tp, &tq)) {
12545 if ((tp > 0) && (tq < 1)) {
12547 if (tp < (
b->
epsilon * 1e+3)) tp = 0.0;
12549 if ((1.0 - tp) < (
b->
epsilon * 1e+3)) tp = 1.0;
12552 if ((tp <= 0) || (tp >= 1))
continue;
12553 if ((tq > 0) && (tq < 1)) {
12555 if (tq < (
b->
epsilon * 1e+3)) tq = 0.0;
12557 if ((1.0 - tq) < (
b->
epsilon * 1e+3)) tq = 1.0;
12560 if ((tq <= 0) || (tq >= 1))
continue;
12567 candseg = *paryseg;
12573 candseg = *paryseg;
12589 point steinerpt, *parypt;
12595 for (i = 0; i < 3; i++) {
12596 steinerpt[i] = startpt[i] + split * (endpt[i] - startpt[i]);
12599 for (i = 0; i < 3; i++) {
12600 P[i] = startpt[i] + split * (endpt[i] - startpt[i]);
12602 pc =
sorg(candseg);
12603 pd =
sdest(candseg);
12604 for (i = 0; i < 3; i++) {
12605 Q[i] = pc[i] + split_q * (pd[i] - pc[i]);
12608 for (i = 0; i < 3; i++) {
12609 steinerpt[i] = 0.5 * (P[i] + Q[i]);
12620 splitseg = *misseg;
12621 spivot(*misseg, splitsh);
12624 splitseg.
sh = NULL;
12638 if (!
insertpoint(steinerpt, &searchtet, &splitsh, &splitseg, &ivf)) {
12647 *parypt = steinerpt;
12652 *paryseg = *misseg;
12671 triface *abtets, searchtet, spintet;
12674 point startpt, endpt;
12675 point pa, pb, pd, steinerpt, *parypt;
12678 int types[2], poss[4];
12679 int n, endi, success;
12683 startpt =
sorg(*misseg);
12686 startpt =
sorg(*misseg);
12688 endpt =
sdest(*misseg);
12699 fsym(searchtet, spintet);
12700 pd =
oppo(spintet);
12701 for (i = 0; i < 3; i++) {
12703 pb =
dest(spintet);
12704 if (
tri_edge_test(pa, pb, pd, startpt, endpt, NULL, 1, types, poss)) {
12713 spintet = searchtet;
12717 if (
apex(spintet) == endpt) {
12722 if (spintet.
tet == searchtet.
tet)
break;
12729 spintet = searchtet;
12730 for (i = 0; i < n; i++) {
12731 abtets[i] = spintet;
12759 if ((n - endi) > 2) {
12781 *paryseg = *misseg;
12786 if (!splitsegflag) {
12791 printf(
" Splitting segment (%d, %d)\n",
pointmark(startpt),
12810 if (steinerpt == NULL) {
12813 for (i = 0; i < 3; i++) {
12814 steinerpt[i] = 0.5 * (startpt[i] + endpt[i]);
12818 spivot(*misseg, splitsh);
12830 if (!
insertpoint(steinerpt, &searchtet, &splitsh, misseg, &ivf)) {
12838 *parypt = steinerpt;
12862 face sseg, *paryseg;
12863 point startpt, endpt;
12871 printf(
" Recover segments [%s level = %2d] #: %ld.\n",
12886 if (searchtet.
tet != NULL) {
12890 startpt =
sorg(sseg);
12891 endpt =
sdest(sseg);
12894 printf(
" Recover segment (%d, %d).\n",
pointmark(startpt),
12909 if (!success && fullsearch) {
12920 spintet = searchtet;
12924 }
while (spintet.
tet != searchtet.
tet);
12926 if (steinerflag > 0) {
12931 if (!success && (steinerflag > 1)) {
12938 if (misseglist != NULL) {
12940 misseglist->
newindex((
void **) &paryseg);
12952 printf(
" Add %ld Steiner points in volume.\n",
12956 printf(
" Add %ld Steiner points in segments.\n",
12988 int types[2], poss[4], intflag;
12997 fc.checkflipeligibility = 1;
13000 for (i = 0; i < 3 && !success; i++) {
13006 spintet = *searchtet;
13008 if (
apex(spintet) ==
fc.fac[(i+2)%3]) {
13010 *searchtet = spintet;
13012 for (j = i; j > 0; j--) {
13019 if (spintet.
tet == searchtet->
tet)
break;
13021 if (success)
break;
13023 flipedge.
tet = NULL;
13025 spintet = *searchtet;
13027 pd =
apex(spintet);
13028 pe =
oppo(spintet);
13031 intflag =
tri_edge_test(pa, pb, pc, pd, pe, NULL, 1, types, poss);
13035 if (intflag == 2) {
13038 if (searchsh != NULL) {
13045 intflag, types, poss);
13052 if (chkface.
tet == flipedge.
tet)
break;
13057 intflag, types, poss);
13062 point touchpt, *parypt;
13063 if (poss[1] == 0) {
13071 face checksh, *parysh;
13072 int siloc = (int)
ONFACE;
13086 spivot(*parysh, checksh);
13088 if (checksh.
sh[3] != NULL) {
13103 searchsh->
sh = NULL;
13108 intflag, types, poss);
13114 intflag, types, poss);
13124 if (spintet.
tet == searchtet->
tet) {
13148 triface searchtet, neightet, spintet;
13149 face searchsh, neighsh, neineish, *parysh;
13151 point startpt, endpt, apexpt, *parypt;
13159 printf(
" Recover subfaces [%s level = %2d] #: %ld.\n",
13170 searchsh = *parysh;
13172 if (searchsh.
sh[3] == NULL)
continue;
13175 if (neightet.
tet != NULL)
continue;
13179 printf(
" Recover subface (%d, %d, %d).\n",
pointmark(
sorg(searchsh)),
13184 for (i = 0; i < 3; i++) {
13185 sspivot(searchsh, bdsegs[i]);
13186 if (bdsegs[i].sh != NULL) {
13189 if (searchtet.
tet == NULL) {
13196 startpt =
sorg(searchsh);
13197 endpt =
sdest(searchsh);
13200 if (
dest(searchtet) == endpt) {
13218 ssbond(searchsh, bdsegs[i]);
13219 spivot(searchsh, neighsh);
13220 if (neighsh.
sh != NULL) {
13221 ssbond(neighsh, bdsegs[i]);
13226 spintet = searchtet;
13230 }
while (spintet.
tet != searchtet.
tet);
13234 for (j = (i - 1); j >= 0; j--) {
13236 spivot(bdsegs[j], neineish);
13238 spivot(neineish, neighsh);
13239 if (neighsh.
sh != NULL) {
13243 spintet = searchtet;
13247 if (spintet.
tet == searchtet.
tet)
break;
13255 printf(
" Add a Steiner point in subedge (%d, %d).\n",
13259 for (j = 0; j < 3; j++) {
13260 steinerpt[j] = 0.5 * (startpt[j] + endpt[j]);
13275 if (!
insertpoint(steinerpt, &searchtet, &searchsh, NULL, &ivf)) {
13281 *parypt = steinerpt;
13294 startpt =
sorg(searchsh);
13295 endpt =
sdest(searchsh);
13296 apexpt =
sapex(searchsh);
13301 for (j = 0; j < 3; j++) {
13303 spivot(bdsegs[j], neineish);
13305 spivot(neineish, neighsh);
13306 if (neighsh.
sh != NULL) {
13310 spintet = neightet;
13314 if (spintet.
tet == neightet.
tet)
break;
13321 if (searchsh.
sh != NULL) {
13323 tsbond(searchtet, searchsh);
13326 tsbond(searchtet, searchsh);
13332 printf(
" Add a Steiner point in subface (%d, %d, %d).\n",
13336 for (j = 0; j < 3; j++) {
13337 steinerpt[j] = (startpt[j] + endpt[j] + apexpt[j]) / 3.0;
13352 if (!
insertpoint(steinerpt, &searchtet, &searchsh, NULL, &ivf)) {
13358 *parypt = steinerpt;
13369 if (misshlist != NULL) {
13371 misshlist->
newindex((
void **) &parysh);
13372 *parysh = searchsh;
13401 triface searchtet, neightet, *parytet;
13402 face checksh, *parysh;
13414 tetlist->
newindex((
void **) &parytet);
13415 *parytet = searchtet;
13416 if (vertlist != NULL) {
13418 j = (searchtet.
ver & 3);
13419 for (i = 1; i < 4; i++) {
13420 pt = (
point) searchtet.
tet[4 + ((j + i) % 4)];
13422 vertlist->
newindex((
void **) &parypt);
13428 esym(searchtet, neightet);
13430 if (shlist != NULL) {
13435 shlist->
newindex((
void **) &parysh);
13448 tetlist->
newindex((
void **) &parytet);
13449 *parytet = neightet;
13450 if (vertlist != NULL) {
13452 pt =
apex(neightet);
13454 vertlist->
newindex((
void **) &parypt);
13460 for (i = 0; i < tetlist->
objects; i++) {
13465 for (j = 0; j < 2; j++) {
13468 esym(searchtet, neightet);
13470 if (shlist != NULL) {
13475 shlist->
newindex((
void **) &parysh);
13488 tetlist->
newindex((
void **) &parytet);
13489 *parytet = neightet;
13490 if (vertlist != NULL) {
13492 pt =
apex(neightet);
13495 vertlist->
newindex((
void **) &parypt);
13506 for (i = 0; i < tetlist->
objects; i++) {
13511 if (vertlist != NULL) {
13512 for (i = 0; i < vertlist->
objects; i++) {
13518 if (shlist != NULL) {
13519 for (i = 0; i < shlist->
objects; i++) {
13525 return (
int) tetlist->
objects;
13542 triface searchtet, neightet, *parytet;
13553 if (
org(*tedge) == e1) {
13554 if (
dest(*tedge) == e2) {
13557 }
else if (
org(*tedge) == e2) {
13558 if (
dest(*tedge) == e1) {
13568 if (
dest(*tedge) == e2) {
13574 if (
dest(*tedge) == e1) {
13587 for (i = 0; i < 3; i++) {
13588 pt =
apex(searchtet);
13598 fnext(searchtet, neightet);
13601 pt =
apex(neightet);
13610 tetlist->
newindex((
void **) &parytet);
13611 *parytet = searchtet;
13613 tetlist->
newindex((
void **) &parytet);
13614 *parytet = neightet;
13618 for (i = 0; (i < tetlist->
objects) && !done; i++) {
13620 searchtet = *parytet;
13621 for (j = 0; (j < 2) && !done; j++) {
13623 fnext(searchtet, neightet);
13626 pt =
apex(neightet);
13633 tetlist->
newindex((
void **) &parytet);
13634 *parytet = neightet;
13641 for (i = 0; i < tetlist->
objects; i++) {
13661 point *pendpt, *parypt;
13669 fc.remvert = startpt;
13670 fc.checkflipeligibility = 1;
13676 for (i = 0; i < endptlist->
objects; i++) {
13684 if (
getedge(startpt, *pendpt, &searchtet)) {
13695 if (
dest(searchtet) == *pendpt) {
13726 return (
int) endptlist->
objects;
13747 triface *fliptets = NULL, wrktets[4];
13748 triface searchtet, spintet, neightet;
13749 face parentsh, spinsh, checksh;
13750 face leftseg, rightseg, checkseg;
13751 point lpt = NULL, rpt = NULL, apexpt;
13755 int valence, removeflag;
13765 if (
sdest(leftseg) == steinerpt) {
13766 senext(leftseg, rightseg);
13768 rightseg.
shver = 0;
13770 rightseg = leftseg;
13775 lpt =
sorg(leftseg);
13776 rpt =
sdest(rightseg);
13778 printf(
" Removing Steiner point %d in segment (%d, %d).\n",
13784 printf(
" Removing Steiner point %d in facet.\n",
13789 printf(
" Removing Steiner point %d in volume.\n",
13794 printf(
" Removing a point %d in volume.\n",
13814 if (valence == 4) {
13820 }
else if (valence == 5) {
13824 if (
org(searchtet) != steinerpt) {
13828 neightet.
tet = NULL;
13829 spintet = searchtet;
13832 if (
apex(spintet) == rpt) {
13834 neightet = spintet;
13837 if (spintet.
tet == searchtet.
tet)
break;
13841 }
else if (i == 4) {
13844 if (
apex(neightet) == rpt) {
13848 esym(neightet, searchtet);
13851 wrktets[0] = searchtet;
13852 for (i = 0; i < 2; i++) {
13853 fnext(wrktets[i], wrktets[i+1]);
13855 if (
apex(wrktets[0]) ==
oppo(wrktets[2])) {
13866 wrktets[0] = searchtet;
13867 wrktets[1] = searchtet;
13870 wrktets[2] = searchtet;
13875 searchtet.
tet = NULL;
13876 for (i = 0; i < 3; i++) {
13877 spintet = wrktets[i];
13882 if (spintet.
tet == wrktets[i].tet)
break;
13884 if (valence == 3) {
13886 searchtet = wrktets[i];
13904 if (
org(searchtet) != steinerpt) {
13907 spintet = searchtet;
13910 eprev(spintet, neightet);
13914 if (
oppo(neightet) != rpt) {
13919 if (spintet.
tet == searchtet.
tet) {
13932 if (
getedge(lpt, rpt, &searchtet)) {
13940 for (i = 0; i < 2; i++) {
13941 checkseg = (i == 0) ? leftseg : rightseg;
13943 spintet = neightet;
13947 if (spintet.
tet == neightet.
tet)
break;
13952 spivot(rightseg, parentsh);
13958 spintet = searchtet;
13962 if (spintet.
tet == searchtet.
tet)
break;
13980 for (i = 0; i < 2; i++) {
13981 checkseg = (i == 0) ? leftseg : rightseg;
13983 spintet = neightet;
13987 if (spintet.
tet == neightet.
tet)
break;
13993 for (i = 0; i < 2; i++) {
13994 checkseg = (i == 0) ? leftseg : rightseg;
13995 spivot(checkseg, parentsh);
13996 if (parentsh.
sh != NULL) {
14000 if (neightet.
tet != NULL) {
14005 if (neightet.
tet != NULL) {
14010 if (spinsh.
sh == parentsh.
sh)
break;
14020 fliptets[0] = searchtet;
14021 for (i = 0; i < 2; i++) {
14022 fnext(fliptets[i], fliptets[i+1]);
14024 eprev(fliptets[0], fliptets[3]);
14032 }
else if (loc ==
ONFACE) {
14037 fliptets[0] = searchtet;
14038 for (i = 0; i < 2; i++) {
14039 fnext(fliptets[i], fliptets[i+1]);
14041 eprev(fliptets[0], fliptets[3]);
14045 for (i = 3; i < 5; i++) {
14046 fnext(fliptets[i], fliptets[i+1]);
14053 for (i = 3; i < 6; i++) {
14059 for (i = 0; i < 3; i++) {
14060 esym(fliptets[i+3], wrktets[i]);
14061 esym(fliptets[i], fliptets[i+3]);
14062 fliptets[i] = wrktets[i];
14065 wrktets[1] = fliptets[1];
14066 fliptets[1] = fliptets[2];
14067 fliptets[2] = wrktets[1];
14068 wrktets[1] = fliptets[4];
14069 fliptets[4] = fliptets[5];
14070 fliptets[5] = wrktets[1];
14087 }
else if (loc ==
ONEDGE) {
14093 spintet = searchtet;
14097 if (spintet.
tet == searchtet.
tet)
break;
14100 fliptets =
new triface[2 * n];
14101 fliptets[0] = searchtet;
14102 for (i = 0; i < (n - 1); i++) {
14103 fnext(fliptets[i], fliptets[i+1]);
14105 eprev(fliptets[0], fliptets[n]);
14109 for (i = n; i < (2 * n - 1); i++) {
14110 fnext(fliptets[i], fliptets[i+1]);
14120 wrktets[0] = fliptets[0];
14124 wrktets[1] = fliptets[n];
14131 fliptets[n] = wrktets[2];
14133 fliptets[0] = wrktets[0];
14142 for (i = 1; i < (n - 1); i++) {
14143 wrktets[0] = wrktets[1];
14147 wrktets[1] = fliptets[n+i];
14149 wrktets[2] = fliptets[i];
14155 fliptets[i] = wrktets[0];
14166 wrktets[3] = wrktets[1];
14167 wrktets[0] = fliptets[n];
14171 wrktets[1] = fliptets[n-1];
14174 wrktets[2] = fliptets[2*n-1];
14181 fliptets[n-1] = wrktets[0];
14185 delete [] fliptets;
14193 spivot(rightseg, parentsh);
14197 rightseg.
shver = 0;
14202 spintet = searchtet;
14206 if (spintet.
tet == searchtet.
tet)
break;
14211 spivot(rightseg, parentsh);
14212 if (parentsh.
sh != NULL) {
14215 if (
sorg(spinsh) != lpt) {
14218 apexpt =
sapex(spinsh);
14220 spintet = searchtet;
14222 if (
apex(spintet) == apexpt) {
14223 tsbond(spintet, spinsh);
14225 fsym(spintet, neightet);
14226 tsbond(neightet, spinsh);
14233 if (spinsh.
sh == parentsh.
sh)
break;
14270 face parentsh, spinsh, *parysh;
14271 face leftseg, rightseg;
14272 point lpt = NULL, rpt = NULL;
14280 if (
sdest(leftseg) == steinerpt) {
14281 senext(leftseg, rightseg);
14283 rightseg.
shver = 0;
14285 rightseg = leftseg;
14290 lpt =
sorg(leftseg);
14291 rpt =
sdest(rightseg);
14293 printf(
" Suppressing Steiner point %d in segment (%d, %d).\n",
14297 spivot(leftseg, parentsh);
14298 if (parentsh.
sh != NULL) {
14307 if (spinsh.
sh == NULL)
break;
14308 if (spinsh.
sh == parentsh.
sh)
break;
14318 printf(
" Suppressing Steiner point %d from facet.\n",
14323 for (i = 0; i < 2; i++) {
14325 *parysh = parentsh;
14332 triface searchtet, neightet, *parytet;
14333 point pa, pb, pc, pd;
14334 REAL v1[3], v2[3], len, u;
14336 REAL startpt[3] = {0,}, samplept[3] = {0,}, candpt[3] = {0,};
14337 REAL ori, minvol, smallvol;
14343 for (i = 0; i < n; i++) newsteiners[i] = NULL;
14357 pa =
sorg(*parysh);
14358 pb =
sdest(*parysh);
14359 pc =
sapex(*parysh);
14361 len = sqrt(
dot(v1, v1));
14367 pd =
sapex(*parysh);
14369 len = sqrt(
dot(v2, v2));
14374 v1[0] = 0.5 * (v1[0] + v2[0]);
14375 v1[1] = 0.5 * (v1[1] + v2[1]);
14376 v1[2] = 0.5 * (v1[2] + v2[2]);
14382 v2[0] = steinerpt[0] + len * v1[0];
14383 v2[1] = steinerpt[1] + len * v1[1];
14384 v2[2] = steinerpt[2] + len * v1[2];
14387 pa =
org(*parytet);
14388 pb =
dest(*parytet);
14389 pc =
apex(*parytet);
14393 ori =
orient3d(steinerpt, pa, pb, v2);
14395 ori =
orient3d(steinerpt, pb, pc, v2);
14397 ori =
orient3d(steinerpt, pc, pa, v2);
14411 *parytet = neightet;
14416 v1[0] = steinerpt[0] - startpt[0];
14417 v1[1] = steinerpt[1] - startpt[1];
14418 v1[2] = steinerpt[2] - startpt[2];
14421 for (j = 1; j < samplesize - 1; j++) {
14422 samplept[0] = startpt[0] + ((
REAL) j / (
REAL) samplesize) * v1[0];
14423 samplept[1] = startpt[1] + ((
REAL) j / (
REAL) samplesize) * v1[1];
14424 samplept[2] = startpt[2] + ((
REAL) j / (
REAL) samplesize) * v1[2];
14429 pa =
org(*parytet);
14430 pb =
dest(*parytet);
14431 pc =
apex(*parytet);
14432 ori =
orient3d(pb, pa, pc, samplept);
14436 if (smallvol == -1) {
14439 if (ori < smallvol) smallvol = ori;
14444 if (minvol == -1.0) {
14445 candpt[0] = samplept[0];
14446 candpt[1] = samplept[1];
14447 candpt[2] = samplept[2];
14450 if (minvol < smallvol) {
14452 candpt[0] = samplept[0];
14453 candpt[1] = samplept[1];
14454 candpt[2] = samplept[2];
14465 if (minvol > 0)
break;
14469 if (minvol == -1.0) {
14477 newsteiners[i][0] = candpt[0];
14478 newsteiners[i][1] = candpt[1];
14479 newsteiners[i][2] = candpt[2];
14484 if (i < cavesegshlist->objects) {
14486 for (; i > 0; i--) {
14487 if (newsteiners[i - 1] != NULL) {
14491 delete [] newsteiners;
14497 triface newtet, newface, spintet;
14498 face newsh, neighsh;
14499 face *splitseg, checkseg;
14507 spintet = neightet;
14511 if (spintet.
tet == neightet.
tet)
break;
14514 spintet = neightet;
14518 if (spintet.
tet == neightet.
tet)
break;
14536 setoppo(*parytet, newsteiners[i]);
14562 spivot(rightseg, parentsh);
14563 splitseg = &rightseg;
14565 if (
sdest(parentsh) == steinerpt) {
14567 }
else if (
sapex(parentsh) == steinerpt) {
14576 rightseg.
shver = 0;
14592 bond(newtet, neightet);
14595 tsbond(neightet, newsh);
14602 spivot(rightseg, parentsh);
14610 spivot(spinsh, neighsh);
14620 bond(newtet, neightet);
14623 if (spinsh.
sh == parentsh.
sh)
break;
14633 for (k = 0; k < 2; k++) {
14635 for (j = 0; j < 3; j++) {
14637 esym(newtet, newface);
14638 if (newface.
tet[newface.
ver & 3] == NULL) {
14641 if (checkseg.
sh != NULL) {
14647 if (neighsh.
sh != NULL) {
14655 spintet = neightet;
14657 esym(spintet, searchtet);
14658 fsym(searchtet, spintet);
14659 if (spintet.
tet == NULL)
break;
14662 neightet = searchtet;
14670 spintet = neightet;
14672 esym(spintet, searchtet);
14673 fsym(searchtet, spintet);
14674 if (spintet.
tet == NULL)
break;
14677 neightet = searchtet;
14679 pc =
apex(newface);
14680 if (
apex(neightet) == steinerpt) {
14684 *parytet = neightet;
14697 bond(newface, neightet);
14735 int steinercount = 0;
14741 for (i = 0; i < n; i++) {
14742 if (newsteiners[i] != NULL) {
14747 *parypt = newsteiners[i];
14756 if (steinercount > 0) {
14758 printf(
" Added %d interior Steiner points.\n", steinercount);
14762 delete [] newsteiners;
14782 printf(
"Suppressing Steiner points ...\n");
14785 point rempt, *parypt;
14789 int suppcount = 0, remcount = 0;
14806 if (suppcount > 0) {
14808 printf(
" Suppressed %d boundary Steiner points.\n", suppcount);
14826 if (remcount > 0) {
14828 printf(
" Removed %d interior Steiner points.\n", remcount);
14840 int smtcount, count, ivcount;
14867 ppt = (
point *) &(parytet->
tet[4]);
14913 printf(
"BUG Report! The mesh contain inverted elements.\n");
14917 if (smtcount > 0) {
14918 printf(
" Smoothed %d Steiner points.\n", smtcount);
14938 face searchsh, *parysh;
14939 face searchseg, *paryseg;
14940 point rempt, *parypt;
14946 long bak_segref_count, bak_facref_count, bak_volref_count;
14949 printf(
"Recovering boundaries...\n");
14954 printf(
" Recovering segments.\n");
14973 *paryseg = searchseg;
14987 if (misseglist->
objects > 0) {
14991 if (misseglist->
objects >= ms) {
15004 for (i = 0; i < misseglist->
objects; i++) {
15018 printf(
" %ld (%ld) segments are recovered (missing).\n",
15022 if (misseglist->
objects > 0) {
15024 while (misseglist->
objects > 0) {
15026 for (i = 0; i < misseglist->
objects; i++) {
15034 if (misseglist->
objects < ms) {
15042 printf(
" %ld (%ld) segments are recovered (missing).\n",
15047 if (misseglist->
objects > 0) {
15050 while (misseglist->
objects > 0) {
15052 for (i = 0; i < misseglist->
objects; i++) {
15060 if (misseglist->
objects < ms) {
15072 if (misseglist->
objects > 0) {
15076 for (i = 0; i < misseglist->
objects; i++) {
15087 printf(
" Added another %ld Steiner points in volume.\n",
15104 bdrysteinerptlist->
newindex((
void **) &parypt);
15111 printf(
" Suppressed %ld Steiner points in segments.\n",
15115 bak_segref_count) {
15116 printf(
" Removed %ld Steiner points in segments.\n",
15129 printf(
" Recovering facets.\n");
15147 *parysh = searchsh;
15160 if (misshlist->
objects > 0) {
15164 if (misshlist->
objects >= ms) {
15177 for (i = 0; i < misshlist->
objects; i++) {
15191 printf(
" %ld (%ld) subfaces are recovered (missing).\n",
15195 if (misshlist->
objects > 0) {
15197 for (i = 0; i < misshlist->
objects; i++) {
15220 bdrysteinerptlist->
newindex((
void **) &parypt);
15226 printf(
" Removed %ld Steiner points in facets.\n",
15234 if (bdrysteinerptlist->
objects > 0) {
15236 printf(
" %ld Steiner points remained in boundary.\n",
15246 delete bdrysteinerptlist;
15270 triface tetloop, neightet, *parytet, *parytet1;
15272 face checksh, *parysh;
15274 point ptloop, *parypt;
15280 printf(
"Marking exterior tetrahedra ...\n");
15282 printf(
"Removing exterior tetrahedra ...\n");
15300 hullarray->
newindex((
void **) &parytet);
15301 *parytet = tetloop;
15306 tetarray->
newindex((
void **) &parytet);
15307 *parytet = neightet;
15314 if (
in->numberofholes > 0) {
15316 for (i = 0; i < 3 *
in->numberofholes; i += 3) {
15318 neightet.
tet = NULL;
15324 tetarray->
newindex((
void **) &parytet);
15325 *parytet = neightet;
15332 hullarray->
newindex((
void **) &parytet);
15334 tetarray->
newindex((
void **) &parytet);
15336 *parytet = tetloop;
15347 hullarray->
newindex((
void **) &parytet);
15348 *parytet = tetloop;
15363 printf(
"Warning: The %d-th hole point ", i/3 + 1);
15364 printf(
"lies outside the convex hull.\n");
15373 regiontets =
new triface[
in->numberofregions];
15375 for (i = 0; i < 5 *
in->numberofregions; i += 5) {
15377 neightet.
tet = NULL;
15380 regiontets[i/5] = neightet;
15383 printf(
"Warning: The %d-th region point ", i/5+1);
15384 printf(
"lies outside the convex hull.\n");
15386 regiontets[i/5].
tet = NULL;
15392 for (i = 0; i < tetarray->
objects; i++) {
15394 j = (parytet->
ver & 3);
15396 for (k = 1; k < 4; k++) {
15397 decode(parytet->
tet[(j + k) % 4], neightet);
15404 tetarray->
newindex((
void **) &parytet1);
15405 *parytet1 = neightet;
15411 hullarray->
newindex((
void **) &parytet1);
15412 *parytet1 = neightet;
15438 for (i = 0; i <
in->numberofregions; i++) {
15441 printf(
"Warning: The %d-th region point ", i+1);
15442 printf(
"lies in the exterior of the domain.\n");
15444 regiontets[i].
tet = NULL;
15455 while (ptloop != NULL) {
15466 (
in->numberofpoints - (
in->firstnumber ? 0 : 1))) {
15479 face segloop, *paryseg;
15481 long delsegcount = 0l;
15487 while (segloop.
sh != NULL) {
15491 *paryseg = segloop;
15499 for (i = 0; i < tetarray->
objects; i++) {
15501 for (j = 0; j < 4; j++) {
15508 pb =
dest(tetloop);
15509 pc =
apex(tetloop);
15511 bond(tetloop, hulltet);
15514 tsbond(hulltet, checksh);
15516 for (k = 0; k < 3; k++) {
15531 newhullfacearray->
newindex((
void **) &parytet1);
15532 parytet1->tet = parytet->
tet;
15539 for (i = 0; i < newhullfacearray->
objects; i++) {
15541 fsym(*parytet, neightet);
15543 fsym(neightet, hulltet);
15544 for (j = 0; j < 3; j++) {
15545 esym(hulltet, casface);
15546 if (casface.
tet[casface.
ver & 3] == NULL) {
15550 neightet = *parytet;
15561 bond(casface, neightet);
15572 face casingout, casingin;
15578 printf(
"Warning: Removed an exterior face (%d, %d, %d) #%d\n",
15584 for (j = 0; j < 3; j++) {
15585 spivot(*parysh, casingout);
15587 if (casingout.
sh != NULL) {
15588 casingin = casingout;
15590 spivot(casingin, checksh);
15591 if (checksh.
sh == parysh->
sh)
break;
15592 casingin = checksh;
15594 if (casingin.
sh != casingout.
sh) {
15596 sbond1(casingin, casingout);
15601 if (checkseg.
sh != NULL) {
15603 ssbond(casingout, checkseg);
15606 if (checkseg.
sh != NULL) {
15608 if (delsegcount == 0) {
15610 printf(
"Warning: Removed an exterior segment (%d, %d) #%d\n",
15633 if (paryseg->
sh && (paryseg->
sh[3] != NULL)) {
15637 printf(
"Warning: Removed an exterior segment (%d, %d) #%d\n",
15649 if (delsegcount > 0) {
15651 printf(
" Deleted %ld segments.\n", delsegcount);
15658 long delsteinercount = 0l;
15666 (
in->numberofpoints - (
in->firstnumber ? 0 : 1))) {
15685 if (delsteinercount > 0l) {
15686 if (
unuverts > (delvertcount + delsteinercount)) {
15687 printf(
" Removed %ld exterior input vertices.\n",
15688 unuverts - delvertcount - delsteinercount);
15690 printf(
" Removed %ld exterior Steiner vertices.\n",
15693 printf(
" Removed %ld exterior input vertices.\n",
15707 for (i = 0; i < tetarray->
objects; i++) {
15713 for (i = 0; i < hullarray->
objects; i++) {
15719 delete newhullfacearray;
15726 for (i = 0; i < tetarray->
objects; i++) {
15732 for (i = 0; i < hullarray->
objects; i++) {
15753 printf(
"Spreading region attributes.\n");
15756 int attr, maxattr = 0;
15760 int regioncount = 0;
15763 if (
in->numberofregions > 0) {
15765 for (i = 0; i < 5 *
in->numberofregions; i += 5) {
15766 if (regiontets[i/5].tet != NULL) {
15767 attr = (int)
in->regionlist[i + 3];
15768 if (attr > maxattr) {
15771 volume =
in->regionlist[i + 4];
15773 infect(regiontets[i/5]);
15774 tetarray->
newindex((
void **) &parytet);
15775 *parytet = regiontets[i/5];
15777 for (j = 0; j < tetarray->
objects; j++) {
15779 tetloop = *parytet;
15784 for (k = 0; k < 4; k++) {
15791 tetarray->
newindex((
void **) &parytet);
15792 *parytet = neightet;
15803 attr = maxattr + 1;
15811 tetarray->
newindex((
void **) &parytet);
15812 *parytet = tetloop;
15814 for (j = 0; j < tetarray->
objects; j++) {
15816 tetloop = *parytet;
15818 for (k = 0; k < 4; k++) {
15825 tetarray->
newindex((
void **) &parytet);
15826 *parytet = neightet;
15848 if (regioncount > 1) {
15849 printf(
" Found %d subdomains.\n", regioncount);
15851 printf(
" Found %d domain.\n", regioncount);
15856 if (regiontets != NULL) {
15857 delete [] regiontets;
15872 fsym(tetloop, neightet);
15882 if (sliver_peel_count > 0l) {
15884 printf(
" Removed %ld hull slivers.\n", sliver_peel_count);
15902 *queface = *chkface;
15933 triface fliptets[5], neightet, hulltet;
15934 face checksh, casingout;
15936 point pd, pe, *pts;
15939 long flipcount, totalcount = 0l;
15940 long sliver_peels = 0l;
15944 size_t loopCount = 0;
15955 fliptets[0] = popface->
tt;
15968 fsym(fliptets[0], fliptets[1]);
15972 tspivot(fliptets[0], checksh);
15973 for (i = 0; i < 3; i++) {
15975 spivot(checksh, casingout);
15978 stpivot(casingout, neightet);
15979 if (neightet.
tet == fliptets[0].
tet) {
15985 pe =
org(neightet);
15994 pd =
dest(neightet);
16001 fliptets[0] = neightet;
16002 fnext(fliptets[0], fliptets[1]);
16003 fnext(fliptets[1], fliptets[2]);
16009 if (
fc->remove_ndelaunay_edge) {
16013 fc->tetprism_vol_sum = 0.0;
16032 pts = (
point *) fliptets[1].tet;
16033 sign =
insphere_s(pts[4], pts[5], pts[6], pts[7],
oppo(fliptets[0]));
16037 pd =
oppo(fliptets[0]);
16038 pe =
oppo(fliptets[1]);
16043 len3 = (len3 * len3 * len3);
16048 for (i = 0; i < 3; i++) {
16053 esym(fliptets[0], fliptets[2]);
16054 esym(fliptets[1], fliptets[3]);
16057 if ((fabs(vol) / len3) <
b->
epsilon) {
16062 if (ori <= 0)
break;
16073 if (
fc->remove_ndelaunay_edge) {
16077 fc->tetprism_vol_sum = 0.0;
16085 if (
issubseg(fliptets[0]))
continue;
16089 for (i = 0; i < 3; i++) {
16090 fnext(fliptets[i], fliptets[i+1]);
16092 if (fliptets[3].tet == fliptets[0].tet) {
16096 if (
fc->remove_ndelaunay_edge) {
16100 fc->tetprism_vol_sum = 0.0;
16105 fnext(fliptets[3], fliptets[4]);
16106 if (fliptets[4].tet == fliptets[0].tet) {
16127 fnext(fliptets[3], fliptets[1]);
16128 fnext(fliptets[1], fliptets[2]);
16141 if (
fc->remove_ndelaunay_edge) {
16145 fc->tetprism_vol_sum = 0.0;
16155 bface->tt = fliptets[0];
16156 bface->forg =
org(fliptets[0]);
16157 bface->fdest =
dest(fliptets[0]);
16158 bface->fapex =
apex(fliptets[0]);
16163 if (flipcount > 0) {
16164 printf(
" Performed %ld flips.\n", flipcount);
16168 totalcount += flipcount;
16173 if (flipcount == 0l)
break;
16189 const size_t surelyHangingLoopCountThreshold = 10000;
16190 if (loopCount > surelyHangingLoopCountThreshold)
16199 if (totalcount > 0) {
16200 printf(
" Performed %ld flips.\n", totalcount);
16202 if (sliver_peels > 0) {
16203 printf(
" Removed %ld hull slivers.\n", sliver_peels);
16210 return totalcount + sliver_peels;
16221 arraypool *flipqueue, *nextflipqueue, *swapqueue;
16222 triface tetloop, neightet, *parytet;
16229 printf(
"Recovering Delaunayness...\n");
16237 while (tetloop.
tet != NULL) {
16238 for (tetloop.
ver = 0; tetloop.
ver < 4; tetloop.
ver++) {
16244 ppt = (
point *) &(tetloop.
tet[4]);
16258 printf(
" Recover Delaunay [Lawson] : %ld\n",
flippool->
items);
16262 fc.remove_ndelaunay_edge = 1;
16276 fc.collectnewtets = 1;
16290 swapqueue = flipqueue;
16294 while (flipqueue->
objects > 0l) {
16297 printf(
" Recover Delaunay [level = %2d] #: %ld.\n",
16301 for (i = 0; i < flipqueue->
objects; i++) {
16306 fc.tetprism_vol_sum = 0.0;
16312 for (parytet->
ver = 0; parytet->
ver < 4; parytet->
ver++) {
16330 flipqueue->
newindex((
void **) &parybface);
16331 *parybface = *bface;
16336 nextflipqueue->
newindex((
void **) &parybface);
16337 *parybface = *bface;
16341 fc.tetprism_vol_sum = 0.0;
16353 swapqueue = flipqueue;
16354 flipqueue = nextflipqueue;
16355 nextflipqueue = swapqueue;
16357 if (flipqueue->
objects > 0l) {
16367 if (flipqueue->
objects > 0l) {
16369 printf(
" %ld non-Delaunay edges remained.\n", flipqueue->
objects);
16379 delete nextflipqueue;
16394 if (
getedge(pa, pb, searchtet)) {
16395 spintet = *searchtet;
16397 if (
apex(spintet) == pc) {
16398 *searchtet = spintet;
16402 if (spintet.
tet == searchtet->
tet)
break;
16404 if (
apex(*searchtet) == pc) {
16405 if (
oppo(*searchtet) == pd) {
16409 if (
oppo(*searchtet) == pd) {
16427 arraypool *flipqueue, *nextflipqueue, *swapqueue;
16432 REAL *cosdd, ncosdd[6], maxdd;
16433 long totalremcount, remcount;
16451 fc.remove_large_angle = 1;
16453 fc.collectnewtets = 1;
16454 fc.checkflipeligibility = 1;
16456 totalremcount = 0l;
16459 swapqueue = flipqueue;
16463 while (flipqueue->
objects > 0l) {
16467 while (flipqueue->
objects > 0l) {
16469 printf(
" Improving mesh qualiy by flips [%d]#: %ld.\n",
16473 for (k = 0; k < flipqueue->
objects; k++) {
16479 if (bface->
tt.
ver != 11) {
16484 &bface->
key, NULL);
16485 bface->
forg = ppt[0];
16486 bface->
fdest = ppt[1];
16487 bface->
fapex = ppt[2];
16488 bface->
foppo = ppt[3];
16489 bface->
tt.
ver = 11;
16491 if (bface->
key == 0) {
16495 &bface->
key, NULL);
16497 cosdd = bface->
cent;
16499 for (i = 0; (i < 6) && !remflag; i++) {
16503 fc.cosdihed_in = cosdd[i];
16504 fc.cosdihed_out = 0.0;
16514 ppt = (
point *) & (parytet->
tet[4]);
16521 nextflipqueue->
newindex((
void **) &parybface);
16522 parybface->tt.tet = parytet->
tet;
16523 parybface->tt.ver = 11;
16524 parybface->forg = ppt[0];
16525 parybface->fdest = ppt[1];
16526 parybface->fapex = ppt[2];
16527 parybface->foppo = ppt[3];
16528 parybface->key = maxdd;
16529 for (n = 0; n < 6; n++) {
16530 parybface->cent[n] = ncosdd[n];
16545 *parybface = *bface;
16553 swapqueue = flipqueue;
16554 flipqueue = nextflipqueue;
16555 nextflipqueue = swapqueue;
16559 printf(
" Removed %ld bad tets.\n", remcount);
16561 totalremcount += remcount;
16573 swapqueue = flipqueue;
16584 delete nextflipqueue;
16586 return totalremcount;
16618 triface *parytet, *parytet1, swaptet;
16620 REAL fcent[3], startpt[3], nextpt[3], bestpt[3];
16621 REAL oldval, minval = 0.0, val;
16628 numdirs = (int) linkfacelist->
objects;
16637 for (i = 0; i < 3; i++) {
16638 bestpt[i] = startpt[i] = smtpt[i];
16647 for (i = 0; i < numdirs; i++) {
16652 pa =
org(*parytet);
16653 pb =
dest(*parytet);
16654 pc =
apex(*parytet);
16655 for (j = 0; j < 3; j++) {
16656 fcent[j] = (pa[j] + pb[j] + pc[j]) / 3.0;
16658 for (j = 0; j < 3; j++) {
16659 nextpt[j] = startpt[j] + opm->
searchstep * (fcent[j] - startpt[j]);
16662 for (j = 0; j < linkfacelist->
objects; j++) {
16665 pa =
org(*parytet);
16666 pb =
dest(*parytet);
16668 pb =
org(*parytet);
16669 pa =
dest(*parytet);
16671 pc =
apex(*parytet);
16672 ori =
orient3d(pa, pb, pc, nextpt);
16682 if (maxcosd < -1) maxcosd = -1.0;
16683 val = maxcosd + 1.0;
16700 if (val <= opm->imprval) {
16707 minval = (val < minval) ? val : minval;
16711 if (j == linkfacelist->
objects) {
16715 for (j = 0; j < 3; j++) bestpt[j] = nextpt[j];
16718 j = linkfacelist->
objects - i - 1;
16721 swaptet = *parytet1;
16722 *parytet1 = *parytet;
16723 *parytet = swaptet;
16726 diff = opm->
imprval - oldval;
16732 if ((diff / oldval) < 1e-3) diff = 0.0;
16745 for (j = 0; j < 3; j++) startpt[j] = bestpt[j];
16761 for (i = 0; i < 3; i++) smtpt[i] = startpt[i];
16780 long totalsmtcount, smtcount;
16788 swapqueue = flipqueue;
16792 totalsmtcount = 0l;
16795 while (flipqueue->
objects > 0l) {
16800 printf(
" Improving mesh quality by smoothing [%d]#: %ld.\n",
16804 for (k = 0; k < flipqueue->
objects; k++) {
16814 &bface->
key, NULL);
16819 for (i = 0; (i < 4) && !smtflag; i++) {
16841 ppt = (
point *) & (parytet->
tet[4]);
16843 bface->
cent, &bface->
key, NULL);
16848 parybface->tt = *parytet;
16849 parybface->forg = ppt[0];
16850 parybface->fdest = ppt[1];
16851 parybface->fapex = ppt[2];
16852 parybface->foppo = ppt[3];
16853 parybface->tt.ver = 11;
16854 parybface->key = 0.0;
16867 parybface->tt = bface->
tt;
16868 parybface->forg = ppt[0];
16869 parybface->fdest = ppt[1];
16870 parybface->fapex = ppt[2];
16871 parybface->foppo = ppt[3];
16872 parybface->tt.ver = 11;
16873 parybface->key = 0.0;
16889 printf(
" Smooth %ld points.\n", smtcount);
16891 totalsmtcount += smtcount;
16893 if (smtcount == 0l) {
16904 swapqueue = flipqueue;
16911 return totalsmtcount;
16923 triface searchtet, spintet, *parytet;
16924 point pa, pb, steinerpt;
16927 REAL smtpt[3], midpt[3];
16943 spintet = searchtet;
16949 if (spintet.
tet == searchtet.
tet)
break;
16957 spintet = searchtet;
16958 for (i = 0; i < n; i++) {
16959 abtets[i] = spintet;
16964 for (i = 0; i < n; i++) {
16965 eprev(abtets[i], searchtet);
16968 *parytet = searchtet;
16969 enext(abtets[i], searchtet);
16972 *parytet = searchtet;
16976 pa =
org(abtets[0]);
16977 pb =
dest(abtets[0]);
16978 for (i = 0; i < 3; i++) {
16979 smtpt[i] = midpt[i] = 0.5 * (pa[i] + pb[i]);
17013 for (i = 0; i < 3; i++) steinerpt[i] = smtpt[i];
17016 for (i = 0; i < n; i++) {
17019 *parytet = abtets[i];
17022 searchtet = abtets[0];
17024 locate(steinerpt, &searchtet);
17034 if (
insertpoint(steinerpt, &searchtet, NULL, NULL, &ivf)) {
17058 REAL cosdd[6], maxcosd;
17059 long totalsptcount, sptcount;
17066 swapqueue = flipqueue;
17070 totalsptcount = 0l;
17078 printf(
" Splitting bad quality tets [%d]#: %ld.\n",
17086 if ((bface->
key == 0) || (bface->
tt.
ver != 11)) {
17091 &bface->
key, NULL);
17097 for (j = 0; j < 6; j++) {
17111 while (parytet != NULL) {
17113 ppt = (
point *) & (parytet->
tet[4]);
17119 parybface->forg = ppt[0];
17120 parybface->fdest = ppt[1];
17121 parybface->fapex = ppt[2];
17122 parybface->foppo = ppt[3];
17123 parybface->tt.tet = parytet->
tet;
17124 parybface->tt.ver = 11;
17125 parybface->key = maxcosd;
17126 for (i = 0; i < 6; i++) {
17127 parybface->cent[i] = cosdd[i];
17136 *parybface = *bface;
17145 printf(
" Split %ld tets.\n", sptcount);
17147 totalsptcount += sptcount;
17149 if (sptcount == 0l) {
17160 swapqueue = flipqueue;
17167 return totalsptcount;
17183 REAL ncosdd[6], maxdd;
17184 long totalremcount, remcount;
17185 long totalsmtcount, smtcount;
17186 long totalsptcount, sptcount;
17192 printf(
"Optimizing mesh...\n");
17198 printf(
" Optimization level = %d.\n",
b->
optlevel);
17199 printf(
" Optimization scheme = %d.\n",
b->
optscheme);
17200 printf(
" Number of iteration = %d.\n", optpasses);
17204 totalsmtcount = totalsptcount = totalremcount = 0l;
17215 while (checktet.
tet != NULL) {
17223 ppt = (
point *) & (checktet.
tet[4]);
17224 tetalldihedral(ppt[0], ppt[1], ppt[2], ppt[3], ncosdd, &maxdd, NULL);
17228 parybface->tt.tet = checktet.
tet;
17229 parybface->tt.ver = 11;
17230 parybface->forg = ppt[0];
17231 parybface->fdest = ppt[1];
17232 parybface->fapex = ppt[2];
17233 parybface->foppo = ppt[3];
17234 parybface->key = maxdd;
17235 for (n = 0; n < 6; n++) {
17236 parybface->cent[n] = ncosdd[n];
17248 sizeof(
void *), 0);
17259 while (iter < optpasses) {
17260 smtcount = sptcount = remcount = 0l;
17263 totalsmtcount += smtcount;
17264 if (smtcount > 0l) {
17266 totalremcount += remcount;
17272 totalsptcount += sptcount;
17273 if (sptcount > 0l) {
17275 totalremcount += remcount;
17280 if (remcount > 0l) {
17302 if (totalremcount > 0l) {
17303 printf(
" Removed %ld edges.\n", totalremcount);
17305 if (totalsmtcount > 0l) {
17306 printf(
" Smoothed %ld points.\n", totalsmtcount);
17308 if (totalsptcount > 0l) {
17309 printf(
" Split %ld slivers.\n", totalsptcount);
17335 int oldidx, newidx;
17339 printf(
"Jettisoning redundant points.\n");
17344 oldidx = newidx = 0;
17346 while (pointloop != (
point) NULL) {
17356 if (
in->pointmarkerlist != (
int *) NULL) {
17357 if (oldidx < in->numberofpoints) {
17359 in->pointmarkerlist[newidx] =
in->pointmarkerlist[oldidx];
17368 printf(
" %ld duplicated vertices are removed.\n",
dupverts);
17369 printf(
" %ld unused vertices are removed.\n",
unuverts);