1 /* ----------------------------------------------------------------------
2 SPARTA - Stochastic PArallel Rarefied-gas Time-accurate Analyzer
3 http://sparta.sandia.gov
4 Steve Plimpton, sjplimp@sandia.gov, Michael Gallis, magalli@sandia.gov
5 Sandia National Laboratories
6
7 Copyright (2014) Sandia Corporation. Under the terms of Contract
8 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
9 certain rights in this software. This software is distributed under
10 the GNU General Public License.
11
12 See the README file in the top-level SPARTA directory.
13 ------------------------------------------------------------------------- */
14
15 #include "math.h"
16 #include "string.h"
17 #include "marching_cubes.h"
18 #include "grid.h"
19 #include "surf.h"
20 #include "irregular.h"
21 #include "lookup_table.h"
22 #include "geometry.h"
23 #include "my_page.h"
24 #include "memory.h"
25 #include "error.h"
26
27 // DEBUG
28 #include "update.h"
29
30 using namespace SPARTA_NS;
31
32 // prototype for non-class function
33
34 int compare_indices(const void *, const void *);
35
36 enum{UNKNOWN,OUTSIDE,INSIDE,OVERLAP}; // several files
37 enum{NCHILD,NPARENT,NUNKNOWN,NPBCHILD,NPBPARENT,NPBUNKNOWN,NBOUND}; // Grid
38
39 #define DELTA 128
40 #define BIG 1.0e20
41 #define EPSILON 1.0e-16
42
43 /* ---------------------------------------------------------------------- */
44
MarchingCubes(SPARTA * sparta,int ggroup_caller,double thresh_caller)45 MarchingCubes::MarchingCubes(SPARTA *sparta, int ggroup_caller,
46 double thresh_caller) :
47 Pointers(sparta)
48 {
49 MPI_Comm_rank(world,&me);
50
51 ggroup = ggroup_caller;
52 thresh = thresh_caller;
53 }
54
55 /* ----------------------------------------------------------------------
56 create 2d implicit surfs from grid point values
57 follows https://en.wikipedia.org/wiki/Marching_squares
58 see 2 sections: Basic algorithm and Disambiguation of saddle points
59 treating open circles as flow volume, solid circles as material
60 NOTE: Wiki page numbers points counter-clockwise
61 SPARTA numbers them in x, then in y
62 so bit2 and bit3 are swapped below
63 this gives case #s here consistent with Wiki page
64 process each grid cell independently
65 4 corner points open/solid -> 2^4 = 16 cases
66 cases infer 0,1,2 line segments in each grid cell
67 order 2 points in each line segment to give normal into flow volume
68 treat two saddle point cases (my 9,6) (Wiki 5,10)
69 based on ave value at cell center
70 ------------------------------------------------------------------------- */
71
invoke(double ** cvalues,int * svalues,int ** mcflags)72 void MarchingCubes::invoke(double **cvalues, int *svalues, int **mcflags)
73 {
74 int i,j,ipt,isurf,nsurf,icase,which;
75 surfint *ptr;
76
77 Grid::ChildCell *cells = grid->cells;
78 Grid::ChildInfo *cinfo = grid->cinfo;
79 MyPage<surfint> *csurfs = grid->csurfs;
80 int nglocal = grid->nlocal;
81 int groupbit = grid->bitmask[ggroup];
82
83 for (int icell = 0; icell < nglocal; icell++) {
84 if (!(cinfo[icell].mask & groupbit)) continue;
85 if (cells[icell].nsplit <= 0) continue;
86 lo = cells[icell].lo;
87 hi = cells[icell].hi;
88
89 // nsurf = # of tris in cell
90 // cvalues[8] = 8 corner point values, each is 0 to 255 inclusive
91 // thresh = value between 0 and 255 to threshhold on
92 // lo[3] = lower left corner pt of grid cell
93 // hi[3] = upper right corner pt of grid cell
94 // pt = list of 3*nsurf points that are the corner pts of each tri
95
96 // cvalues are ordered
97 // bottom-lower-left, bottom-lower-right,
98 // bottom-upper-left, bottom-upper-right
99 // top-lower-left, top-lower-right, top-upper-left, top-upper-right
100 // Vzyx encodes this as 0/1 in each dim
101
102 v000 = cvalues[icell][0];
103 v001 = cvalues[icell][1];
104 v010 = cvalues[icell][2];
105 v011 = cvalues[icell][3];
106 v100 = cvalues[icell][4];
107 v101 = cvalues[icell][5];
108 v110 = cvalues[icell][6];
109 v111 = cvalues[icell][7];
110
111 v000iso = v000 - thresh;
112 v001iso = v001 - thresh;
113 v010iso = v010 - thresh;
114 v011iso = v011 - thresh;
115 v100iso = v100 - thresh;
116 v101iso = v101 - thresh;
117 v110iso = v110 - thresh;
118 v111iso = v111 - thresh;
119
120 // make bits 2, 3, 6 and 7 consistent with Lewiner paper (see NOTE above)
121
122 bit0 = v000 <= thresh ? 0 : 1;
123 bit1 = v001 <= thresh ? 0 : 1;
124 bit2 = v011 <= thresh ? 0 : 1;
125 bit3 = v010 <= thresh ? 0 : 1;
126 bit4 = v100 <= thresh ? 0 : 1;
127 bit5 = v101 <= thresh ? 0 : 1;
128 bit6 = v111 <= thresh ? 0 : 1;
129 bit7 = v110 <= thresh ? 0 : 1;
130
131 which = (bit7 << 7) + (bit6 << 6) + (bit5 << 5) + (bit4 << 4) +
132 (bit3 << 3) + (bit2 << 2) + (bit1 << 1) + bit0;
133
134 // icase = case of the active cube in [0..15]
135
136 icase = cases[which][0];
137 config = cases[which][1];
138 subconfig = 0;
139
140 switch (icase) {
141 case 0:
142 nsurf = 0;
143 break;
144
145 case 1:
146 nsurf = add_triangle(tiling1[config], 1);
147 break;
148
149 case 2:
150 nsurf = add_triangle(tiling2[config], 2);
151 break;
152
153 case 3:
154 if (test_face(test3[config]))
155 nsurf = add_triangle(tiling3_2[config], 4); // 3.2
156 else
157 nsurf = add_triangle(tiling3_1[config], 2); // 3.1
158 break;
159
160 case 4:
161 if (modified_test_interior(test4[config],icase))
162 nsurf = add_triangle(tiling4_1[config], 2); // 4.1.1
163 else
164 nsurf = add_triangle(tiling4_2[config], 6); // 4.1.2
165 break;
166
167 case 5:
168 nsurf = add_triangle(tiling5[config], 3);
169 break;
170
171 case 6:
172 if (test_face(test6[config][0]))
173 nsurf = add_triangle(tiling6_2[config], 5); // 6.2
174 else {
175 if (modified_test_interior(test6[config][1],icase))
176 nsurf = add_triangle(tiling6_1_1[config], 3); // 6.1.1
177 else {
178 nsurf = add_triangle(tiling6_1_2[config], 9); // 6.1.2
179 }
180 }
181 break;
182
183 case 7:
184 if (test_face(test7[config][0])) subconfig += 1;
185 if (test_face(test7[config][1])) subconfig += 2;
186 if (test_face(test7[config][2])) subconfig += 4;
187 switch (subconfig) {
188 case 0:
189 nsurf = add_triangle(tiling7_1[config], 3); break;
190 case 1:
191 nsurf = add_triangle(tiling7_2[config][0], 5); break;
192 case 2:
193 nsurf = add_triangle(tiling7_2[config][1], 5); break;
194 case 3:
195 nsurf = add_triangle(tiling7_3[config][0], 9); break;
196 case 4:
197 nsurf = add_triangle(tiling7_2[config][2], 5); break;
198 case 5:
199 nsurf = add_triangle(tiling7_3[config][1], 9); break;
200 case 6:
201 nsurf = add_triangle(tiling7_3[config][2], 9); break;
202 case 7:
203 if (test_interior(test7[config][3],icase))
204 nsurf = add_triangle(tiling7_4_2[config], 9);
205 else
206 nsurf = add_triangle(tiling7_4_1[config], 5);
207 break;
208 };
209 break;
210
211 case 8:
212 nsurf = add_triangle(tiling8[config], 2);
213 break;
214
215 case 9:
216 nsurf = add_triangle(tiling9[config], 4);
217 break;
218
219 case 10:
220 if (test_face(test10[config][0])) {
221 if (test_face(test10[config][1]))
222 nsurf = add_triangle(tiling10_1_1_[config], 4); // 10.1.1
223 else {
224 nsurf = add_triangle(tiling10_2[config], 8); // 10.2
225 }
226 } else {
227 if (test_face(test10[config][1])) {
228 nsurf = add_triangle(tiling10_2_[config], 8); // 10.2
229 } else {
230 if (test_interior(test10[config][2],icase))
231 nsurf = add_triangle(tiling10_1_1[config], 4); // 10.1.1
232 else
233 nsurf = add_triangle(tiling10_1_2[config], 8); // 10.1.2
234 }
235 }
236 break;
237
238 case 11:
239 nsurf = add_triangle(tiling11[config], 4);
240 break;
241
242 case 12:
243 if (test_face(test12[config][0])) {
244 if (test_face(test12[config][1]))
245 nsurf = add_triangle(tiling12_1_1_[config], 4); // 12.1.1
246 else {
247 nsurf = add_triangle(tiling12_2[config], 8); // 12.2
248 }
249 } else {
250 if (test_face(test12[config][1])) {
251 nsurf = add_triangle(tiling12_2_[config], 8); // 12.2
252 } else {
253 if (test_interior(test12[config][2],icase))
254 nsurf = add_triangle(tiling12_1_1[config], 4); // 12.1.1
255 else
256 nsurf = add_triangle(tiling12_1_2[config], 8); // 12.1.2
257 }
258 }
259 break;
260
261 case 13:
262 if (test_face(test13[config][0])) subconfig += 1;
263 if (test_face(test13[config][1])) subconfig += 2;
264 if (test_face(test13[config][2])) subconfig += 4;
265 if (test_face(test13[config][3])) subconfig += 8;
266 if (test_face(test13[config][4])) subconfig += 16;
267 if (test_face(test13[config][5])) subconfig += 32;
268
269 switch (subconfig13[subconfig]) {
270 case 0:/* 13.1 */
271 nsurf = add_triangle(tiling13_1[config], 4); break;
272
273 case 1:/* 13.2 */
274 nsurf = add_triangle(tiling13_2[config][0], 6); break;
275 case 2:/* 13.2 */
276 nsurf = add_triangle(tiling13_2[config][1], 6); break;
277 case 3:/* 13.2 */
278 nsurf = add_triangle(tiling13_2[config][2], 6); break;
279 case 4:/* 13.2 */
280 nsurf = add_triangle(tiling13_2[config][3], 6); break;
281 case 5:/* 13.2 */
282 nsurf = add_triangle(tiling13_2[config][4], 6); break;
283 case 6:/* 13.2 */
284 nsurf = add_triangle(tiling13_2[config][5], 6); break;
285
286 case 7:/* 13.3 */
287 nsurf = add_triangle(tiling13_3[config][0], 10); break;
288 case 8:/* 13.3 */
289 nsurf = add_triangle(tiling13_3[config][1], 10); break;
290 case 9:/* 13.3 */
291 nsurf = add_triangle(tiling13_3[config][2], 10); break;
292 case 10:/* 13.3 */
293 nsurf = add_triangle(tiling13_3[config][3], 10); break;
294 case 11:/* 13.3 */
295 nsurf = add_triangle(tiling13_3[config][4], 10); break;
296 case 12:/* 13.3 */
297 nsurf = add_triangle(tiling13_3[config][5], 10); break;
298 case 13:/* 13.3 */
299 nsurf = add_triangle(tiling13_3[config][6], 10); break;
300 case 14:/* 13.3 */
301 nsurf = add_triangle(tiling13_3[config][7], 10); break;
302 case 15:/* 13.3 */
303 nsurf = add_triangle(tiling13_3[config][8], 10); break;
304 case 16:/* 13.3 */
305 nsurf = add_triangle(tiling13_3[config][9], 10); break;
306 case 17:/* 13.3 */
307 nsurf = add_triangle(tiling13_3[config][10], 10); break;
308 case 18:/* 13.3 */
309 nsurf = add_triangle(tiling13_3[config][11], 10); break;
310
311 case 19:/* 13.4 */
312 nsurf = add_triangle(tiling13_4[config][0], 12); break;
313 case 20:/* 13.4 */
314 nsurf = add_triangle(tiling13_4[config][1], 12); break;
315 case 21:/* 13.4 */
316 nsurf = add_triangle(tiling13_4[config][2], 12); break;
317 case 22:/* 13.4 */
318 nsurf = add_triangle(tiling13_4[config][3], 12); break;
319
320 case 23:/* 13.5 */
321 subconfig = 0;
322 if (interior_test_case13())
323 nsurf = add_triangle(tiling13_5_1[config][0], 6);
324 else
325 nsurf = add_triangle(tiling13_5_2[config][0], 10);
326 break;
327
328 case 24:/* 13.5 */
329 subconfig = 1;
330 if (interior_test_case13())
331 nsurf = add_triangle(tiling13_5_1[config][1], 6);
332 else
333 nsurf = add_triangle(tiling13_5_2[config][1], 10);
334 break;
335
336 case 25:/* 13.5 */
337 subconfig = 2;
338 if (interior_test_case13())
339 nsurf = add_triangle(tiling13_5_1[config][2], 6);
340 else
341 nsurf = add_triangle(tiling13_5_2[config][2], 10);
342 break;
343
344 case 26:/* 13.5 */
345 subconfig = 3;
346 if (interior_test_case13())
347 nsurf = add_triangle(tiling13_5_1[config][3], 6);
348 else
349 nsurf = add_triangle(tiling13_5_2[config][3], 10);
350 break;
351
352 case 27:/* 13.3 */
353 nsurf = add_triangle(tiling13_3_[config][0], 10); break;
354 case 28:/* 13.3 */
355 nsurf = add_triangle(tiling13_3_[config][1], 10); break;
356 case 29:/* 13.3 */
357 nsurf = add_triangle(tiling13_3_[config][2], 10); break;
358 case 30:/* 13.3 */
359 nsurf = add_triangle(tiling13_3_[config][3], 10); break;
360 case 31:/* 13.3 */
361 nsurf = add_triangle(tiling13_3_[config][4], 10); break;
362 case 32:/* 13.3 */
363 nsurf = add_triangle(tiling13_3_[config][5], 10); break;
364 case 33:/* 13.3 */
365 nsurf = add_triangle(tiling13_3_[config][6], 10); break;
366 case 34:/* 13.3 */
367 nsurf = add_triangle(tiling13_3_[config][7], 10); break;
368 case 35:/* 13.3 */
369 nsurf = add_triangle(tiling13_3_[config][8], 10); break;
370 case 36:/* 13.3 */
371 nsurf = add_triangle(tiling13_3_[config][9], 10); break;
372 case 37:/* 13.3 */
373 nsurf = add_triangle(tiling13_3_[config][10], 10); break;
374 case 38:/* 13.3 */
375 nsurf = add_triangle(tiling13_3_[config][11], 10); break;
376
377 case 39:/* 13.2 */
378 nsurf = add_triangle(tiling13_2_[config][0], 6); break;
379 case 40:/* 13.2 */
380 nsurf = add_triangle(tiling13_2_[config][1], 6); break;
381 case 41:/* 13.2 */
382 nsurf = add_triangle(tiling13_2_[config][2], 6); break;
383 case 42:/* 13.2 */
384 nsurf = add_triangle(tiling13_2_[config][3], 6); break;
385 case 43:/* 13.2 */
386 nsurf = add_triangle(tiling13_2_[config][4], 6); break;
387 case 44:/* 13.2 */
388 nsurf = add_triangle(tiling13_2_[config][5], 6); break;
389
390 case 45:/* 13.1 */
391 nsurf = add_triangle(tiling13_1_[config], 4); break;
392
393 default:
394 print_cube();
395 error->one(FLERR,"Marching cubes - impossible case 13");
396 }
397 break;
398
399 case 14:
400 nsurf = add_triangle(tiling14[config], 4);
401 break;
402 };
403
404 // store 4 MC labels for FixAblate caller
405
406 mcflags[icell][0] = icase;
407 mcflags[icell][1] = config;
408 mcflags[icell][2] = subconfig;
409 mcflags[icell][3] = nsurf;
410
411 // populate Grid and Surf data structs
412 // points will be duplicated, not unique
413 // surf ID = cell ID for all surfs in cell
414
415 ptr = csurfs->get(nsurf);
416
417 ipt = 0;
418 for (i = 0; i < nsurf; i++) {
419 if (svalues) surf->add_tri(cells[icell].id,svalues[icell],
420 pt[ipt+2],pt[ipt+1],pt[ipt]);
421 else surf->add_tri(cells[icell].id,1,pt[ipt+2],pt[ipt+1],pt[ipt]);
422 ipt += 3;
423 isurf = surf->nlocal - 1;
424 ptr[i] = isurf;
425 }
426
427 cells[icell].nsurf = nsurf;
428 if (nsurf) {
429 cells[icell].csurfs = ptr;
430 cinfo[icell].type = OVERLAP;
431 }
432 }
433 }
434
435 /* ----------------------------------------------------------------------
436 interpolate function used by both marching squares and cubes
437 lo/hi = coordinates of end points of edge of square
438 v0/v1 = values at lo/hi end points
439 value = interpolated coordinate for thresh value
440 ------------------------------------------------------------------------- */
441
interpolate(double v0,double v1,double lo,double hi)442 double MarchingCubes::interpolate(double v0, double v1, double lo, double hi)
443 {
444 double value = lo + (hi-lo)*(thresh-v0)/(v1-v0);
445 value = MAX(value,lo);
446 value = MIN(value,hi);
447 return value;
448 }
449
450 /* ----------------------------------------------------------------------
451 clean up issues that marching cubes occasionally generates
452 that cause problems for SPARTA
453 what MC does:
454 may generate 0 or 2 triangles on the face of a cell
455 the cell sharing the face may also generate 0 or 2 triangles
456 the normals for the 2 triangles may be into or out of the owning cell
457 what SPARTA needs:
458 let cell1 and cell2 be two cells that share a face
459 if cell1 has 2 tris on face and cell2 has none:
460 if norm is into cell1: keep them in cell1
461 if norm is into cell2: assign both tris to cell2
462 if both cell1 and cell2 have 2 tris on face: delete all 4 tris
463 algorithm to do this:
464 loop over all my cells with implicit tris:
465 count how many surfs on each face
466 loop over all my cells with implicit tris:
467 loop over faces with 2 tris:
468 if I own adjoining cell:
469 check its tally on shared face
470 reassign and/or delete triangles as necessary
471 if I do not own adjoining cell:
472 add 2 tris to send list for this proc
473 irregular comm of send list to nearby procs (share faces of my cells)
474 each proc loops over its recv list:
475 if my cell face has 2 tris: delete them
476 if my cell face has 0 tris: skip or add 2 tris depending on norm
477 ------------------------------------------------------------------------- */
478
cleanup()479 void MarchingCubes::cleanup()
480 {
481 int i,j,k,m,icell,iface,nsurf,idim,nflag,inwardnorm;
482 int ntri_other,othercell,otherface,otherproc,otherlocal,othernsurf;
483 surfint *oldcsurfs;
484 surfint *ptr,*csurfs_other;
485 cellint cellID;
486 double *lo,*hi;
487 double *norm;
488
489 Surf::Tri *tris = surf->tris;
490 Grid::ChildCell *cells = grid->cells;
491 MyPage<surfint> *csurfs = grid->csurfs;
492 int nglocal = grid->nlocal;
493
494 Surf::Tri *tlist = NULL;
495 int nlist = 0;
496 int maxlist = 0;
497
498 // DEBUG
499
500 //int nstotal;
501 //MPI_Allreduce(&surf->nlocal,&nstotal,1,MPI_INT,MPI_SUM,world);
502 //if (me == 0) printf("TOTAL TRI before count: %d\n",nstotal);
503
504 // END of DEBUG
505
506 // count # of tris on each face of every cell I own
507
508 int **nfacetri;
509 int ***facetris;
510 memory->create(nfacetri,nglocal,6,"readisurf:nfacetri");
511 memory->create(facetris,nglocal,6,2,"readisurf:facetris");
512
513 for (icell = 0; icell < nglocal; icell++) {
514 nfacetri[icell][0] = nfacetri[icell][1] = nfacetri[icell][2] =
515 nfacetri[icell][3] = nfacetri[icell][4] = nfacetri[icell][5] = 0;
516
517 if (cells[icell].nsplit <= 0) continue;
518 nsurf = cells[icell].nsurf;
519 if (nsurf == 0) continue;
520
521 lo = cells[icell].lo;
522 hi = cells[icell].hi;
523
524 for (j = 0; j < nsurf; j++) {
525 m = cells[icell].csurfs[j];
526 iface = Geometry::tri_on_hex_face(tris[m].p1,tris[m].p2,tris[m].p3,lo,hi);
527 if (iface < 0) continue;
528 if (nfacetri[icell][iface] < 2)
529 facetris[icell][iface][nfacetri[icell][iface]] = m;
530 nfacetri[icell][iface]++;
531 }
532 }
533
534 // check that every face has 0 or 2 tris
535
536 int flag = 0;
537 for (icell = 0; icell < nglocal; icell++)
538 for (iface = 0; iface < 6; iface++)
539 if (nfacetri[icell][iface] != 0 && nfacetri[icell][iface] != 2)
540 flag++;
541
542 int flagall;
543 MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
544 if (flagall)
545 error->all(FLERR,"Some cell faces do not have zero or 2 triangles");
546
547 // loop over all cell faces
548 // check tri count for that face for both adjoining cells
549
550 int *proclist = NULL;
551 SendDatum *bufsend = NULL;
552 int nsend = 0;
553 int maxsend = 0;
554
555 int *dellist = NULL;
556 int ndelete = 0;
557 int maxdelete = 0;
558
559 // DEBUG
560 //int ntotal = 0;
561 //int nadd = 0;
562 //int ndel = 0;
563
564 for (icell = 0; icell < nglocal; icell++) {
565 if (cells[icell].nsplit <= 0) continue;
566 nsurf = cells[icell].nsurf;
567 if (nsurf == 0) continue;
568
569 for (iface = 0; iface < 6; iface++) {
570 if (nfacetri[icell][iface] != 2) continue;
571 //ntotal += 2;
572
573 // other cell/face/proc = info for matching face in adjacent cell
574
575 nflag = grid->neigh_decode(cells[icell].nmask,iface);
576 if (nflag != NCHILD && nflag != NPBCHILD)
577 error->one(FLERR,"Invalid neighbor cell in cleanup_MC()");
578
579 norm = tris[facetris[icell][iface][0]].norm;
580 idim = iface/2;
581 if (iface % 2 && norm[idim] < 0.0) inwardnorm = 1;
582 else if (iface % 2 == 0 && norm[idim] > 0.0) inwardnorm = 1;
583 else inwardnorm = 0;
584 if (iface % 2) otherface = iface-1;
585 else otherface = iface+1;
586 othercell = (int) cells[icell].neigh[iface];
587 otherproc = cells[othercell].proc;
588 otherlocal = cells[othercell].ilocal;
589
590 // if I own the adjacent cell, make decision about shared tris
591 // if both cells have 2 tris on face, delete all of them
592 // otherwise cell that matches inward normal is assigned the 2 tris
593
594 if (otherproc == me) {
595 ntri_other = nfacetri[othercell][otherface];
596
597 // icell keeps the 2 tris
598
599 if (ntri_other == 0 && inwardnorm) continue;
600
601 // add 2 tris to othercell
602 // reset tri IDs to new owning cell
603
604 if (ntri_other == 0) {
605 othernsurf = cells[othercell].nsurf;
606 oldcsurfs = cells[othercell].csurfs;
607 ptr = csurfs->get(othernsurf+2);
608 for (k = 0; k < othernsurf; k++)
609 ptr[k] = oldcsurfs[k];
610 ptr[othernsurf] = facetris[icell][iface][0];
611 ptr[othernsurf+1] = facetris[icell][iface][1];
612 cells[othercell].nsurf += 2;
613 cells[othercell].csurfs = ptr;
614 tris[facetris[icell][iface][0]].id = cells[othercell].id;
615 tris[facetris[icell][iface][1]].id = cells[othercell].id;
616 //printf("MC add1 %d %d\n",cells[icell].id,cells[othercell].id);
617 //nadd += 2;
618 }
619
620 // delete 2 tris from othercell
621 // set nfacetri[othercell] = 0, so won't delete again when it is icell
622
623 if (ntri_other == 2) {
624 nfacetri[othercell][otherface] = 0;
625 othernsurf = cells[othercell].nsurf;
626 ptr = cells[othercell].csurfs;
627 m = facetris[othercell][otherface][0];
628 for (k = 0; k < othernsurf; k++)
629 if (ptr[k] == m) break;
630 if (k == othernsurf)
631 error->one(FLERR,"Could not find surf in cleanup_MC");
632 cells[othercell].csurfs[k] = cells[othercell].csurfs[othernsurf-1];
633 othernsurf--;
634 m = facetris[othercell][otherface][1];
635 for (k = 0; k < othernsurf; k++)
636 if (ptr[k] == m) break;
637 if (k == othernsurf)
638 error->one(FLERR,"Could not find surf in cleanup_MC");
639 cells[othercell].csurfs[k] = cells[othercell].csurfs[othernsurf-1];
640 othernsurf--;
641 cells[othercell].nsurf -= 2;
642 //printf("MC del1 %d %d\n",cells[icell].id,cells[othercell].id);
643 //ndel += 2;
644 }
645
646 // delete 2 tris from icell
647
648 ptr = cells[icell].csurfs;
649 m = facetris[icell][iface][0];
650 for (k = 0; k < nsurf; k++)
651 if (ptr[k] == m) break;
652 if (k == nsurf) error->one(FLERR,"Could not find surf in cleanup_MC");
653 cells[icell].csurfs[k] = cells[icell].csurfs[nsurf-1];
654 nsurf--;
655 m = facetris[icell][iface][1];
656 for (k = 0; k < nsurf; k++)
657 if (ptr[k] == m) break;
658 if (k == nsurf) error->one(FLERR,"Could not find surf in cleanup_MC");
659 cells[icell].csurfs[k] = cells[icell].csurfs[nsurf-1];
660 nsurf--;
661 cells[icell].nsurf -= 2;
662 //printf("MC dele %d %d\n",cells[icell].id,cells[othercell].id);
663 //ndel += 2;
664
665 // add 4 tris to delete list if both cells deleted them
666
667 if (ntri_other == 2) {
668 if (ndelete+4 > maxdelete) {
669 maxdelete += DELTA;
670 memory->grow(dellist,maxdelete,"readisurf:dellist");
671 }
672 dellist[ndelete++] = facetris[icell][iface][0];
673 dellist[ndelete++] = facetris[icell][iface][1];
674 dellist[ndelete++] = facetris[othercell][otherface][0];
675 dellist[ndelete++] = facetris[othercell][otherface][1];
676 }
677
678 // cell face is shared with another proc
679 // send it the cell/face indices and the 2 tris,
680 // in case they need to be assigned to the other cell based on norm
681
682 } else {
683 if (nsend == maxsend) {
684 maxsend += DELTA;
685 proclist = (int *)
686 memory->srealloc(proclist,maxsend*sizeof(int),
687 "readisurf:proclist");
688 bufsend = (SendDatum *)
689 memory->srealloc(bufsend,maxsend*sizeof(SendDatum),
690 "readisurf:bufsend");
691 }
692 proclist[nsend] = otherproc;
693 bufsend[nsend].sendcell = icell;
694 bufsend[nsend].sendface = iface;
695 bufsend[nsend].othercell = otherlocal;
696 bufsend[nsend].otherface = otherface;
697 bufsend[nsend].inwardnorm = inwardnorm;
698 memcpy(&bufsend[nsend].tri1,&tris[facetris[icell][iface][0]],
699 sizeof(Surf::Tri));
700 memcpy(&bufsend[nsend].tri2,&tris[facetris[icell][iface][1]],
701 sizeof(Surf::Tri));
702 nsend++;
703
704 // if not inwardnorm, delete 2 tris from this cell
705 // also add them to delete list
706
707 if (!inwardnorm) {
708 ptr = cells[icell].csurfs;
709 m = facetris[icell][iface][0];
710 for (k = 0; k < nsurf; k++)
711 if (ptr[k] == m) break;
712 if (k == nsurf) error->one(FLERR,"Could not find surf in cleanup_MC");
713 cells[icell].csurfs[k] = cells[icell].csurfs[nsurf-1];
714 nsurf--;
715 m = facetris[icell][iface][1];
716 for (k = 0; k < nsurf; k++)
717 if (ptr[k] == m) break;
718 if (k == nsurf) error->one(FLERR,"Could not find surf in cleanup_MC");
719 cells[icell].csurfs[k] = cells[icell].csurfs[nsurf-1];
720 nsurf--;
721 cells[icell].nsurf -= 2;
722 //ndel += 2;
723
724 if (ndelete+2 > maxdelete) {
725 maxdelete += DELTA;
726 memory->grow(dellist,maxdelete,"readisurf:dellist");
727 }
728 dellist[ndelete++] = facetris[icell][iface][0];
729 dellist[ndelete++] = facetris[icell][iface][1];
730 }
731 }
732 }
733 }
734
735 // perform irregular communication of list of cell faces and tri pairs
736
737 Irregular *irregular = new Irregular(sparta);
738 int nrecv = irregular->create_data_uniform(nsend,proclist,1);
739
740 SendDatum *bufrecv = (SendDatum *)
741 memory->smalloc(nrecv*sizeof(SendDatum),"readisurf:bufrecv");
742
743 irregular->exchange_uniform((char *) bufsend,sizeof(SendDatum),
744 (char *) bufrecv);
745 delete irregular;
746 memory->sfree(proclist);
747 memory->sfree(bufsend);
748
749 // loop over list of received face/tri info
750 // if my matching face has 2 tris, delete them
751 // if my matching face has 0 tris, skip or add 2 tris depending on norm
752
753 for (i = 0; i < nrecv; i++) {
754 icell = bufrecv[i].othercell;
755 iface = bufrecv[i].otherface;
756
757 // my icell is not affected, sender cell keeps its 2 tris
758
759 if (nfacetri[icell][iface] == 0 && bufrecv[i].inwardnorm) continue;
760
761 // add 2 tris to icell and this processor's Surf::tris list
762 // set tri IDs to new owning cell, must be done after memcpy()
763 // NOTE: what about tri types?
764
765 if (nfacetri[icell][iface] == 0) {
766 int nslocal = surf->nlocal;
767 surf->add_tri(cells[icell].id,1,
768 bufrecv[i].tri1.p1,bufrecv[i].tri1.p2,bufrecv[i].tri1.p3);
769 memcpy(&surf->tris[nslocal],&bufrecv[i].tri1,sizeof(Surf::Tri));
770 surf->tris[nslocal].id = cells[icell].id;
771 surf->add_tri(cells[icell].id,1,
772 bufrecv[i].tri2.p1,bufrecv[i].tri2.p2,bufrecv[i].tri2.p3);
773 memcpy(&surf->tris[nslocal+1],&bufrecv[i].tri2,sizeof(Surf::Tri));
774 surf->tris[nslocal+1].id = cells[icell].id;
775
776 nsurf = cells[icell].nsurf;
777 oldcsurfs = cells[icell].csurfs;
778 ptr = csurfs->get(nsurf+2);
779 for (k = 0; k < nsurf; k++)
780 ptr[k] = oldcsurfs[k];
781 ptr[nsurf] = nslocal;
782 ptr[nsurf+1] = nslocal+1;
783 cells[icell].nsurf += 2;
784 cells[icell].csurfs = ptr;
785 //nadd += 2;
786 }
787
788 // both cells have 2 tris on common face
789 // need to delete my 2 tris from icell
790 // sender will get similar message from me and delete
791 // inwardnorm check to see if I already deleted when sent a message,
792 // else delete now and add 2 tris to delete list
793
794 if (nfacetri[icell][iface] == 2) {
795 norm = tris[facetris[icell][iface][0]].norm;
796 idim = iface/2;
797 if (iface % 2 && norm[idim] < 0.0) inwardnorm = 1;
798 else if (iface % 2 == 0 && norm[idim] > 0.0) inwardnorm = 1;
799 else inwardnorm = 0;
800 if (!inwardnorm) continue;
801
802 nsurf = cells[icell].nsurf;
803 ptr = cells[icell].csurfs;
804 m = facetris[icell][iface][0];
805 for (k = 0; k < nsurf; k++)
806 if (ptr[k] == m) break;
807 if (k == nsurf) error->one(FLERR,"Could not find surf in cleanup_MC");
808 cells[icell].csurfs[k] = cells[icell].csurfs[nsurf-1];
809 nsurf--;
810 m = facetris[icell][iface][1];
811 for (k = 0; k < nsurf; k++)
812 if (ptr[k] == m) break;
813 if (k == nsurf) error->one(FLERR,"Could not find surf in cleanup_MC");
814 cells[icell].csurfs[k] = cells[icell].csurfs[nsurf-1];
815 nsurf--;
816 cells[icell].nsurf -= 2;
817 //ndel += 2;
818
819 if (ndelete+2 > maxdelete) {
820 maxdelete += DELTA;
821 memory->grow(dellist,maxdelete,"readisurf:dellist");
822 }
823 dellist[ndelete++] = facetris[icell][iface][0];
824 dellist[ndelete++] = facetris[icell][iface][1];
825 }
826 }
827
828 memory->sfree(bufrecv);
829 memory->destroy(nfacetri);
830 memory->destroy(facetris);
831
832 // compress Surf::tris list to remove deleted tris
833 // must sort dellist, so as to compress tris in DESCENDING index order
834 // descending, not ascending, so that a surf is not moved from end-of-list
835 // that is flagged for later deletion
836 // must repoint one location in cells->csurfs to moved surf
837 // requires grid hash to find owning cell of moved surf
838 // note that ghost surfs exist at this point, but caller will clear them
839
840 if (!grid->hashfilled) grid->rehash();
841
842 qsort(dellist,ndelete,sizeof(int),compare_indices);
843
844 tris = surf->tris;
845 int nslocal = surf->nlocal;
846 cellint celllID;
847 for (i = 0; i < ndelete; i++) {
848 m = dellist[i];
849 if (m != nslocal-1) memcpy(&tris[m],&tris[nslocal-1],sizeof(Surf::Tri));
850 nslocal--;
851
852 icell = (*grid->hash)[tris[m].id];
853 nsurf = cells[icell].nsurf;
854 ptr = cells[icell].csurfs;
855 for (k = 0; k < nsurf; k++)
856 if (ptr[k] == nslocal) {
857 ptr[k] = m;
858 break;
859 }
860 if (k == nsurf) error->one(FLERR,"Did not find moved tri in cleanup_MC()");
861 }
862
863 surf->nlocal = nslocal;
864 memory->destroy(dellist);
865
866 // DEBUG
867
868 /*
869 MPI_Allreduce(&surf->nlocal,&nstotal,1,MPI_INT,MPI_SUM,world);
870 if (me == 0) printf("TOTAL TRI after count: %d\n",nstotal);
871
872 int alltotal,alladd,alldel,allsend,allrecv;
873 MPI_Allreduce(&ntotal,&alltotal,1,MPI_INT,MPI_SUM,world);
874 MPI_Allreduce(&nadd,&alladd,1,MPI_INT,MPI_SUM,world);
875 MPI_Allreduce(&ndel,&alldel,1,MPI_INT,MPI_SUM,world);
876 MPI_Allreduce(&nsend,&allsend,1,MPI_INT,MPI_SUM,world);
877 MPI_Allreduce(&nrecv,&allrecv,1,MPI_INT,MPI_SUM,world);
878 if (me == 0)
879 printf("CLEANUP counts: total %d add %d del %d send %d recv %d\n",
880 alltotal,alladd,alldel,allsend,allrecv);
881
882 ntotal = 0;
883 int nbad = 0;
884 int nonface = 0;
885
886 for (icell = 0; icell < nglocal; icell++) {
887 if (cells[icell].nsplit <= 0) continue;
888 nsurf = cells[icell].nsurf;
889 if (nsurf == 0) continue;
890 ntotal += nsurf;
891
892 lo = cells[icell].lo;
893 hi = cells[icell].hi;
894
895 for (j = 0; j < nsurf; j++) {
896 m = cells[icell].csurfs[j];
897 iface = Geometry::tri_on_hex_face(tris[m].p1,tris[m].p2,tris[m].p3,lo,hi);
898 if (iface < 0) continue;
899
900 norm = tris[m].norm;
901 idim = iface/2;
902 if (iface % 2 && norm[idim] < 0.0) inwardnorm = 1;
903 else if (iface % 2 == 0 && norm[idim] > 0.0) inwardnorm = 1;
904 else inwardnorm = 0;
905
906 nonface++;
907 if (!inwardnorm) nbad++;
908 }
909 }
910
911 int nbadall;
912 MPI_Allreduce(&nbad,&nbadall,1,MPI_INT,MPI_SUM,world);
913 if (me == 0) printf("BAD NORM %d\n",nbadall);
914
915 int nonfaceall;
916 MPI_Allreduce(&nonface,&nonfaceall,1,MPI_INT,MPI_SUM,world);
917 if (me == 0) printf("Total onface %d\n",nonfaceall);
918
919 if (ntotal != surf->nlocal) error->one(FLERR,"Bad surf total");
920 */
921
922 // END of DEBUG
923 }
924
925 /* ----------------------------------------------------------------------
926 adding triangles
927 ------------------------------------------------------------------------- */
928
add_triangle(int * trig,int n)929 int MarchingCubes::add_triangle(int *trig, int n)
930 {
931 for(int t = 0; t < 3*n; t++) {
932 switch (trig[t]) {
933 case 0:
934 pt[t][0] = interpolate(v000,v001,lo[0],hi[0]);
935 pt[t][1] = lo[1];
936 pt[t][2] = lo[2];
937 break;
938 case 1:
939 pt[t][0] = hi[0];
940 pt[t][1] = interpolate(v001,v011,lo[1],hi[1]);
941 pt[t][2] = lo[2];
942 break;
943 case 2:
944 pt[t][0] = interpolate(v010,v011,lo[0],hi[0]);
945 pt[t][1] = hi[1];
946 pt[t][2] = lo[2];
947 break;
948 case 3:
949 pt[t][0] = lo[0];
950 pt[t][1] = interpolate(v000,v010,lo[1],hi[1]);
951 pt[t][2] = lo[2];
952 break;
953 case 4:
954 pt[t][0] = interpolate(v100,v101,lo[0],hi[0]);
955 pt[t][1] = lo[1];
956 pt[t][2] = hi[2];
957 break;
958 case 5:
959 pt[t][0] = hi[0];
960 pt[t][1] = interpolate(v101,v111,lo[1],hi[1]);
961 pt[t][2] = hi[2];
962 break;
963 case 6:
964 pt[t][0] = interpolate(v110,v111,lo[0],hi[0]);
965 pt[t][1] = hi[1];
966 pt[t][2] = hi[2];
967 break;
968 case 7:
969 pt[t][0] = lo[0];
970 pt[t][1] = interpolate(v100,v110,lo[1],hi[1]);
971 pt[t][2] = hi[2];
972 break;
973 case 8:
974 pt[t][0] = lo[0];
975 pt[t][1] = lo[1];
976 pt[t][2] = interpolate(v000,v100,lo[2],hi[2]);
977 break;
978 case 9:
979 pt[t][0] = hi[0];
980 pt[t][1] = lo[1];
981 pt[t][2] = interpolate(v001,v101,lo[2],hi[2]);
982 break;
983 case 10:
984 pt[t][0] = hi[0];
985 pt[t][1] = hi[1];
986 pt[t][2] = interpolate(v011,v111,lo[2],hi[2]);
987 break;
988 case 11:
989 pt[t][0] = lo[0];
990 pt[t][1] = hi[1];
991 pt[t][2] = interpolate(v010,v110,lo[2],hi[2]);
992 break;
993 case 12: {
994 int u = 0;
995 pt[t][0] = pt[t][1] = pt[t][2] = 0.0;
996 if (bit0 ^ bit1) {
997 ++u;
998 pt[t][0] += interpolate(v000,v001,lo[0],hi[0]);
999 pt[t][1] += lo[1];
1000 pt[t][2] += lo[2];
1001 }
1002 if (bit1 ^ bit2) {
1003 ++u;
1004 pt[t][0] += hi[0];
1005 pt[t][1] += interpolate(v001,v011,lo[1],hi[1]);
1006 pt[t][2] += lo[2];
1007 }
1008 if (bit2 ^ bit3) {
1009 ++u;
1010 pt[t][0] += interpolate(v010,v011,lo[0],hi[0]);
1011 pt[t][1] += hi[1];
1012 pt[t][2] += lo[2];
1013 }
1014 if (bit3 ^ bit0) {
1015 ++u;
1016 pt[t][0] += lo[0];
1017 pt[t][1] += interpolate(v000,v010,lo[1],hi[1]);
1018 pt[t][2] += lo[2];
1019 }
1020 if (bit4 ^ bit5) {
1021 ++u;
1022 pt[t][0] += interpolate(v100,v101,lo[0],hi[0]);
1023 pt[t][1] += lo[1];
1024 pt[t][2] += hi[2];
1025 }
1026 if (bit5 ^ bit6) {
1027 ++u;
1028 pt[t][0] += hi[0];
1029 pt[t][1] += interpolate(v101,v111,lo[1],hi[1]);
1030 pt[t][2] += hi[2];
1031 }
1032 if (bit6 ^ bit7) {
1033 ++u;
1034 pt[t][0] += interpolate(v110,v111,lo[0],hi[0]);
1035 pt[t][1] += hi[1];
1036 pt[t][2] += hi[2];
1037 }
1038 if (bit7 ^ bit4) {
1039 ++u;
1040 pt[t][0] += lo[0];
1041 pt[t][1] += interpolate(v100,v110,lo[1],hi[1]);
1042 pt[t][2] += hi[2];
1043 }
1044 if (bit0 ^ bit4) {
1045 ++u;
1046 pt[t][0] += lo[0];
1047 pt[t][1] += lo[1];
1048 pt[t][2] += interpolate(v000,v100,lo[2],hi[2]);
1049 }
1050 if (bit1 ^ bit5) {
1051 ++u;
1052 pt[t][0] += hi[0];
1053 pt[t][1] += lo[1];
1054 pt[t][2] += interpolate(v001,v101,lo[2],hi[2]);
1055 }
1056 if (bit2 ^ bit6) {
1057 ++u;
1058 pt[t][0] += hi[0];
1059 pt[t][1] += hi[1];
1060 pt[t][2] += interpolate(v011,v111,lo[2],hi[2]);
1061 }
1062 if (bit3 ^ bit7) {
1063 ++u;
1064 pt[t][0] += lo[0];
1065 pt[t][1] += hi[1];
1066 pt[t][2] += interpolate(v010,v110,lo[2],hi[2]);
1067 }
1068
1069 pt[t][0] /= static_cast<double> (u);
1070 pt[t][1] /= static_cast<double> (u);
1071 pt[t][2] /= static_cast<double> (u);
1072 break;
1073 }
1074
1075 default:
1076 break;
1077 }
1078 }
1079
1080 return n;
1081 }
1082
1083 /* ----------------------------------------------------------------------
1084 test a face
1085 if face > 0 return true if the face contains a part of the surface
1086 ------------------------------------------------------------------------- */
1087
test_face(int face)1088 bool MarchingCubes::test_face(int face)
1089 {
1090 double A,B,C,D;
1091
1092 switch (face) {
1093 case -1:
1094 case 1:
1095 A = v000iso;
1096 B = v100iso;
1097 C = v101iso;
1098 D = v001iso;
1099 break;
1100 case -2:
1101 case 2:
1102 A = v001iso;
1103 B = v101iso;
1104 C = v111iso;
1105 D = v011iso;
1106 break;
1107 case -3:
1108 case 3:
1109 A = v011iso;
1110 B = v111iso;
1111 C = v110iso;
1112 D = v010iso;
1113 break;
1114 case -4:
1115 case 4:
1116 A = v010iso;
1117 B = v110iso;
1118 C = v100iso;
1119 D = v000iso;
1120 break;
1121 case -5:
1122 case 5:
1123 A = v000iso;
1124 B = v010iso;
1125 C = v011iso;
1126 D = v001iso;
1127 break;
1128 case -6:
1129 case 6:
1130 A = v100iso;
1131 B = v110iso;
1132 C = v111iso;
1133 D = v101iso;
1134 break;
1135
1136 default:
1137 A = B = C = D = 0.0;
1138 print_cube();
1139 error->one(FLERR,"Invalid face code");
1140 };
1141
1142 if (fabs(A*C - B*D) < EPSILON) return face >= 0;
1143 return face * A * (A*C - B*D) >= 0 ; // face and A invert signs
1144 }
1145
1146 /* ----------------------------------------------------------------------
1147 test the interior of a cube
1148 icase = case of the active cube in [0..15]
1149 if s == 7, return true if the interior is empty
1150 if s == -7, return false if the interior is empty
1151 ------------------------------------------------------------------------- */
1152
test_interior(int s,int icase)1153 bool MarchingCubes::test_interior(int s, int icase)
1154 {
1155 double t,a,b,At=0.0,Bt=0.0,Ct=0.0,Dt=0.0;
1156 int test = 0;
1157 int edge = -1; // reference edge of the triangulation
1158
1159 switch (icase) {
1160 case 4 :
1161 case 10 :
1162 a = ( v100iso - v000iso ) * ( v111iso - v011iso ) -
1163 ( v110iso - v010iso ) * ( v101iso - v001iso ) ;
1164 b = v011iso * ( v100iso - v000iso ) + v000iso * ( v111iso - v011iso ) -
1165 v001iso * ( v110iso - v010iso ) - v010iso * ( v101iso - v001iso ) ;
1166 t = - b / (2*a) ;
1167 if (t < 0 || t > 1) return s>0 ;
1168
1169 At = v000iso + ( v100iso - v000iso ) * t ;
1170 Bt = v010iso + ( v110iso - v010iso ) * t ;
1171 Ct = v011iso + ( v111iso - v011iso ) * t ;
1172 Dt = v001iso + ( v101iso - v001iso ) * t ;
1173 break ;
1174
1175 case 6 :
1176 case 7 :
1177 case 12 :
1178 case 13 :
1179 switch( icase ) {
1180 case 6 : edge = test6 [config][2] ; break ;
1181 case 7 : edge = test7 [config][4] ; break ;
1182 case 12 : edge = test12[config][3] ; break ;
1183 case 13 : edge = tiling13_5_1[config][subconfig][0] ; break ;
1184 }
1185 switch( edge ) {
1186 case 0 :
1187 t = v000iso / ( v000iso - v001iso ) ;
1188 At = 0.0 ;
1189 Bt = v010iso + ( v011iso - v010iso ) * t ;
1190 Ct = v110iso + ( v111iso - v110iso ) * t ;
1191 Dt = v100iso + ( v101iso - v100iso ) * t ;
1192 break ;
1193 case 1 :
1194 t = v001iso / ( v001iso - v011iso ) ;
1195 At = 0.0 ;
1196 Bt = v000iso + ( v010iso - v000iso ) * t ;
1197 Ct = v100iso + ( v110iso - v100iso ) * t ;
1198 Dt = v101iso + ( v111iso - v101iso ) * t ;
1199 break ;
1200 case 2 :
1201 t = v011iso / ( v011iso - v010iso ) ;
1202 At = 0.0 ;
1203 Bt = v001iso + ( v000iso - v001iso ) * t ;
1204 Ct = v101iso + ( v100iso - v101iso ) * t ;
1205 Dt = v111iso + ( v110iso - v111iso ) * t ;
1206 break ;
1207 case 3 :
1208 t = v010iso / ( v010iso - v000iso ) ;
1209 At = 0.0 ;
1210 Bt = v011iso + ( v001iso - v011iso ) * t ;
1211 Ct = v111iso + ( v101iso - v111iso ) * t ;
1212 Dt = v110iso + ( v100iso - v110iso ) * t ;
1213 break ;
1214 case 4 :
1215 t = v100iso / ( v100iso - v101iso ) ;
1216 At = 0.0 ;
1217 Bt = v110iso + ( v111iso - v110iso ) * t ;
1218 Ct = v010iso + ( v011iso - v010iso ) * t ;
1219 Dt = v000iso + ( v001iso - v000iso ) * t ;
1220 break ;
1221 case 5 :
1222 t = v101iso / ( v101iso - v111iso ) ;
1223 At = 0.0 ;
1224 Bt = v100iso + ( v110iso - v100iso ) * t ;
1225 Ct = v000iso + ( v010iso - v000iso ) * t ;
1226 Dt = v001iso + ( v011iso - v001iso ) * t ;
1227 break ;
1228 case 6 :
1229 t = v111iso / ( v111iso - v110iso ) ;
1230 At = 0.0 ;
1231 Bt = v101iso + ( v100iso - v101iso ) * t ;
1232 Ct = v001iso + ( v000iso - v001iso ) * t ;
1233 Dt = v011iso + ( v010iso - v011iso ) * t ;
1234 break ;
1235 case 7 :
1236 t = v110iso / ( v110iso - v100iso ) ;
1237 At = 0.0 ;
1238 Bt = v111iso + ( v101iso - v111iso ) * t ;
1239 Ct = v011iso + ( v001iso - v011iso ) * t ;
1240 Dt = v010iso + ( v000iso - v010iso ) * t ;
1241 break ;
1242 case 8 :
1243 t = v000iso / ( v000iso - v100iso ) ;
1244 At = 0.0 ;
1245 Bt = v010iso + ( v110iso - v010iso ) * t ;
1246 Ct = v011iso + ( v111iso - v011iso ) * t ;
1247 Dt = v001iso + ( v101iso - v001iso ) * t ;
1248 break ;
1249 case 9 :
1250 t = v001iso / ( v001iso - v101iso ) ;
1251 At = 0.0 ;
1252 Bt = v000iso + ( v100iso - v000iso ) * t ;
1253 Ct = v010iso + ( v110iso - v010iso ) * t ;
1254 Dt = v011iso + ( v111iso - v011iso ) * t ;
1255 break ;
1256 case 10 :
1257 t = v011iso / ( v011iso - v111iso ) ;
1258 At = 0.0 ;
1259 Bt = v001iso + ( v101iso - v001iso ) * t ;
1260 Ct = v000iso + ( v100iso - v000iso ) * t ;
1261 Dt = v010iso + ( v110iso - v010iso ) * t ;
1262 break ;
1263 case 11 :
1264 t = v010iso / ( v010iso - v110iso ) ;
1265 At = 0.0 ;
1266 Bt = v011iso + ( v111iso - v011iso ) * t ;
1267 Ct = v001iso + ( v101iso - v001iso ) * t ;
1268 Dt = v000iso + ( v100iso - v000iso ) * t ;
1269 break ;
1270
1271 default:
1272 print_cube();
1273 error->one(FLERR,"Marching cubes - invalid edge");
1274 }
1275 break;
1276
1277 default:
1278 print_cube();
1279 error->one(FLERR,"Marching cubes - invalid ambiguous case");
1280 }
1281
1282 if (At >= 0.0) test ++;
1283 if (Bt >= 0.0) test += 2;
1284 if (Ct >= 0.0) test += 4;
1285 if (Dt >= 0.0) test += 8;
1286 switch (test) {
1287 case 0: return s>0;
1288 case 1: return s>0;
1289 case 2: return s>0;
1290 case 3: return s>0;
1291 case 4: return s>0;
1292 case 5:
1293 if (At * Ct - Bt * Dt < EPSILON) return s>0;
1294 break;
1295 case 6: return s>0;
1296 case 7: return s<0;
1297 case 8: return s>0;
1298 case 9: return s>0;
1299 case 10:
1300 if (At * Ct - Bt * Dt >= EPSILON) return s>0;
1301 break;
1302 case 11: return s<0;
1303 case 12: return s>0;
1304 case 13: return s<0;
1305 case 14: return s<0;
1306 case 15: return s<0;
1307 }
1308
1309 return s<0;
1310 }
1311
1312 /* ---------------------------------------------------------------------- */
1313
modified_test_interior(int s,int icase)1314 bool MarchingCubes::modified_test_interior(int s, int icase)
1315 {
1316 int edge = -1;
1317 int amb_face;
1318
1319 int inter_amb = 0;
1320
1321 switch (icase) {
1322 case 4:
1323 amb_face = 1;
1324 edge = interior_ambiguity(amb_face, s);
1325 inter_amb += interior_ambiguity_verification(edge);
1326
1327 amb_face = 2;
1328 edge = interior_ambiguity(amb_face, s);
1329 inter_amb += interior_ambiguity_verification(edge);
1330
1331 amb_face = 5;
1332 edge = interior_ambiguity(amb_face, s);
1333 inter_amb += interior_ambiguity_verification(edge);
1334
1335 if (inter_amb == 0) return false;
1336 else return true;
1337 break;
1338
1339 case 6:
1340 amb_face = abs(test6[config][0]);
1341
1342 edge = interior_ambiguity(amb_face, s);
1343 inter_amb = interior_ambiguity_verification(edge);
1344
1345 if (inter_amb == 0) return false;
1346 else return true;
1347
1348 break;
1349
1350 case 7:
1351 s = s * -1;
1352
1353 amb_face = 1;
1354 edge = interior_ambiguity(amb_face, s);
1355 inter_amb += interior_ambiguity_verification(edge);
1356
1357 amb_face = 2;
1358 edge = interior_ambiguity(amb_face, s);
1359 inter_amb += interior_ambiguity_verification(edge);
1360
1361 amb_face = 5;
1362 edge = interior_ambiguity(amb_face, s);
1363 inter_amb += interior_ambiguity_verification(edge);
1364
1365 if (inter_amb == 0) return false;
1366 else return true;
1367 break;
1368
1369 case 10:
1370 amb_face = abs(test10[config][0]);
1371
1372 edge = interior_ambiguity(amb_face, s);
1373 inter_amb = interior_ambiguity_verification(edge);
1374
1375 if (inter_amb == 0) return false;
1376 else return true;
1377 break;
1378
1379 case 12:
1380 amb_face = abs(test12[config][0]);
1381 edge = interior_ambiguity(amb_face, s);
1382 inter_amb += interior_ambiguity_verification(edge);
1383
1384
1385 amb_face = abs(test12[config][1]);
1386 edge = interior_ambiguity(amb_face, s);
1387 inter_amb += interior_ambiguity_verification(edge);
1388
1389 if (inter_amb == 0) return false;
1390 else return true;
1391 break;
1392 }
1393
1394 // should never reach here
1395
1396 return true;
1397 }
1398
1399 /* ---------------------------------------------------------------------- */
1400
interior_ambiguity(int amb_face,int s)1401 int MarchingCubes::interior_ambiguity(int amb_face, int s)
1402 {
1403 int edge;
1404
1405 switch (amb_face) {
1406 case 1:
1407 case 3:
1408 if (((v001iso * s) > 0) && ((v110iso * s) > 0)) edge = 4;
1409 if (((v000iso * s) > 0) && ((v111iso * s) > 0)) edge = 5;
1410 if (((v010iso * s) > 0) && ((v101iso * s) > 0)) edge = 6;
1411 if (((v011iso * s) > 0) && ((v100iso * s) > 0)) edge = 7;
1412 break;
1413
1414 case 2:
1415 case 4:
1416 if (((v001iso * s) > 0) && ((v110iso * s) > 0)) edge = 0;
1417 if (((v011iso * s) > 0) && ((v100iso * s) > 0)) edge = 1;
1418 if (((v010iso * s) > 0) && ((v101iso * s) > 0)) edge = 2;
1419 if (((v000iso * s) > 0) && ((v111iso * s) > 0)) edge = 3;
1420 break;
1421
1422 case 5:
1423 case 6:
1424 case 0:
1425 if (((v000iso * s) > 0) && ((v111iso * s) > 0)) edge = 8;
1426 if (((v001iso * s) > 0) && ((v110iso * s) > 0)) edge = 9;
1427 if (((v011iso * s) > 0) && ((v100iso * s) > 0)) edge = 10;
1428 if (((v010iso * s) > 0) && ((v101iso * s) > 0)) edge = 11;
1429 break;
1430 }
1431
1432 return edge;
1433 }
1434
1435 /* ---------------------------------------------------------------------- */
1436
interior_ambiguity_verification(int edge)1437 int MarchingCubes::interior_ambiguity_verification(int edge)
1438 {
1439 double t, At = 0.0, Bt = 0.0, Ct = 0.0, Dt = 0.0, a = 0.0, b = 0.0;
1440 double verify;
1441
1442 switch (edge) {
1443
1444 case 0:
1445 a = (v000iso - v001iso) * (v110iso - v111iso)
1446 - (v100iso - v101iso) * (v010iso - v011iso);
1447 b = v111iso * (v000iso - v001iso) + v001iso * (v110iso - v111iso)
1448 - v011iso * (v100iso - v101iso)
1449 - v101iso * (v010iso - v011iso);
1450
1451 if (a > 0)
1452 return 1;
1453
1454 t = -b / (2 * a);
1455 if (t < 0 || t > 1)
1456 return 1;
1457
1458 At = v001iso + (v000iso - v001iso) * t;
1459 Bt = v101iso + (v100iso - v101iso) * t;
1460 Ct = v111iso + (v110iso - v111iso) * t;
1461 Dt = v011iso + (v010iso - v011iso) * t;
1462
1463 verify = At * Ct - Bt * Dt;
1464
1465 if (verify > 0)
1466 return 0;
1467 if (verify < 0)
1468 return 1;
1469
1470 break;
1471
1472 case 1:
1473 a = (v010iso - v011iso) * (v100iso - v101iso)
1474 - (v000iso - v001iso) * (v110iso - v111iso);
1475 b = v101iso * (v010iso - v011iso) + v011iso * (v100iso - v101iso)
1476 - v111iso * (v000iso - v001iso)
1477 - v001iso * (v110iso - v111iso);
1478
1479 if (a > 0)
1480 return 1;
1481
1482 t = -b / (2 * a);
1483 if (t < 0 || t > 1)
1484 return 1;
1485
1486 At = v011iso + (v010iso - v011iso) * t;
1487 Bt = v001iso + (v000iso - v001iso) * t;
1488 Ct = v101iso + (v100iso - v101iso) * t;
1489 Dt = v111iso + (v110iso - v111iso) * t;
1490
1491 verify = At * Ct - Bt * Dt;
1492
1493 if (verify > 0)
1494 return 0;
1495 if (verify < 0)
1496 return 1;
1497 break;
1498
1499 case 2:
1500 a = (v011iso - v010iso) * (v101iso - v100iso)
1501 - (v111iso - v110iso) * (v001iso - v000iso);
1502 b = v100iso * (v011iso - v010iso) + v010iso * (v101iso - v100iso)
1503 - v000iso * (v111iso - v110iso)
1504 - v110iso * (v001iso - v000iso);
1505 if (a > 0)
1506 return 1;
1507
1508 t = -b / (2 * a);
1509 if (t < 0 || t > 1)
1510 return 1;
1511
1512 At = v010iso + (v011iso - v010iso) * t;
1513 Bt = v110iso + (v111iso - v110iso) * t;
1514 Ct = v100iso + (v101iso - v100iso) * t;
1515 Dt = v000iso + (v001iso - v000iso) * t;
1516
1517 verify = At * Ct - Bt * Dt;
1518
1519 if (verify > 0)
1520 return 0;
1521 if (verify < 0)
1522 return 1;
1523 break;
1524
1525 case 3:
1526 a = (v001iso - v000iso) * (v111iso - v110iso)
1527 - (v011iso - v010iso) * (v101iso - v100iso);
1528 b = v110iso * (v001iso - v000iso) + v000iso * (v111iso - v110iso)
1529 - v100iso * (v011iso - v010iso)
1530 - v010iso * (v101iso - v100iso);
1531 if (a > 0)
1532 return 1;
1533
1534 t = -b / (2 * a);
1535 if (t < 0 || t > 1)
1536 return 1;
1537
1538 At = v000iso + (v001iso - v000iso) * t;
1539 Bt = v010iso + (v011iso - v010iso) * t;
1540 Ct = v110iso + (v111iso - v110iso) * t;
1541 Dt = v100iso + (v101iso - v100iso) * t;
1542
1543 verify = At * Ct - Bt * Dt;
1544
1545 if (verify > 0)
1546 return 0;
1547 if (verify < 0)
1548 return 1;
1549 break;
1550
1551 case 4:
1552
1553 a = (v011iso - v001iso) * (v110iso - v100iso)
1554 - (v010iso - v000iso) * (v111iso - v101iso);
1555 b = v100iso * (v011iso - v001iso) + v001iso * (v110iso - v100iso)
1556 - v101iso * (v010iso - v000iso)
1557 - v000iso * (v111iso - v101iso);
1558
1559 if (a > 0)
1560 return 1;
1561
1562 t = -b / (2 * a);
1563 if (t < 0 || t > 1)
1564 return 1;
1565
1566 At = v001iso + (v011iso - v001iso) * t;
1567 Bt = v000iso + (v010iso - v000iso) * t;
1568 Ct = v100iso + (v110iso - v100iso) * t;
1569 Dt = v101iso + (v111iso - v101iso) * t;
1570
1571 verify = At * Ct - Bt * Dt;
1572
1573 if (verify > 0)
1574 return 0;
1575 if (verify < 0)
1576 return 1;
1577 break;
1578
1579 case 5:
1580
1581 a = (v010iso - v000iso) * (v111iso - v101iso)
1582 - (v011iso - v001iso) * (v110iso - v100iso);
1583 b = v101iso * (v010iso - v000iso) + v000iso * (v111iso - v101iso)
1584 - v100iso * (v011iso - v001iso)
1585 - v001iso * (v110iso - v100iso);
1586 if (a > 0)
1587 return 1;
1588
1589 t = -b / (2 * a);
1590 if (t < 0 || t > 1)
1591 return 1;
1592
1593 At = v000iso + (v010iso - v000iso) * t;
1594 Bt = v001iso + (v011iso - v001iso) * t;
1595 Ct = v101iso + (v111iso - v101iso) * t;
1596 Dt = v100iso + (v110iso - v100iso) * t;
1597
1598 verify = At * Ct - Bt * Dt;
1599
1600 if (verify > 0)
1601 return 0;
1602 if (verify < 0)
1603 return 1;
1604 break;
1605
1606 case 6:
1607 a = (v000iso - v010iso) * (v101iso - v111iso)
1608 - (v100iso - v110iso) * (v001iso - v011iso);
1609 b = v111iso * (v000iso - v010iso) + v010iso * (v101iso - v111iso)
1610 - v011iso * (v100iso - v110iso)
1611 - v110iso * (v001iso - v011iso);
1612 if (a > 0)
1613 return 1;
1614
1615 t = -b / (2 * a);
1616 if (t < 0 || t > 1)
1617 return 1;
1618
1619 At = v010iso + (v000iso - v010iso) * t;
1620 Bt = v110iso + (v100iso - v110iso) * t;
1621 Ct = v111iso + (v101iso - v111iso) * t;
1622 Dt = v011iso + (v001iso - v011iso) * t;
1623
1624 verify = At * Ct - Bt * Dt;
1625
1626 if (verify > 0)
1627 return 0;
1628 if (verify < 0)
1629 return 1;
1630 break;
1631
1632 case 7:
1633 a = (v001iso - v011iso) * (v100iso - v110iso)
1634 - (v000iso - v010iso) * (v101iso - v111iso);
1635 b = v110iso * (v001iso - v011iso) + v011iso * (v100iso - v110iso)
1636 - v111iso * (v000iso - v010iso)
1637 - v010iso * (v101iso - v111iso);
1638 if (a > 0)
1639 return 1;
1640
1641 t = -b / (2 * a);
1642 if (t < 0 || t > 1)
1643 return 1;
1644
1645 At = v011iso + (v001iso - v011iso) * t;
1646 Bt = v010iso + (v000iso - v010iso) * t;
1647 Ct = v110iso + (v100iso - v110iso) * t;
1648 Dt = v111iso + (v101iso - v111iso) * t;
1649
1650 verify = At * Ct - Bt * Dt;
1651
1652 if (verify > 0)
1653 return 0;
1654 if (verify < 0)
1655 return 1;
1656 break;
1657
1658 case 8:
1659 a = (v100iso - v000iso) * (v111iso - v011iso)
1660 - (v110iso - v010iso) * (v101iso - v001iso);
1661 b = v011iso * (v100iso - v000iso) + v000iso * (v111iso - v011iso)
1662 - v001iso * (v110iso - v010iso)
1663 - v010iso * (v101iso - v001iso);
1664 if (a > 0)
1665 return 1;
1666
1667 t = -b / (2 * a);
1668 if (t < 0 || t > 1)
1669 return 1;
1670
1671 At = v000iso + (v100iso - v000iso) * t;
1672 Bt = v010iso + (v110iso - v010iso) * t;
1673 Ct = v011iso + (v111iso - v011iso) * t;
1674 Dt = v001iso + (v101iso - v001iso) * t;
1675
1676 verify = At * Ct - Bt * Dt;
1677
1678 if (verify > 0)
1679 return 0;
1680 if (verify < 0)
1681 return 1;
1682 break;
1683
1684 case 9:
1685 a = (v101iso - v001iso) * (v110iso - v010iso)
1686 - (v100iso - v000iso) * (v111iso - v011iso);
1687 b = v010iso * (v101iso - v001iso) + v001iso * (v110iso - v010iso)
1688 - v011iso * (v100iso - v000iso)
1689 - v000iso * (v111iso - v011iso);
1690 if (a > 0)
1691 return 1;
1692
1693 t = -b / (2 * a);
1694 if (t < 0 || t > 1)
1695 return 1;
1696
1697 At = v001iso + (v101iso - v001iso) * t;
1698 Bt = v000iso + (v100iso - v000iso) * t;
1699 Ct = v010iso + (v110iso - v010iso) * t;
1700 Dt = v011iso + (v111iso - v011iso) * t;
1701
1702 verify = At * Ct - Bt * Dt;
1703
1704 if (verify > 0)
1705 return 0;
1706 if (verify < 0)
1707 return 1;
1708 break;
1709
1710 case 10:
1711 a = (v111iso - v011iso) * (v100iso - v000iso)
1712 - (v101iso - v001iso) * (v110iso - v010iso);
1713 b = v000iso * (v111iso - v011iso) + v011iso * (v100iso - v000iso)
1714 - v010iso * (v101iso - v001iso)
1715 - v001iso * (v110iso - v010iso);
1716 if (a > 0)
1717 return 1;
1718
1719 t = -b / (2 * a);
1720 if (t < 0 || t > 1)
1721 return 1;
1722
1723 At = v011iso + (v111iso - v011iso) * t;
1724 Bt = v001iso + (v101iso - v001iso) * t;
1725 Ct = v000iso + (v100iso - v000iso) * t;
1726 Dt = v010iso + (v110iso - v010iso) * t;
1727
1728 verify = At * Ct - Bt * Dt;
1729
1730 if (verify > 0)
1731 return 0;
1732 if (verify < 0)
1733 return 1;
1734 break;
1735
1736 case 11:
1737 a = (v110iso - v010iso) * (v101iso - v001iso)
1738 - (v111iso - v011iso) * (v100iso - v000iso);
1739 b = v001iso * (v110iso - v010iso) + v010iso * (v101iso - v001iso)
1740 - v000iso * (v111iso - v011iso)
1741 - v011iso * (v100iso - v000iso);
1742 if (a > 0)
1743 return 1;
1744
1745 t = -b / (2 * a);
1746 if (t < 0 || t > 1)
1747 return 1;
1748
1749 At = v010iso + (v110iso - v010iso) * t;
1750 Bt = v011iso + (v111iso - v011iso) * t;
1751 Ct = v001iso + (v101iso - v001iso) * t;
1752 Dt = v000iso + (v100iso - v000iso) * t;
1753
1754 verify = At * Ct - Bt * Dt;
1755
1756 if (verify > 0)
1757 return 0;
1758 if (verify < 0)
1759 return 1;
1760 break;
1761 }
1762
1763 // should never reach here
1764
1765 return 1;
1766 }
1767
1768 /* ----------------------------------------------------------------------
1769 return true if the interior is empty (two faces)
1770 ------------------------------------------------------------------------- */
1771
interior_test_case13()1772 bool MarchingCubes::interior_test_case13()
1773 {
1774 double t1, t2, At1 = 0.0, Bt1 = 0.0, Ct1 = 0.0, Dt1 = 0.0;
1775 double At2 = 0.0, Bt2 = 0.0, Ct2 = 0.0, Dt2 = 0.0, a = 0.0, b = 0.0, c = 0.0;
1776
1777 a = (v000iso - v001iso) * (v110iso - v111iso)
1778 - (v100iso - v101iso) * (v010iso - v011iso);
1779 b = v111iso * (v000iso - v001iso) + v001iso * (v110iso - v111iso)
1780 - v011iso * (v100iso - v101iso)
1781 - v101iso * (v010iso - v011iso);
1782 c = v001iso*v111iso - v101iso*v011iso;
1783
1784 double delta = b*b - 4*a*c;
1785
1786 t1 = (-b + sqrt(delta))/(2*a);
1787 t2 = (-b - sqrt(delta))/(2*a);
1788
1789 // DEBUG
1790 // printf("delta = %f, t1 = %f, t2 = %f\n", delta, t1, t2);
1791
1792 if ((t1 < 1)&&(t1>0) &&(t2 < 1)&&(t2 > 0)) {
1793 At1 = v001iso + (v000iso - v001iso) * t1;
1794 Bt1 = v101iso + (v100iso - v101iso) * t1;
1795 Ct1 = v111iso + (v110iso - v111iso) * t1;
1796 Dt1 = v011iso + (v010iso - v011iso) * t1;
1797
1798 double x1 = (At1 - Dt1)/(At1 + Ct1 - Bt1 - Dt1);
1799 double y1 = (At1 - Bt1)/(At1 + Ct1 - Bt1 - Dt1);
1800
1801 At2 = v001iso + (v000iso - v001iso) * t2;
1802 Bt2 = v101iso + (v100iso - v101iso) * t2;
1803 Ct2 = v111iso + (v110iso - v111iso) * t2;
1804 Dt2 = v011iso + (v010iso - v011iso) * t2;
1805
1806 double x2 = (At2 - Dt2)/(At2 + Ct2 - Bt2 - Dt2);
1807 double y2 = (At2 - Bt2)/(At2 + Ct2 - Bt2 - Dt2);
1808
1809 if ((x1 < 1)&&(x1>0) &&(x2 < 1)&&(x2 > 0) &&
1810 (y1 < 1)&&(y1>0) &&(y2 < 1)&&(y2 > 0)) return false;
1811 }
1812
1813 return true;
1814 }
1815
1816 /* ----------------------------------------------------------------------
1817 comparison function invoked by qsort() called by cleanup()
1818 used to sort the dellist of removed tris into DESCENDING order
1819 this is not a class method
1820 ------------------------------------------------------------------------- */
1821
compare_indices(const void * iptr,const void * jptr)1822 int compare_indices(const void *iptr, const void *jptr)
1823 {
1824 int i = *((int *) iptr);
1825 int j = *((int *) jptr);
1826 if (i < j) return 1;
1827 if (i > j) return -1;
1828 return 0;
1829 }
1830
1831 /* ----------------------------------------------------------------------
1832 print cube for debugging
1833 ------------------------------------------------------------------------- */
1834
print_cube()1835 void MarchingCubes::print_cube()
1836 {
1837 fprintf(screen,"\t %d %d %d %d %d %d %d %d\n",
1838 v000,v001,v011,v010,v100,v101,v111,v110);
1839 }
1840