1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <string.h>
4 #include <ctype.h>
5 #include <math.h>
6 #ifdef _WIN32
7 # define WIN32_LEAN_AND_MEAN
8 # include <windows.h>
9 #endif
10 #include "gl_config.h"
11
12 #include "tk.h"
13 #include "cgnslib.h"
14 #include "hash.h"
15
16 /* define this to exclude structured mesh boundaries */
17 /* imin, imax, jmin, jmax, kmin, kmax */
18
19 /*#define NO_MESH_BOUNDARIES*/
20
21 /* define this to disable a cutting plane */
22 /* the cutting plane requires more memory since all
23 nodes and elements are saved */
24
25 /*#define NO_CUTTING_PLANE*/
26
27 /* this is the cosine of the angle between faces to define an interior edge */
28 /* it's set to 60 degrees - comment to remove interior edges */
29
30 #define EDGE_ANGLE 0.5
31
32 typedef float Node[3];
33
34 typedef struct {
35 cgsize_t id;
36 cgsize_t nodes[2];
37 } Edge;
38
39 typedef struct {
40 cgsize_t id;
41 int flags;
42 int nnodes;
43 cgsize_t *nodes;
44 float normal[3];
45 } Face;
46
47 typedef struct {
48 Face *face;
49 cgsize_t num;
50 int flags;
51 } PolyFace;
52
53 #ifndef NO_CUTTING_PLANE
54
55 typedef struct {
56 int id;
57 cgsize_t nodes[3];
58 float ratio;
59 } CutNode;
60
61 typedef struct {
62 cgsize_t nelems;
63 cgsize_t *elems;
64 cgsize_t nedges;
65 Edge *edges;
66 } CutData;
67
68 typedef struct {
69 float plane[4];
70 cgsize_t nelems;
71 cgsize_t nedges;
72 cgsize_t nnodes;
73 Node *nodes;
74 } CutPlane;
75
76 static float cutcolor[4] = {(float)0.8, (float)0.4, (float)0.8, (float)0.5};
77 static int usecutclr = 0;
78 static int ignorevis = 0;
79 static CutPlane cutplane;
80 static int CutDL = 0;
81 static int PlaneDL = 0;
82
83 #endif
84
85 typedef struct {
86 char name[33];
87 int type;
88 int dim;
89 cgsize_t data[10];
90 char d_name[33];
91 cgsize_t nedges;
92 Edge *edges;
93 cgsize_t nfaces;
94 Face **faces;
95 #ifndef NO_CUTTING_PLANE
96 CGNS_ENUMT(ElementType_t) elemtype;
97 cgsize_t nelems;
98 cgsize_t *elems;
99 cgsize_t *elem_offsets;
100 int npoly;
101 Face **poly;
102 CutData cut;
103 #endif
104 float bbox[3][2];
105 int dlist;
106 int mode;
107 float color[4];
108 char errmsg[81];
109 } Regn;
110
111 typedef struct {
112 char name[33];
113 cgsize_t nnodes;
114 Node *nodes;
115 int nregs;
116 Regn *regs;
117 } Zone;
118
119 static int cgnsfn = 0;
120 static int cgnsbase = 0;
121 static int cgnszone = 0;
122 static int nbases = 0;
123 static int nzones = 0;
124 static Zone *zones;
125 static int AxisDL = 0;
126
127 static char BaseName[33];
128 static int CellDim, PhyDim;
129
130 static Tcl_Interp *global_interp;
131
132 enum {
133 REG_MESH,
134 REG_ELEM,
135 REG_1TO1,
136 REG_CONN,
137 REG_HOLE,
138 REG_BOCO,
139 REG_BNDS
140 };
141
142 /* mapping from elements to faces */
143
144 static int facenodes[22][5] = {
145 /* tri */
146 {3, 0, 1, 2, 0},
147 /* quad */
148 {4, 0, 1, 2, 3},
149 /* tet */
150 {3, 0, 2, 1, 0},
151 {3, 0, 1, 3, 0},
152 {3, 1, 2, 3, 0},
153 {3, 2, 0, 3, 0},
154 /* pyramid */
155 {4, 0, 3, 2, 1},
156 {3, 0, 1, 4, 0},
157 {3, 1, 2, 4, 0},
158 {3, 2, 3, 4, 0},
159 {3, 3, 0, 4, 0},
160 /* wedge */
161 {4, 0, 1, 4, 3},
162 {4, 1, 2, 5, 4},
163 {4, 2, 0, 3, 5},
164 {3, 0, 2, 1, 0},
165 {3, 3, 4, 5, 0},
166 /* hex */
167 {4, 0, 3, 2, 1},
168 {4, 0, 1, 5, 4},
169 {4, 1, 2, 6, 5},
170 {4, 2, 3, 7, 6},
171 {4, 0, 4, 7, 3},
172 {4, 4, 5, 6, 7}
173 };
174
175 /*===================================================================
176 * local routines
177 *===================================================================*/
178
179 /*-------------------------------------------------------------------*/
180
FATAL(char * errmsg)181 static void FATAL (char *errmsg)
182 {
183 char cmd[129];
184
185 sprintf (cmd, "error_exit {%s}", errmsg);
186 Tcl_Eval (global_interp, cmd);
187 exit (1);
188 }
189
190 /*-------------------------------------------------------------------*/
191
MALLOC(char * funcname,size_t bytes)192 static void *MALLOC (char *funcname, size_t bytes)
193 {
194 void *data = calloc (bytes, 1);
195 if (NULL == data) {
196 char msg[128];
197 if (funcname != NULL)
198 sprintf (msg, "%s:malloc failed for %lu bytes", funcname,
199 (unsigned long)bytes);
200 else
201 sprintf (msg, "malloc failed for %lu bytes",
202 (unsigned long)bytes);
203 FATAL (msg);
204 }
205 return data;
206 }
207
208 /*-------------------------------------------------------------------*/
209
REALLOC(char * funcname,size_t bytes,void * old_data)210 static void *REALLOC (char *funcname, size_t bytes, void *old_data)
211 {
212 void *new_data;
213
214 if (NULL == old_data)
215 return MALLOC (funcname, bytes);
216 new_data = realloc (old_data, bytes);
217 if (NULL == new_data) {
218 char msg[128];
219 if (funcname != NULL)
220 sprintf (msg, "%s:realloc failed for %lu bytes", funcname,
221 (unsigned long)bytes);
222 else
223 sprintf (msg, "realloc failed for %lu bytes",
224 (unsigned long)bytes);
225 FATAL (msg);
226 }
227 return new_data;
228 }
229
230 /*-------------------------------------------------------------------*/
231
new_face(char * funcname,int nnodes)232 static Face *new_face(char *funcname, int nnodes)
233 {
234 Face *f = (Face *) calloc (nnodes * sizeof(cgsize_t) + sizeof(Face), 1);
235 if (f == NULL) {
236 char msg[128];
237 if (funcname != NULL)
238 sprintf (msg, "%s:malloc failed for face with %d nodes",
239 funcname, nnodes);
240 else
241 sprintf (msg, "malloc failed for face with %d nodes",
242 nnodes);
243 FATAL (msg);
244 }
245 f->nnodes = nnodes;
246 f->nodes = (cgsize_t *)(f + 1);
247 return f;
248 }
249
250 /*-------------------------------------------------------------------*/
251
copy_face(char * funcname,Face * face)252 static Face *copy_face(char *funcname, Face *face)
253 {
254 int n;
255 Face *f = new_face(funcname, face->nnodes);
256
257 f->id = face->id;
258 f->flags = face->flags;
259 for (n = 0; n < face->nnodes; n++)
260 f->nodes[n] = face->nodes[n];
261 for (n = 0; n < 3; n++)
262 f->normal[n] = face->normal[n];
263 return f;
264 }
265
266 /*-------------------------------------------------------------------*/
267
zone_message(char * msg,char * name)268 static void zone_message (char *msg, char *name)
269 {
270 char cmd[129];
271
272 if (name != NULL)
273 sprintf (cmd, "display_message {Zone %d : %s %s...}",
274 cgnszone, msg, name);
275 else
276 sprintf (cmd, "display_message {Zone %d : %s...}", cgnszone, msg);
277 Tcl_Eval (global_interp, cmd);
278 }
279
280 /*-------------------------------------------------------------------*/
281
free_all(void)282 static void free_all (void)
283 {
284 int nz, nr, nf;
285
286 if (!nzones) return;
287 for (nz = 0; nz < nzones; nz++) {
288 for (nr = 0; nr < zones[nz].nregs; nr++) {
289 if (zones[nz].regs[nr].dlist)
290 glDeleteLists (zones[nz].regs[nr].dlist, 1);
291 if (zones[nz].regs[nr].nedges)
292 free (zones[nz].regs[nr].edges);
293 if (zones[nz].regs[nr].nfaces) {
294 for (nf = 0; nf < zones[nz].regs[nr].nfaces; nf++) {
295 if (zones[nz].regs[nr].faces[nf])
296 free (zones[nz].regs[nr].faces[nf]);
297 }
298 free (zones[nz].regs[nr].faces);
299 }
300 #ifndef NO_CUTTING_PLANE
301 if (zones[nz].regs[nr].nelems)
302 free (zones[nz].regs[nr].elems);
303 if (zones[nz].regs[nr].nelems && zones[nz].regs[nr].elem_offsets)
304 free (zones[nz].regs[nr].elem_offsets);
305 if (zones[nz].regs[nr].npoly)
306 free (zones[nz].regs[nr].poly);
307 if (zones[nz].regs[nr].cut.nelems)
308 free (zones[nz].regs[nr].cut.elems);
309 if (zones[nz].regs[nr].cut.nedges)
310 free (zones[nz].regs[nr].cut.edges);
311 #endif
312 }
313 if (zones[nz].nregs) free (zones[nz].regs);
314 if (zones[nz].nnodes) free(zones[nz].nodes);
315 }
316 free (zones);
317 nzones = 0;
318 #ifndef NO_CUTTING_PLANE
319 if (cutplane.nnodes)
320 free (cutplane.nodes);
321 cutplane.nelems = 0;
322 cutplane.nedges = 0;
323 cutplane.nnodes = 0;
324 #endif
325 }
326
327 /*-------------------------------------------------------------------*/
328
compare_ints(const void * v1,const void * v2)329 static int compare_ints (const void *v1, const void *v2)
330 {
331 return (int)(*((cgsize_t *)v1) - *((cgsize_t *)v2));
332 }
333
334 /*-------------------------------------------------------------------*/
335
find_int(cgsize_t value,cgsize_t nlist,cgsize_t * list)336 static int find_int (cgsize_t value, cgsize_t nlist, cgsize_t *list)
337 {
338 cgsize_t lo = 0, hi = nlist - 1, mid;
339
340 if (!nlist || value < list[0]) return 0;
341 if (value == list[0]) return 1;
342 if (!hi || value > list[hi]) return 0;
343 if (value == list[hi]) return 1;
344
345 while (lo <= hi) {
346 mid = (lo + hi) >> 1;
347 if (value == list[mid]) return 1;
348 if (value < list[mid])
349 hi = mid - 1;
350 else
351 lo = mid + 1;
352 }
353 return 0;
354 }
355
356 /*====================================================================
357 * structured grid regions
358 *====================================================================*/
359
360 #define NODE_INDEX(I,J,K) ((I)+dim[0]*((J)+dim[1]*(K)))
361
362 /*-------------------------------------------------------------------*/
363
structured_range(Regn * reg,cgsize_t * dim,cgsize_t * ptrng,CGNS_ENUMT (GridLocation_t)location)364 static int structured_range (Regn *reg, cgsize_t *dim, cgsize_t *ptrng,
365 CGNS_ENUMT(GridLocation_t) location)
366 {
367 int n, i, j, nf;
368 cgsize_t ii, jj, kk, nfaces, rng[3][2];
369 Face *f;
370 static char *funcname = "structured_range";
371
372 if (location != CGNS_ENUMV(Vertex) &&
373 location != CGNS_ENUMV(IFaceCenter) &&
374 location != CGNS_ENUMV(JFaceCenter) &&
375 location != CGNS_ENUMV(KFaceCenter)) {
376 i = j = 0;
377 for (n = 0; n < CellDim; n++) {
378 if (ptrng[n] == ptrng[n+CellDim] &&
379 (ptrng[n] == 1 || ptrng[n] == dim[n])) {
380 if (ptrng[n] == 1)
381 i++;
382 else if (j) {
383 j = 4;
384 break;
385 }
386 else
387 j = n + 1;
388 }
389 }
390 if (!j && i == 1) {
391 for (n = 0; n < CellDim; n++) {
392 if (ptrng[n] == ptrng[n+CellDim] && ptrng[n] == 1) {
393 j = n + 1;
394 break;
395 }
396 }
397 }
398 if (j == 1)
399 location = CGNS_ENUMV(IFaceCenter);
400 else if (j == 2)
401 location = CGNS_ENUMV(JFaceCenter);
402 else if (j == 3)
403 location = CGNS_ENUMV(KFaceCenter);
404 else {
405 strcpy (reg->errmsg,
406 "unable to determine boundary - use [IJK]FaceCenter");
407 return 0;
408 }
409 }
410
411 nfaces = 1;
412 if (location == CGNS_ENUMV(Vertex)) {
413 for (n = 0, i = 0; i < CellDim; i++) {
414 if (ptrng[i] < 1 || ptrng[i] > dim[i]) return 0;
415 if (ptrng[i] == ptrng[i+CellDim]) {
416 if (n || (ptrng[i] != 1 && ptrng[i] != dim[i]))
417 return 0;
418 n = i + 1;
419 rng[i][0] = rng[i][1] = ptrng[i] - 1;
420 } else {
421 if (ptrng[i] < ptrng[i+CellDim]) {
422 rng[i][0] = ptrng[i] - 1;
423 rng[i][1] = ptrng[i+CellDim] - 1;
424 }
425 else {
426 rng[i][0] = ptrng[i+CellDim] - 1;
427 rng[i][1] = ptrng[i] - 1;
428 }
429 nfaces *= (rng[i][1] - rng[i][0]);
430 }
431 }
432 }
433 else {
434 if (location == CGNS_ENUMV(IFaceCenter))
435 n = 0;
436 else if (location == CGNS_ENUMV(JFaceCenter))
437 n = 1;
438 else
439 n = 2;
440 for (i = 0; i < CellDim; i++) {
441 if (i == n) {
442 if (ptrng[i] != ptrng[i+CellDim] ||
443 (ptrng[i] != 1 && ptrng[i] != dim[i])) return 0;
444 rng[i][0] = rng[i][1] = ptrng[i] - 1;
445 }
446 else {
447 if (ptrng[i] < ptrng[i+CellDim]) {
448 rng[i][0] = ptrng[i] - 1;
449 rng[i][1] = ptrng[i+CellDim];
450 }
451 else {
452 rng[i][0] = ptrng[i+CellDim] - 1;
453 rng[i][1] = ptrng[i];
454 }
455 if (rng[i][0] < 0 || rng[i][1] >= dim[i]) return 0;
456 nfaces *= (rng[i][1] - rng[i][0]);
457 }
458 }
459 n++;
460 }
461 if (!nfaces || n < 1 || n > CellDim) {
462 strcpy (reg->errmsg, "couldn't find any exterior faces");
463 return 0;
464 }
465
466 if (CellDim == 2) {
467 reg->nedges = nfaces;
468 reg->edges = (Edge *) MALLOC (funcname, (size_t)nfaces * sizeof(Edge));
469 nf = 0;
470 kk = 0;
471
472 if (n == 1) {
473 ii = rng[0][0];
474 for (jj = rng[1][0]; jj < rng[1][1]; jj++) {
475 reg->edges[nf].nodes[0] = NODE_INDEX(ii, jj, kk);
476 reg->edges[nf].nodes[1] = NODE_INDEX(ii, jj+1, kk);
477 nf++;
478 }
479 }
480 else {
481 jj = rng[1][0];
482 for (ii = rng[0][0]; ii < rng[0][1]; ii++) {
483 reg->edges[nf].nodes[0] = NODE_INDEX(ii, jj, kk);
484 reg->edges[nf].nodes[1] = NODE_INDEX(ii+1, jj, kk);
485 nf++;
486 }
487 }
488 return 1;
489 }
490
491 reg->nfaces = nfaces;
492 reg->faces = (Face **) MALLOC (funcname, (size_t)nfaces * sizeof(Face *));
493 for (nf = 0; nf < nfaces; nf++)
494 reg->faces[nf] = new_face (funcname, 4);
495 nf = 0;
496
497 if (n == 1) {
498 if ((ii = rng[0][0]) == 0) {
499 for (kk = rng[2][0]; kk < rng[2][1]; kk++) {
500 for (jj = rng[1][0]; jj < rng[1][1]; jj++) {
501 f = reg->faces[nf++];
502 f->nodes[0] = NODE_INDEX (ii, jj, kk);
503 f->nodes[1] = NODE_INDEX (ii, jj, kk+1);
504 f->nodes[2] = NODE_INDEX (ii, jj+1, kk+1);
505 f->nodes[3] = NODE_INDEX (ii, jj+1, kk);
506 }
507 }
508 }
509 else {
510 for (kk = rng[2][0]; kk < rng[2][1]; kk++) {
511 for (jj = rng[1][0]; jj < rng[1][1]; jj++) {
512 f = reg->faces[nf++];
513 f->nodes[0] = NODE_INDEX (ii, jj, kk);
514 f->nodes[1] = NODE_INDEX (ii, jj+1, kk);
515 f->nodes[2] = NODE_INDEX (ii, jj+1, kk+1);
516 f->nodes[3] = NODE_INDEX (ii, jj, kk+1);
517 }
518 }
519 }
520 }
521 else if (n == 2) {
522 if ((jj = rng[1][0]) == 0) {
523 for (kk = rng[2][0]; kk < rng[2][1]; kk++) {
524 for (ii = rng[0][0]; ii < rng[0][1]; ii++) {
525 f = reg->faces[nf++];
526 f->nodes[0] = NODE_INDEX (ii, jj, kk);
527 f->nodes[1] = NODE_INDEX (ii+1, jj, kk);
528 f->nodes[2] = NODE_INDEX (ii+1, jj, kk+1);
529 f->nodes[3] = NODE_INDEX (ii, jj, kk+1);
530 }
531 }
532 }
533 else {
534 for (kk = rng[2][0]; kk < rng[2][1]; kk++) {
535 for (ii = rng[0][0]; ii < rng[0][1]; ii++) {
536 f = reg->faces[nf++];
537 f->nodes[0] = NODE_INDEX (ii, jj, kk);
538 f->nodes[1] = NODE_INDEX (ii, jj, kk+1);
539 f->nodes[2] = NODE_INDEX (ii+1, jj, kk+1);
540 f->nodes[3] = NODE_INDEX (ii+1, jj, kk);
541 }
542 }
543 }
544 }
545 else {
546 if ((kk = rng[2][0]) == 0) {
547 for (jj = rng[1][0]; jj < rng[1][1]; jj++) {
548 for (ii = rng[0][0]; ii < rng[0][1]; ii++) {
549 f = reg->faces[nf++];
550 f->nodes[0] = NODE_INDEX (ii, jj, kk);
551 f->nodes[1] = NODE_INDEX (ii, jj+1, kk);
552 f->nodes[2] = NODE_INDEX (ii+1, jj+1, kk);
553 f->nodes[3] = NODE_INDEX (ii+1, jj, kk);
554 }
555 }
556 }
557 else {
558 for (jj = rng[1][0]; jj < rng[1][1]; jj++) {
559 for (ii = rng[0][0]; ii < rng[0][1]; ii++) {
560 f = reg->faces[nf++];
561 f->nodes[0] = NODE_INDEX (ii, jj, kk);
562 f->nodes[1] = NODE_INDEX (ii+1, jj, kk);
563 f->nodes[2] = NODE_INDEX (ii+1, jj+1, kk);
564 f->nodes[3] = NODE_INDEX (ii, jj+1, kk);
565 }
566 }
567 }
568 }
569 return 2;
570 }
571
572 /*-------------------------------------------------------------------*/
573
structured_list(Regn * mesh,Regn * reg,cgsize_t * dim,cgsize_t npts,cgsize_t * pts,CGNS_ENUMT (GridLocation_t)location)574 static int structured_list (Regn *mesh, Regn *reg, cgsize_t *dim, cgsize_t npts,
575 cgsize_t *pts, CGNS_ENUMT(GridLocation_t) location)
576 {
577 cgsize_t n, nn, nf, nfaces = 0, noff, nmax;
578 Face **faces, *f;
579 static char *funcname = "structured_list";
580
581 if (location != CGNS_ENUMV(Vertex) &&
582 location != CGNS_ENUMV(IFaceCenter) &&
583 location != CGNS_ENUMV(JFaceCenter) &&
584 location != CGNS_ENUMV(KFaceCenter)) {
585 int i, j, k;
586 cgsize_t rng[3][2];
587 for (i = 0; i < CellDim; i++)
588 rng[i][0] = rng[i][1] = pts[i];
589 for (nf = 1; nf < npts; nf++) {
590 nn = nf * CellDim;
591 for (i = 0; i < CellDim; i++) {
592 if (rng[i][0] > pts[nn+i]) rng[i][0] = pts[nn+i];
593 if (rng[i][1] < pts[nn+i]) rng[i][1] = pts[nn+i];
594 }
595 }
596 j = k = 0;
597 for (i = 0; i < CellDim; i++) {
598 if (rng[i][0] == rng[i][1] &&
599 (rng[i][0] == 1 || rng[i][0] == dim[i])) {
600 if (rng[i][0] == 1)
601 j++;
602 else if (k) {
603 k = 4;
604 break;
605 }
606 else
607 k = i + 1;
608 }
609 }
610 if (!k && j == 1) {
611 for (i = 0; i < CellDim; i++) {
612 if (rng[i][0] == rng[i][1] && rng[i][0] == 1) {
613 k = i + 1;
614 break;
615 }
616 }
617 }
618 if (k == 1)
619 location = CGNS_ENUMV(IFaceCenter);
620 else if (k == 2)
621 location = CGNS_ENUMV(JFaceCenter);
622 else if (k == 3)
623 location = CGNS_ENUMV(KFaceCenter);
624 else {
625 strcpy (reg->errmsg,
626 "unable to determine boundary - use [IJK]FaceCenter");
627 return 0;
628 }
629 }
630 nmax = npts;
631
632 if (CellDim == 2) {
633 cgsize_t ii, jj, n0, n1, ne;
634 Edge *edges = (Edge *) MALLOC (funcname, (size_t)nmax * sizeof(Edge));
635
636 ne = 0;
637 if (location == CGNS_ENUMV(Vertex)) {
638 for (nn = 0, n = 0; n < npts; n++) {
639 pts[n] = NODE_INDEX (pts[nn]-1, pts[nn+1]-1, 0);
640 nn += 2;
641 }
642 for (n = 1; n < npts; n++) {
643 if (pts[n] < pts[n-1]) {
644 qsort (pts, (size_t)npts, sizeof(cgsize_t), compare_ints);
645 break;
646 }
647 }
648
649 ne = 0;
650 jj = 0;
651 for (ii = 1; ii < dim[0]; ii++) {
652 n0 = NODE_INDEX(ii-1, jj, 0);
653 n1 = NODE_INDEX(ii, jj, 0);
654 if (find_int(n0, npts, pts) &&
655 find_int(n1, npts, pts)) {
656 if (ne == nmax) {
657 nmax += 100;
658 edges = (Edge *) REALLOC (funcname,
659 (size_t)nmax * sizeof(Edge), edges);
660 }
661 edges[ne].nodes[0] = n0;
662 edges[ne].nodes[1] = n1;
663 ne++;
664 }
665 }
666
667 jj = dim[1] - 1;
668 for (ii = 1; ii < dim[0]; ii++) {
669 n0 = NODE_INDEX(ii-1, jj, 0);
670 n1 = NODE_INDEX(ii, jj, 0);
671 if (find_int(n0, npts, pts) &&
672 find_int(n1, npts, pts)) {
673 if (ne == nmax) {
674 nmax += 100;
675 edges = (Edge *) REALLOC (funcname,
676 (size_t)nmax * sizeof(Edge), edges);
677 }
678 edges[ne].nodes[0] = n0;
679 edges[ne].nodes[1] = n1;
680 ne++;
681 }
682 }
683
684 ii = 0;
685 for (jj = 1; jj < dim[1]; jj++) {
686 n0 = NODE_INDEX(ii, jj-1, 0);
687 n1 = NODE_INDEX(ii, jj, 0);
688 if (find_int(n0, npts, pts) &&
689 find_int(n1, npts, pts)) {
690 if (ne == nmax) {
691 nmax += 100;
692 edges = (Edge *) REALLOC (funcname,
693 (size_t)nmax * sizeof(Edge), edges);
694 }
695 edges[ne].nodes[0] = n0;
696 edges[ne].nodes[1] = n1;
697 ne++;
698 }
699 }
700
701 ii = dim[0] - 1;
702 for (jj = 1; jj < dim[1]; jj++) {
703 n0 = NODE_INDEX(ii, jj-1, 0);
704 n1 = NODE_INDEX(ii, jj, 0);
705 if (find_int(n0, npts, pts) &&
706 find_int(n1, npts, pts)) {
707 if (ne == nmax) {
708 nmax += 100;
709 edges = (Edge *) REALLOC (funcname,
710 (size_t)nmax * sizeof(Edge), edges);
711 }
712 edges[ne].nodes[0] = n0;
713 edges[ne].nodes[1] = n1;
714 ne++;
715 }
716 }
717 }
718
719 else if (location == CGNS_ENUMV(IFaceCenter)) {
720 for (nn = 0, n = 0; n < npts; n++) {
721 if ((pts[nn] == 1 || pts[nn] == dim[0]) &&
722 pts[nn+1] > 0 && pts[nn+1] < dim[1]) {
723 ii = pts[nn] - 1;
724 jj = pts[nn+1] - 1;
725 n0 = NODE_INDEX(ii, jj, 0);
726 n1 = NODE_INDEX(ii, jj+1, 0);
727 edges[ne].nodes[0] = n0;
728 edges[ne].nodes[1] = n1;
729 ne++;
730 }
731 nn += 2;
732 }
733 }
734
735 else if (location == CGNS_ENUMV(JFaceCenter)) {
736 for (nn = 0, n = 0; n < npts; n++) {
737 if ((pts[nn+1] == 1 || pts[nn+1] == dim[1]) &&
738 pts[nn] > 0 && pts[nn] < dim[0]) {
739 ii = pts[nn] - 1;
740 jj = pts[nn+1] - 1;
741 n0 = NODE_INDEX(ii, jj, 0);
742 n1 = NODE_INDEX(ii+1, jj, 0);
743 edges[ne].nodes[0] = n0;
744 edges[ne].nodes[1] = n1;
745 ne++;
746 }
747 nn += 2;
748 }
749 }
750
751 if (ne == 0) {
752 free(edges);
753 strcpy (reg->errmsg, "couldn't find any exterior edges");
754 return 0;
755 }
756
757 reg->nedges = ne;
758 reg->edges = edges;
759 return 1;
760 }
761
762 faces = (Face **) MALLOC (funcname, (size_t)nmax * sizeof(Face *));
763
764 if (location == CGNS_ENUMV(Vertex)) {
765 for (nn = 0, n = 0; n < npts; n++) {
766 pts[n] = NODE_INDEX (pts[nn]-1, pts[nn+1]-1, pts[nn+2]-1);
767 nn += 3;
768 }
769 for (n = 1; n < npts; n++) {
770 if (pts[n] < pts[n-1]) {
771 qsort (pts, (size_t)npts, sizeof(cgsize_t), compare_ints);
772 break;
773 }
774 }
775
776 for (nf = 0; nf < mesh->nfaces; nf++) {
777 f = mesh->faces[nf];
778 for (nn = 0; nn < f->nnodes; nn++) {
779 if (!find_int (f->nodes[nn], npts, pts))
780 break;
781 }
782 if (nn == f->nnodes) {
783 if (nfaces == nmax) {
784 nmax += 100;
785 faces = (Face **) REALLOC (funcname,
786 (size_t)nmax * sizeof(Face *), faces);
787 }
788 faces[nfaces++] = copy_face (funcname, f);
789 }
790 }
791 }
792
793 else if (location == CGNS_ENUMV(IFaceCenter)) {
794 for (n = 0; n < npts; n++) {
795 nn = 3 * n;
796 if ((pts[nn] == 1 || pts[nn] == dim[0]) &&
797 pts[nn+1] > 0 && pts[nn+1] < dim[1] &&
798 pts[nn+2] > 0 && pts[nn+2] < dim[2]) {
799 nf = pts[nn+1]-1 + (pts[nn+2]-1) * (dim[1]-1);
800 if (pts[nn] == dim[0])
801 nf += (dim[1]-1) * (dim[2]-1);
802 faces[nfaces++] = copy_face (funcname, mesh->faces[nf]);
803 }
804 }
805 }
806
807 else if (location == CGNS_ENUMV(JFaceCenter)) {
808 noff = 2 * (dim[1]-1) * (dim[2]-1);
809 for (n = 0; n < npts; n++) {
810 nn = 3 * n;
811 if ((pts[nn+1] == 1 || pts[nn+1] == dim[1]) &&
812 pts[nn] > 0 && pts[nn] < dim[0] &&
813 pts[nn+2] > 0 && pts[nn+2] < dim[2]) {
814 nf = noff + pts[nn]-1 + (pts[nn+2]-1) * (dim[0]-1);
815 if (pts[nn+1] == dim[1])
816 nf += (dim[0]-1) * (dim[2]-1);
817 faces[nfaces++] = copy_face (funcname, mesh->faces[nf]);
818 }
819 }
820 }
821
822 else {
823 noff = 2 * ((dim[1]-1) * (dim[2]-1) + (dim[0]-1) * (dim[2]-1));
824 for (n = 0; n < npts; n++) {
825 nn = 3 * n;
826 if ((pts[nn+2] == 1 || pts[nn+2] == dim[2]) &&
827 pts[nn] > 0 && pts[nn] < dim[0] &&
828 pts[nn+1] > 0 && pts[nn+1] < dim[1]) {
829 nf = noff + pts[nn]-1 + (pts[nn+1]-1) * (dim[0]-1);
830 if (pts[nn+2] == dim[2])
831 nf += (dim[0]-1) * (dim[1]-1);
832 faces[nfaces++] = copy_face (funcname, mesh->faces[nf]);
833 }
834 }
835 }
836
837 if (nfaces == 0) {
838 free (faces);
839 strcpy (reg->errmsg, "couldn't find any exterior faces");
840 return 0;
841 }
842
843 reg->nfaces = nfaces;
844 reg->faces = faces;
845 return 2;
846 }
847
848 /*-------------------------------------------------------------------*/
849
structured_zone(Tcl_Interp * interp,cgsize_t * dim)850 static int structured_zone (Tcl_Interp *interp, cgsize_t *dim)
851 {
852 char name[33], d_name[33];
853 int nints, nconns, nbocos, nholes, ii, nn, nr, nsets;
854 cgsize_t i, j, k, n, ni, nj, nk, nf, ne, fn;
855 cgsize_t npts, *pts, ndpts;
856 cgsize_t range[6], d_range[6];
857 int transform[3];
858 CGNS_ENUMT(GridLocation_t) location;
859 CGNS_ENUMT(GridConnectivityType_t) type;
860 CGNS_ENUMT(PointSetType_t) ptype, d_ptype;
861 CGNS_ENUMT(ZoneType_t) d_ztype;
862 CGNS_ENUMT(DataType_t) datatype;
863 CGNS_ENUMT(BCType_t) bctype;
864 Face *f;
865 Zone *z = &zones[cgnszone-1];
866 static char *funcname = "structured_zone";
867
868 zone_message ("finding exterior faces", NULL);
869 if (cg_n1to1 (cgnsfn, cgnsbase, cgnszone, &nints) ||
870 cg_nconns (cgnsfn, cgnsbase, cgnszone, &nconns) ||
871 cg_nholes (cgnsfn, cgnsbase, cgnszone, &nholes) ||
872 cg_nbocos (cgnsfn, cgnsbase, cgnszone, &nbocos)) {
873 Tcl_SetResult (interp, (char *)cg_get_error(), TCL_STATIC);
874 return 1;
875 }
876 z->nregs = nints + nconns + nholes + nbocos + 1;
877 #ifndef NO_MESH_BOUNDARIES
878 z->nregs += (2 * CellDim);
879 #endif
880 z->regs = (Regn *) MALLOC (funcname, z->nregs * sizeof(Regn));
881 ni = dim[0] - 1;
882 nj = dim[1] - 1;
883 nk = dim[2] - 1;
884 nr = 1;
885
886 /* mesh boundaries */
887
888 strcpy (z->regs[0].name, "<mesh>");
889 z->regs[0].type = REG_MESH;
890 for (n = 0; n < CellDim; n++) {
891 z->regs[0].data[n] = dim[n];
892 range[n] = 1;
893 range[n+CellDim] = dim[n];
894 }
895
896 if (CellDim == 2) {
897 z->regs[0].dim = 2;
898 z->regs[0].nfaces = ni * nj;
899 z->regs[0].faces = (Face **) MALLOC (funcname,
900 (size_t)z->regs[0].nfaces * sizeof(Face *));
901 fn = 0;
902 k = 0;
903 for (j = 0; j < nj; j++) {
904 for (i = 0; i < ni; i++) {
905 f = z->regs[0].faces[fn++] = new_face (funcname, 4);
906 f->nodes[0] = NODE_INDEX (i, j, k);
907 f->nodes[1] = NODE_INDEX (i, j+1, k);
908 f->nodes[2] = NODE_INDEX (i+1, j+1, k);
909 f->nodes[3] = NODE_INDEX (i+1, j, k);
910 }
911 }
912
913 #ifndef NO_CUTTING_PLANE
914 nf = z->regs[0].nfaces;
915 z->regs[0].elemtype = CGNS_ENUMV(QUAD_4);
916 z->regs[0].nelems = nf;
917 z->regs[0].elems = (cgsize_t *) MALLOC (funcname,
918 (size_t)(4 * nf) * sizeof(cgsize_t));
919 for (n = 0, j = 0; j < nf; j++) {
920 f = z->regs[0].faces[j];
921 for (i = 0; i < 4; i++)
922 z->regs[0].elems[n++] = f->nodes[i];
923 }
924 #endif
925
926 #ifndef NO_MESH_BOUNDARIES
927 strcpy (z->regs[nr].name, "<imin>");
928 z->regs[nr].dim = 1;
929 z->regs[nr].type = REG_BNDS;
930 for (nn = 0; nn < 4; nn++)
931 z->regs[nr].data[nn] = range[nn];
932 z->regs[nr].data[2] = 1;
933 z->regs[nr].nedges = nj;
934 z->regs[nr].edges = (Edge *) MALLOC (funcname,
935 (size_t)nj * sizeof(Edge));
936 for (i = 0, j = 0; j < nj; j++) {
937 z->regs[nr].edges[j].nodes[0] = NODE_INDEX(i, j, k);
938 z->regs[nr].edges[j].nodes[1] = NODE_INDEX(i, j+1, k);
939 }
940 nr++;
941
942 strcpy (z->regs[nr].name, "<imax>");
943 z->regs[nr].dim = 1;
944 z->regs[nr].type = REG_BNDS;
945 for (nn = 0; nn < 4; nn++)
946 z->regs[nr].data[nn] = range[nn];
947 z->regs[nr].data[0] = dim[0];
948 z->regs[nr].nedges = nj;
949 z->regs[nr].edges = (Edge *) MALLOC (funcname,
950 (size_t)nj * sizeof(Edge));
951 for (i = ni, j = 0; j < nj; j++) {
952 z->regs[nr].edges[j].nodes[0] = NODE_INDEX(i, j, k);
953 z->regs[nr].edges[j].nodes[1] = NODE_INDEX(i, j+1, k);
954 }
955 nr++;
956
957 strcpy (z->regs[nr].name, "<jmin>");
958 z->regs[nr].dim = 1;
959 z->regs[nr].type = REG_BNDS;
960 for (nn = 0; nn < 4; nn++)
961 z->regs[nr].data[nn] = range[nn];
962 z->regs[nr].data[3] = 1;
963 z->regs[nr].nedges = ni;
964 z->regs[nr].edges = (Edge *) MALLOC (funcname,
965 (size_t)ni * sizeof(Edge));
966 for (j = 0, i = 0; i < ni; i++) {
967 z->regs[nr].edges[i].nodes[0] = NODE_INDEX(i, j, k);
968 z->regs[nr].edges[i].nodes[1] = NODE_INDEX(i+1, j, k);
969 }
970 nr++;
971
972 strcpy (z->regs[nr].name, "<jmax>");
973 z->regs[nr].dim = 1;
974 z->regs[nr].type = REG_BNDS;
975 for (nn = 0; nn < 4; nn++)
976 z->regs[nr].data[nn] = range[nn];
977 z->regs[nr].data[1] = dim[1];
978 z->regs[nr].nedges = ni;
979 z->regs[nr].edges = (Edge *) MALLOC (funcname,
980 (size_t)ni * sizeof(Edge));
981 for (j = nj, i = 0; i < ni; i++) {
982 z->regs[nr].edges[i].nodes[0] = NODE_INDEX(i, j, k);
983 z->regs[nr].edges[i].nodes[1] = NODE_INDEX(i+1, j, k);
984 }
985 nr++;
986 #endif
987 }
988
989 else {
990 z->regs[0].dim = 3;
991
992 #ifndef NO_CUTTING_PLANE
993 z->regs[0].elemtype = CGNS_ENUMV(HEXA_8);
994 z->regs[0].nelems = ni * nj * nk;
995 z->regs[0].elems = pts = (cgsize_t *) MALLOC (funcname,
996 (size_t)(8 * z->regs[0].nelems) * sizeof(cgsize_t));
997
998 for (n = 0, k = 0; k < nk; k++) {
999 for (j = 0; j < nj; j++) {
1000 for (i = 0; i < ni; i++) {
1001 pts[n++] = NODE_INDEX (i, j, k);
1002 pts[n++] = NODE_INDEX (i+1, j, k);
1003 pts[n++] = NODE_INDEX (i+1, j+1, k);
1004 pts[n++] = NODE_INDEX (i, j+1, k);
1005 pts[n++] = NODE_INDEX (i, j, k+1);
1006 pts[n++] = NODE_INDEX (i+1, j, k+1);
1007 pts[n++] = NODE_INDEX (i+1, j+1, k+1);
1008 pts[n++] = NODE_INDEX (i, j+1, k+1);
1009 }
1010 }
1011 }
1012 #endif
1013
1014 z->regs[0].nfaces = 2 * (nj * nk + ni * nk + ni * nj);
1015 z->regs[0].faces = (Face **) MALLOC (funcname,
1016 (size_t)z->regs[0].nfaces * sizeof(Face *));
1017 fn = 0;
1018
1019 for (i = 0, k = 0; k < nk; k++) {
1020 for (j = 0; j < nj; j++) {
1021 f = z->regs[0].faces[fn++] = new_face (funcname, 4);
1022 f->nodes[0] = NODE_INDEX (i, j, k);
1023 f->nodes[1] = NODE_INDEX (i, j, k+1);
1024 f->nodes[2] = NODE_INDEX (i, j+1, k+1);
1025 f->nodes[3] = NODE_INDEX (i, j+1, k);
1026 }
1027 }
1028 for (i = ni, k = 0; k < nk; k++) {
1029 for (j = 0; j < nj; j++) {
1030 f = z->regs[0].faces[fn++] = new_face (funcname, 4);
1031 f->nodes[0] = NODE_INDEX (i, j, k);
1032 f->nodes[1] = NODE_INDEX (i, j+1, k);
1033 f->nodes[2] = NODE_INDEX (i, j+1, k+1);
1034 f->nodes[3] = NODE_INDEX (i, j, k+1);
1035 }
1036 }
1037 for (j = 0, k = 0; k < nk; k++) {
1038 for (i = 0; i < ni; i++) {
1039 f = z->regs[0].faces[fn++] = new_face (funcname, 4);
1040 f->nodes[0] = NODE_INDEX (i, j, k);
1041 f->nodes[1] = NODE_INDEX (i+1, j, k);
1042 f->nodes[2] = NODE_INDEX (i+1, j, k+1);
1043 f->nodes[3] = NODE_INDEX (i, j, k+1);
1044 }
1045 }
1046 for (j = nj, k = 0; k < nk; k++) {
1047 for (i = 0; i < ni; i++) {
1048 f = z->regs[0].faces[fn++] = new_face (funcname, 4);
1049 f->nodes[0] = NODE_INDEX (i, j, k);
1050 f->nodes[1] = NODE_INDEX (i, j, k+1);
1051 f->nodes[2] = NODE_INDEX (i+1, j, k+1);
1052 f->nodes[3] = NODE_INDEX (i+1, j, k);
1053 }
1054 }
1055 for (k = 0, j = 0; j < nj; j++) {
1056 for (i = 0; i < ni; i++) {
1057 f = z->regs[0].faces[fn++] = new_face (funcname, 4);
1058 f->nodes[0] = NODE_INDEX (i, j, k);
1059 f->nodes[1] = NODE_INDEX (i, j+1, k);
1060 f->nodes[2] = NODE_INDEX (i+1, j+1, k);
1061 f->nodes[3] = NODE_INDEX (i+1, j, k);
1062 }
1063 }
1064 for (k = nk, j = 0; j < nj; j++) {
1065 for (i = 0; i < ni; i++) {
1066 f = z->regs[0].faces[fn++] = new_face (funcname, 4);
1067 f->nodes[0] = NODE_INDEX (i, j, k);
1068 f->nodes[1] = NODE_INDEX (i+1, j, k);
1069 f->nodes[2] = NODE_INDEX (i+1, j+1, k);
1070 f->nodes[3] = NODE_INDEX (i, j+1, k);
1071 }
1072 }
1073
1074 #ifndef NO_MESH_BOUNDARIES
1075 fn = 0;
1076 nf = nj * nk;
1077 strcpy (z->regs[nr].name, "<imin>");
1078 z->regs[nr].dim = 2;
1079 z->regs[nr].type = REG_BNDS;
1080 for (nn = 0; nn < 6; nn++)
1081 z->regs[nr].data[nn] = range[nn];
1082 z->regs[nr].data[3] = 1;
1083 z->regs[nr].nfaces = nf;
1084 z->regs[nr].faces = (Face **) MALLOC (funcname,
1085 (size_t)nf * sizeof(Face *));
1086 for (n = 0; n < nf; n++)
1087 z->regs[nr].faces[n] = copy_face (funcname, z->regs->faces[fn++]);
1088 nr++;
1089
1090 strcpy (z->regs[nr].name, "<imax>");
1091 z->regs[nr].dim = 2;
1092 z->regs[nr].type = REG_BNDS;
1093 for (nn = 0; nn < 6; nn++)
1094 z->regs[nr].data[nn] = range[nn];
1095 z->regs[nr].data[0] = dim[0];
1096 z->regs[nr].nfaces = nf;
1097 z->regs[nr].faces = (Face **) MALLOC (funcname,
1098 (size_t)nf * sizeof(Face *));
1099 for (n = 0; n < nf; n++)
1100 z->regs[nr].faces[n] = copy_face (funcname, z->regs->faces[fn++]);
1101 nr++;
1102
1103 nf = ni * nk;
1104 strcpy (z->regs[nr].name, "<jmin>");
1105 z->regs[nr].dim = 2;
1106 z->regs[nr].type = REG_BNDS;
1107 for (nn = 0; nn < 6; nn++)
1108 z->regs[nr].data[nn] = range[nn];
1109 z->regs[nr].data[4] = 1;
1110 z->regs[nr].nfaces = nf;
1111 z->regs[nr].faces = (Face **) MALLOC (funcname,
1112 (size_t)nf * sizeof(Face *));
1113 for (n = 0; n < nf; n++)
1114 z->regs[nr].faces[n] = copy_face (funcname, z->regs->faces[fn++]);
1115 nr++;
1116
1117 strcpy (z->regs[nr].name, "<jmax>");
1118 z->regs[nr].dim = 2;
1119 z->regs[nr].type = REG_BNDS;
1120 for (nn = 0; nn < 6; nn++)
1121 z->regs[nr].data[nn] = range[nn];
1122 z->regs[nr].data[1] = dim[1];
1123 z->regs[nr].nfaces = nf;
1124 z->regs[nr].faces = (Face **) MALLOC (funcname,
1125 (size_t)nf * sizeof(Face *));
1126 for (n = 0; n < nf; n++)
1127 z->regs[nr].faces[n] = copy_face (funcname, z->regs->faces[fn++]);
1128 nr++;
1129
1130 nf = ni * nj;
1131 strcpy (z->regs[nr].name, "<kmin>");
1132 z->regs[nr].dim = 2;
1133 z->regs[nr].type = REG_BNDS;
1134 for (nn = 0; nn < 6; nn++)
1135 z->regs[nr].data[nn] = range[nn];
1136 z->regs[nr].data[5] = 1;
1137 z->regs[nr].nfaces = nf;
1138 z->regs[nr].faces = (Face **) MALLOC (funcname,
1139 (size_t)nf * sizeof(Face *));
1140 for (n = 0; n < nf; n++)
1141 z->regs[nr].faces[n] = copy_face (funcname, z->regs->faces[fn++]);
1142 nr++;
1143
1144 strcpy (z->regs[nr].name, "<kmax>");
1145 z->regs[nr].dim = 2;
1146 z->regs[nr].type = REG_BNDS;
1147 for (nn = 0; nn < 6; nn++)
1148 z->regs[nr].data[nn] = range[nn];
1149 z->regs[nr].data[2] = dim[2];
1150 z->regs[nr].nfaces = nf;
1151 z->regs[nr].faces = (Face **) MALLOC (funcname,
1152 (size_t)nf * sizeof(Face *));
1153 for (n = 0; n < nf; n++)
1154 z->regs[nr].faces[n] = copy_face (funcname, z->regs->faces[fn++]);
1155 nr++;
1156 #endif
1157 }
1158
1159 /* 1 to 1 interfaces */
1160
1161 for (nn = 1; nn <= nints; nn++) {
1162 if (cg_1to1_read (cgnsfn, cgnsbase, cgnszone, nn,
1163 name, d_name, range, d_range, transform)) {
1164 Tcl_SetResult (interp, (char *)cg_get_error(), TCL_STATIC);
1165 return 1;
1166 }
1167 strcpy (z->regs[nr].name, name);
1168 z->regs[nr].type = REG_1TO1;
1169 z->regs[nr].data[0] = 2 * CellDim;
1170 for (ii = 0; ii < 2*CellDim; ii++)
1171 z->regs[nr].data[ii+1] = range[ii];
1172 strcpy (z->regs[nr].d_name, d_name);
1173 k = structured_range (&z->regs[nr], dim, range, CGNS_ENUMV(Vertex));
1174 z->regs[nr].dim = k;
1175 nr++;
1176 }
1177
1178 /* general connectivities */
1179
1180 for (nn = 1; nn <= nconns; nn++) {
1181 if (cg_conn_info (cgnsfn, cgnsbase, cgnszone, nn, name,
1182 &location, &type, &ptype, &npts, d_name, &d_ztype,
1183 &d_ptype, &datatype, &ndpts)) {
1184 Tcl_SetResult (interp, (char *)cg_get_error(), TCL_STATIC);
1185 return 1;
1186 }
1187 pts = (cgsize_t *) MALLOC (funcname, (size_t)(3 * npts) * sizeof(cgsize_t));
1188 if (cg_conn_read_short (cgnsfn, cgnsbase, cgnszone, nn, pts)) {
1189 free (pts);
1190 Tcl_SetResult (interp, (char *)cg_get_error(), TCL_STATIC);
1191 return 1;
1192 }
1193 strcpy (z->regs[nr].name, name);
1194 z->regs[nr].type = REG_CONN;
1195 z->regs[nr].data[0] = type;
1196 z->regs[nr].data[1] = location;
1197 z->regs[nr].data[2] = ptype;
1198 z->regs[nr].data[3] = npts;
1199 if (ptype == CGNS_ENUMV(PointRange)) {
1200 z->regs[nr].data[3] = 2 * CellDim;
1201 for (ii = 0; ii < 2*CellDim; ii++)
1202 z->regs[nr].data[ii+4] = pts[ii];
1203 }
1204 strcpy (z->regs[nr].d_name, d_name);
1205
1206 if (type == CGNS_ENUMV(Abutting1to1) || type == CGNS_ENUMV(Abutting)) {
1207 if (ptype == CGNS_ENUMV(PointRange))
1208 k = structured_range (&z->regs[nr], dim, pts, location);
1209 else if (ptype == CGNS_ENUMV(PointList))
1210 k = structured_list (z->regs, &z->regs[nr], dim, npts, pts, location);
1211 else {
1212 k = 0;
1213 strcpy (z->regs[nr].errmsg, "invalid point set type");
1214 }
1215 }
1216 else if (type == CGNS_ENUMV(Overset)) {
1217 k = 0;
1218 strcpy (z->regs[nr].errmsg, "Overset connectivity not implemented");
1219 }
1220 else {
1221 k = 0;
1222 strcpy (z->regs[nr].errmsg, "invalid connectivity type");
1223 }
1224 z->regs[nr].dim = k;
1225 free (pts);
1226 nr++;
1227 }
1228
1229 /* holes */
1230
1231 for (nn = 1; nn <= nholes; nn++) {
1232 if (cg_hole_info (cgnsfn, cgnsbase, cgnszone, nn, name,
1233 &location, &ptype, &nsets, &npts)) {
1234 Tcl_SetResult (interp, (char *)cg_get_error(), TCL_STATIC);
1235 return 1;
1236 }
1237 pts = (cgsize_t *) MALLOC (funcname, (size_t)(3 * npts * nsets) * sizeof(cgsize_t));
1238 if (cg_hole_read (cgnsfn, cgnsbase, cgnszone, nn, pts)) {
1239 free (pts);
1240 Tcl_SetResult (interp, (char *)cg_get_error(), TCL_STATIC);
1241 return 1;
1242 }
1243 strcpy (z->regs[nr].name, name);
1244 z->regs[nr].type = REG_HOLE;
1245 z->regs[nr].data[0] = nsets;
1246 z->regs[nr].data[1] = location;
1247 z->regs[nr].data[2] = ptype;
1248 z->regs[nr].data[3] = npts;
1249
1250 if (ptype == CGNS_ENUMV(PointRange)) {
1251 z->regs[nr].data[3] = 2 * CellDim;
1252 for (ii = 0; ii < 2*CellDim; ii++)
1253 z->regs[nr].data[ii+4] = pts[ii];
1254 z->regs[nr].dim = structured_range (&z->regs[nr], dim, pts, location);
1255
1256 if (z->regs[nr].dim && nsets > 1) {
1257 Edge *edges = z->regs[nr].edges;
1258 Face **faces = z->regs[nr].faces;
1259 ne = z->regs[nr].nedges;
1260 nf = z->regs[nr].nfaces;
1261 for (ii = 1; ii < nsets; ii++) {
1262 z->regs[nr].nedges = z->regs[nr].nfaces = 0;
1263 k = structured_range (&z->regs[nr], dim,
1264 &pts[ii*2*CellDim], location);
1265 if (k && z->regs[nr].nedges) {
1266 edges = (Edge *) REALLOC (funcname,
1267 (ne + z->regs[nr].nedges) * sizeof(Edge), edges);
1268 for (i = 0; i < z->regs[nr].nedges; i++) {
1269 edges[ne].nodes[0] = z->regs[nr].edges[i].nodes[0];
1270 edges[ne].nodes[1] = z->regs[nr].edges[i].nodes[1];
1271 ne++;
1272 }
1273 free(z->regs[nr].edges);
1274 }
1275 if (k && z->regs[nr].nfaces) {
1276 faces = (Face **) REALLOC (funcname,
1277 (nf + z->regs[nr].nfaces) * sizeof(Face *), faces);
1278 for (i = 0; i < z->regs[nr].nfaces; i++)
1279 faces[nf++] = z->regs[nr].faces[i];
1280 free(z->regs[nr].faces);
1281 }
1282 }
1283 z->regs[nr].nedges = ne;
1284 z->regs[nr].edges = edges;
1285 z->regs[nr].nfaces = nf;
1286 z->regs[nr].faces = faces;
1287 }
1288 }
1289 else if (ptype == CGNS_ENUMV(PointList)) {
1290 z->regs[nr].dim = structured_list (z->regs, &z->regs[nr],
1291 dim, npts, pts, location);
1292 }
1293 else {
1294 strcpy (z->regs[nr].errmsg, "invalid Point Set Type");
1295 }
1296 free(pts);
1297 nr++;
1298 }
1299
1300 /* boundary conditions */
1301
1302 for (nn = 1; nn <= nbocos; nn++) {
1303 if (cg_boco_info (cgnsfn, cgnsbase, cgnszone, nn, name,
1304 &bctype, &ptype, &npts, transform, &j, &datatype, &ii) ||
1305 cg_boco_gridlocation_read(cgnsfn, cgnsbase, cgnszone, nn,
1306 &location)) {
1307 Tcl_SetResult (interp, (char *)cg_get_error(), TCL_STATIC);
1308 return 1;
1309 }
1310 pts = (cgsize_t *) MALLOC (funcname, (size_t)(3 * npts) * sizeof(cgsize_t));
1311 if (cg_boco_read (cgnsfn, cgnsbase, cgnszone, nn, pts, 0)) {
1312 free (pts);
1313 Tcl_SetResult (interp, (char *)cg_get_error(), TCL_STATIC);
1314 return 1;
1315 }
1316 strcpy (z->regs[nr].name, name);
1317 z->regs[nr].type = REG_BOCO;
1318 z->regs[nr].data[0] = bctype;
1319 z->regs[nr].data[1] = location;
1320 z->regs[nr].data[2] = ptype;
1321 z->regs[nr].data[3] = i;
1322 if (ptype == CGNS_ENUMV(PointRange) || ptype == CGNS_ENUMV(ElementRange)) {
1323 z->regs[nr].data[3] = 2 * CellDim;
1324 for (ii = 0; ii < 2*CellDim; ii++)
1325 z->regs[nr].data[ii+4] = pts[ii];
1326 }
1327
1328 if (ptype == CGNS_ENUMV(PointRange) || ptype == CGNS_ENUMV(ElementRange))
1329 k = structured_range (&z->regs[nr], dim, pts, location);
1330 else if (ptype == CGNS_ENUMV(PointList) || ptype == CGNS_ENUMV(ElementList))
1331 k = structured_list (z->regs, &z->regs[nr], dim, npts, pts, location);
1332 else {
1333 k = 0;
1334 strcpy (z->regs[nr].errmsg, "invalid point set type");
1335 }
1336 z->regs[nr].dim = k;
1337 free (pts);
1338 nr++;
1339 }
1340
1341 z->nregs = nr;
1342
1343 #ifndef NO_CUTTING_PLANE
1344 for (nr = 0; nr < z->nregs; nr++) {
1345 if (z->regs[nr].dim == 2 && z->regs[nr].nfaces) {
1346 nf = z->regs[nr].nfaces;
1347 z->regs[nr].elemtype = CGNS_ENUMV(QUAD_4);
1348 z->regs[nr].nelems = nf;
1349 z->regs[nr].elems = (cgsize_t *) MALLOC (funcname,
1350 (size_t)(4 * nf) * sizeof(cgsize_t));
1351 for (n = 0, j = 0; j < nf; j++, f++) {
1352 f = z->regs[nr].faces[j];
1353 for (ii = 0; ii < 4; ii++)
1354 z->regs[nr].elems[n++] = f->nodes[ii];
1355 }
1356 }
1357 }
1358 #endif
1359
1360 return 0;
1361 }
1362
1363 /*====================================================================
1364 * unstructured grid regions
1365 *====================================================================*/
1366
1367 static int max_face_nodes = 0;
1368 static cgsize_t *sort_face_nodes = NULL;
1369
1370 /*-------------------------------------------------------------------*/
1371
compare_faces(void * v1,void * v2)1372 static int compare_faces (void *v1, void *v2)
1373 {
1374 Face *f1 = (Face *)v1;
1375 Face *f2 = (Face *)v2;
1376 int i, k;
1377 cgsize_t id, nn, *n1, *n2;
1378
1379 if (f1->nnodes != f2->nnodes)
1380 return (f1->nnodes - f2->nnodes);
1381
1382 if (f1->nnodes > max_face_nodes) {
1383 max_face_nodes += 10;
1384 sort_face_nodes = (cgsize_t *) REALLOC ("compare_faces",
1385 2 * max_face_nodes * sizeof(cgsize_t), sort_face_nodes);
1386 }
1387 n1 = sort_face_nodes;
1388 n2 = sort_face_nodes + max_face_nodes;
1389
1390 for (i = 0; i < f1->nnodes; i++) {
1391 id = f1->nodes[i];
1392 for (k = 0; k < i; k++) {
1393 if (n1[k] > id) {
1394 nn = n1[k];
1395 n1[k] = id;
1396 id = nn;
1397 }
1398 }
1399 n1[i] = id;
1400 }
1401 for (i = 0; i < f2->nnodes; i++) {
1402 id = f2->nodes[i];
1403 for (k = 0; k < i; k++) {
1404 if (n2[k] > id) {
1405 nn = n2[k];
1406 n2[k] = id;
1407 id = nn;
1408 }
1409 }
1410 n2[i] = id;
1411 }
1412
1413 for (i = 0; i < f1->nnodes; i++) {
1414 if (n1[i] != n2[i])
1415 return (int)(n1[i] - n2[i]);
1416 }
1417 return 0;
1418 }
1419
1420 /*-------------------------------------------------------------------*/
1421
hash_face(void * v)1422 static size_t hash_face (void *v)
1423 {
1424 Face *f = (Face *)v;
1425 int n;
1426 size_t hash = 0;
1427
1428 for (n = 0; n < f->nnodes; n++)
1429 hash += (size_t)f->nodes[n];
1430 return hash;
1431 }
1432
1433 /*-------------------------------------------------------------------*/
1434
get_faces(void * vf,void * vr)1435 static size_t get_faces (void *vf, void *vr)
1436 {
1437 Face *f = (Face *)vf;
1438 Regn *r = (Regn *)vr;
1439
1440 r->faces[r->nfaces] = f;
1441 (r->nfaces)++;
1442 return 1;
1443 }
1444
1445 /*-------------------------------------------------------------------*/
1446
compare_poly(void * v1,void * v2)1447 static int compare_poly (void *v1, void *v2)
1448 {
1449 PolyFace *p1 = (PolyFace *)v1;
1450 PolyFace *p2 = (PolyFace *)v2;
1451
1452 return compare_faces(p1->face, p2->face);
1453 }
1454
1455 /*-------------------------------------------------------------------*/
1456
hash_poly(void * v)1457 static size_t hash_poly (void *v)
1458 {
1459 PolyFace *p = (PolyFace *)v;
1460
1461 return hash_face(p->face);
1462 }
1463
1464 /*-------------------------------------------------------------------*/
1465
1466 static cgsize_t nPolyFaces;
1467
poly_faces(void * vp,void * vl)1468 static size_t poly_faces (void *vp, void *vl)
1469 {
1470 PolyFace *pf = (PolyFace *)vp;
1471 PolyFace **pl = (PolyFace **)vl;
1472
1473 pl[nPolyFaces++] = pf;
1474 return 1;
1475 }
1476
1477 /*-------------------------------------------------------------------*/
1478
poly_sort(const void * v1,const void * v2)1479 static int poly_sort (const void *v1, const void *v2)
1480 {
1481 const PolyFace **p1 = (const PolyFace **)v1;
1482 const PolyFace **p2 = (const PolyFace **)v2;
1483
1484 return (int)((*p1)->num - (*p2)->num);
1485 }
1486
1487 /*-------------------------------------------------------------------*/
1488
find_face(Zone * z,cgsize_t fnum)1489 static Face *find_face (Zone *z, cgsize_t fnum)
1490 {
1491 int nr;
1492 cgsize_t nf;
1493
1494 for (nr = 0; nr < z->nregs; nr++) {
1495 if (z->regs[nr].type == REG_ELEM &&
1496 z->regs[nr].dim == 2 &&
1497 z->regs[nr].data[1] <= fnum &&
1498 z->regs[nr].data[2] >= fnum) {
1499 nf = fnum - z->regs[nr].data[1];
1500 return z->regs[nr].faces[nf];
1501 }
1502 }
1503 return NULL;
1504 }
1505
1506 /*-------------------------------------------------------------------*/
1507
element_dimension(CGNS_ENUMT (ElementType_t)elemtype)1508 static int element_dimension (CGNS_ENUMT(ElementType_t) elemtype)
1509 {
1510 switch (elemtype) {
1511 case CGNS_ENUMV(NODE):
1512 return 0;
1513 case CGNS_ENUMV(BAR_2):
1514 case CGNS_ENUMV(BAR_3):
1515 case CGNS_ENUMV(BAR_4):
1516 case CGNS_ENUMV(BAR_5):
1517 return 1;
1518 case CGNS_ENUMV(TRI_3):
1519 case CGNS_ENUMV(TRI_6):
1520 case CGNS_ENUMV(TRI_9):
1521 case CGNS_ENUMV(TRI_10):
1522 case CGNS_ENUMV(TRI_12):
1523 case CGNS_ENUMV(TRI_15):
1524 case CGNS_ENUMV(QUAD_4):
1525 case CGNS_ENUMV(QUAD_8):
1526 case CGNS_ENUMV(QUAD_9):
1527 case CGNS_ENUMV(QUAD_12):
1528 case CGNS_ENUMV(QUAD_16):
1529 case CGNS_ENUMV(QUAD_P4_16):
1530 case CGNS_ENUMV(QUAD_25):
1531 case CGNS_ENUMV(NGON_n):
1532 return 2;
1533 case CGNS_ENUMV(TETRA_4):
1534 case CGNS_ENUMV(TETRA_10):
1535 case CGNS_ENUMV(TETRA_16):
1536 case CGNS_ENUMV(TETRA_20):
1537 case CGNS_ENUMV(TETRA_22):
1538 case CGNS_ENUMV(TETRA_34):
1539 case CGNS_ENUMV(TETRA_35):
1540 case CGNS_ENUMV(PYRA_5):
1541 case CGNS_ENUMV(PYRA_13):
1542 case CGNS_ENUMV(PYRA_14):
1543 case CGNS_ENUMV(PYRA_21):
1544 case CGNS_ENUMV(PYRA_29):
1545 case CGNS_ENUMV(PYRA_30):
1546 case CGNS_ENUMV(PYRA_P4_29):
1547 case CGNS_ENUMV(PYRA_50):
1548 case CGNS_ENUMV(PYRA_55):
1549 case CGNS_ENUMV(PENTA_6):
1550 case CGNS_ENUMV(PENTA_15):
1551 case CGNS_ENUMV(PENTA_18):
1552 case CGNS_ENUMV(PENTA_24):
1553 case CGNS_ENUMV(PENTA_38):
1554 case CGNS_ENUMV(PENTA_40):
1555 case CGNS_ENUMV(PENTA_33):
1556 case CGNS_ENUMV(PENTA_66):
1557 case CGNS_ENUMV(PENTA_75):
1558 case CGNS_ENUMV(HEXA_8):
1559 case CGNS_ENUMV(HEXA_20):
1560 case CGNS_ENUMV(HEXA_27):
1561 case CGNS_ENUMV(HEXA_32):
1562 case CGNS_ENUMV(HEXA_56):
1563 case CGNS_ENUMV(HEXA_64):
1564 case CGNS_ENUMV(HEXA_44):
1565 case CGNS_ENUMV(HEXA_98):
1566 case CGNS_ENUMV(HEXA_125):
1567 case CGNS_ENUMV(NFACE_n):
1568 return 3;
1569 default:
1570 break;
1571 }
1572 return -1;
1573 }
1574
1575 /*-------------------------------------------------------------------*/
1576
edge_elements(Regn * r,cgsize_t * conn)1577 static void edge_elements (Regn *r, cgsize_t *conn)
1578 {
1579 int ip;
1580 cgsize_t istart, n, ne, nelems;
1581
1582 istart = r->data[1];
1583 nelems = r->data[2] - istart + 1;
1584 cg_npe ((CGNS_ENUMT(ElementType_t))r->data[0], &ip);
1585
1586 r->nedges = nelems;
1587 r->edges = (Edge *) MALLOC ("edge_elements", nelems * sizeof(Edge));
1588 for (n = 0, ne = 0; ne < nelems; ne++) {
1589 r->edges[ne].id = istart + ne;
1590 r->edges[ne].nodes[0] = conn[n] - 1;
1591 r->edges[ne].nodes[1] = conn[n+1] - 1;
1592 n += ip;
1593 }
1594 }
1595
1596 /*-------------------------------------------------------------------*/
1597
face_elements(Regn * r,cgsize_t * conn,cgsize_t * conn_offsets)1598 static void face_elements (Regn *r, cgsize_t *conn, cgsize_t *conn_offsets)
1599 {
1600 int i, ip;
1601 cgsize_t ne, nn, istart, nelems;
1602 cgsize_t rind0, rind1;
1603 CGNS_ENUMT(ElementType_t) elemtype, type;
1604 static char *funcname = "face_elements";
1605
1606 elemtype = type = (CGNS_ENUMT(ElementType_t))r->data[0];
1607 istart = r->data[1];
1608 nelems = r->data[2] - istart + 1;
1609 rind0 = r->data[3];
1610 rind1 = nelems - r->data[4];
1611
1612 r->nfaces = nelems;
1613 r->faces = (Face **) MALLOC (funcname, nelems * sizeof(Face *));
1614
1615 for (nn = 0, ne = 0; ne < nelems; ne++) {
1616 if (elemtype == CGNS_ENUMV(MIXED))
1617 type = (CGNS_ENUMT(ElementType_t))conn[nn++];
1618 switch (type) {
1619 case CGNS_ENUMV(TRI_3):
1620 case CGNS_ENUMV(TRI_6):
1621 case CGNS_ENUMV(TRI_9):
1622 case CGNS_ENUMV(TRI_10):
1623 case CGNS_ENUMV(TRI_12):
1624 case CGNS_ENUMV(TRI_15):
1625 ip = 3;
1626 break;
1627 case CGNS_ENUMV(QUAD_4):
1628 case CGNS_ENUMV(QUAD_8):
1629 case CGNS_ENUMV(QUAD_9):
1630 case CGNS_ENUMV(QUAD_12):
1631 case CGNS_ENUMV(QUAD_16):
1632 case CGNS_ENUMV(QUAD_P4_16):
1633 case CGNS_ENUMV(QUAD_25):
1634 ip = 4;
1635 break;
1636 case CGNS_ENUMV(NGON_n):
1637 ip = (int)(conn_offsets[ne+1]-conn_offsets[ne]);
1638 break;
1639 default:
1640 if (type < CGNS_ENUMV(NODE) ||
1641 type >= NofValidElementTypes)
1642 FATAL ("face_elements:unknown element type");
1643 ip = 0;
1644 break;
1645 }
1646 if (ip) {
1647 r->faces[ne] = new_face (funcname, ip);
1648 r->faces[ne]->id = istart + ne;
1649 for (i = 0; i < ip; i++)
1650 r->faces[ne]->nodes[i] = conn[nn+i] - 1;
1651 if (ne < rind0 || ne >= rind1)
1652 r->faces[ne]->flags = 1;
1653 }
1654 else {
1655 r->faces[ne] = NULL;
1656 }
1657 if (type == CGNS_ENUMV(NGON_n))
1658 nn += ip;
1659 else {
1660 cg_npe (type, &i);
1661 nn += i;
1662 }
1663 }
1664 }
1665
1666 /*-------------------------------------------------------------------*/
1667
exterior_faces(Zone * z,Regn * r,cgsize_t * conn,cgsize_t * conn_offsets)1668 static void exterior_faces (Zone *z, Regn *r, cgsize_t *conn, cgsize_t *conn_offsets)
1669 {
1670 int i, j, nf, ip, flag;
1671 cgsize_t ne, nn, istart, nelems;
1672 cgsize_t rind0, rind1;
1673 CGNS_ENUMT(ElementType_t) elemtype;
1674 HASH *facehash;
1675 Face *face, *pf;
1676 static char *funcname = "exterior_faces";
1677
1678 elemtype = (CGNS_ENUMT(ElementType_t))r->data[0];
1679 istart = r->data[1];
1680 nelems = r->data[2] - istart + 1;
1681 rind0 = r->data[3];
1682 rind1 = nelems - r->data[4];
1683
1684 facehash = HashCreate (nelems > 1024 ? (size_t)nelems / 3 : 127,
1685 compare_faces, hash_face);
1686 if (NULL == facehash)
1687 FATAL ("exterior_faces:face hash table creation failed");
1688
1689 if (elemtype == CGNS_ENUMV(NFACE_n)) {
1690 for (ne = 0; ne < nelems; ne++) {
1691 flag = (ne < rind0 || ne >= rind1) ? 1 : 0;
1692 for (j = conn_offsets[ne]; j < conn_offsets[ne+1]; j++) {
1693 face = find_face (z, abs(conn[j]));
1694 if (face != NULL) {
1695 pf = (Face *) HashFind (facehash, face);
1696 if (NULL == pf) {
1697 pf = copy_face (funcname, face);
1698 pf->id = 0;
1699 pf->flags = flag;
1700 (void) HashAdd (facehash, pf);
1701 }
1702 else if (flag == pf->flags) {
1703 HashDelete (facehash, pf);
1704 free (pf);
1705 }
1706 else {
1707 pf->flags = 0;
1708 }
1709 }
1710 }
1711 }
1712 }
1713 else {
1714 CGNS_ENUMT(ElementType_t) type = elemtype;
1715 face = new_face(funcname, 4);
1716 for (nn = 0, ne = 0; ne < nelems; ne++) {
1717 flag = (ne < rind0 || ne >= rind1) ? 1 : 0;
1718 if (elemtype == CGNS_ENUMV(MIXED))
1719 type = (CGNS_ENUMT(ElementType_t))conn[nn++];
1720 switch (type) {
1721 case CGNS_ENUMV(TETRA_4):
1722 case CGNS_ENUMV(TETRA_10):
1723 case CGNS_ENUMV(TETRA_16):
1724 case CGNS_ENUMV(TETRA_20):
1725 case CGNS_ENUMV(TETRA_22):
1726 case CGNS_ENUMV(TETRA_34):
1727 case CGNS_ENUMV(TETRA_35):
1728 ip = 2;
1729 nf = 4;
1730 break;
1731 case CGNS_ENUMV(PYRA_5):
1732 case CGNS_ENUMV(PYRA_13):
1733 case CGNS_ENUMV(PYRA_14):
1734 case CGNS_ENUMV(PYRA_21):
1735 case CGNS_ENUMV(PYRA_29):
1736 case CGNS_ENUMV(PYRA_30):
1737 case CGNS_ENUMV(PYRA_P4_29):
1738 case CGNS_ENUMV(PYRA_50):
1739 case CGNS_ENUMV(PYRA_55):
1740 ip = 6;
1741 nf = 5;
1742 break;
1743 case CGNS_ENUMV(PENTA_6):
1744 case CGNS_ENUMV(PENTA_15):
1745 case CGNS_ENUMV(PENTA_18):
1746 case CGNS_ENUMV(PENTA_24):
1747 case CGNS_ENUMV(PENTA_38):
1748 case CGNS_ENUMV(PENTA_40):
1749 case CGNS_ENUMV(PENTA_33):
1750 case CGNS_ENUMV(PENTA_66):
1751 case CGNS_ENUMV(PENTA_75):
1752 ip = 11;
1753 nf = 5;
1754 break;
1755 case CGNS_ENUMV(HEXA_8):
1756 case CGNS_ENUMV(HEXA_20):
1757 case CGNS_ENUMV(HEXA_27):
1758 case CGNS_ENUMV(HEXA_32):
1759 case CGNS_ENUMV(HEXA_56):
1760 case CGNS_ENUMV(HEXA_64):
1761 case CGNS_ENUMV(HEXA_44):
1762 case CGNS_ENUMV(HEXA_98):
1763 case CGNS_ENUMV(HEXA_125):
1764 ip = 16;
1765 nf = 6;
1766 break;
1767 default:
1768 if (type < CGNS_ENUMV(NODE) ||
1769 type >= NofValidElementTypes)
1770 FATAL ("exterior_faces:unknown element type");
1771 ip = 0;
1772 nf = 0;
1773 break;
1774 }
1775 for (j = 0; j < nf; j++) {
1776 face->nnodes = facenodes[ip+j][0];
1777 for (i = 0; i < face->nnodes; i++)
1778 face->nodes[i] = conn[nn + facenodes[ip+j][i+1]] - 1;
1779 pf = (Face *) HashFind (facehash, face);
1780 if (NULL == pf) {
1781 pf = copy_face (funcname, face);
1782 pf->id = 0;
1783 pf->flags = flag;
1784 (void) HashAdd (facehash, pf);
1785 }
1786 else if (flag == pf->flags) {
1787 HashDelete (facehash, pf);
1788 free (pf);
1789 }
1790 else {
1791 pf->flags = 0;
1792 }
1793 }
1794 cg_npe (type, &j);
1795 nn += j;
1796 }
1797 free(face);
1798 }
1799
1800 r->nfaces = 0;
1801 ne = (cgsize_t) HashSize (facehash);
1802 if (ne) {
1803 r->faces = (Face **) MALLOC (funcname, ne * sizeof(Face *));
1804 HashList (facehash, get_faces, r);
1805 }
1806 else {
1807 strcpy (r->errmsg, "couldn't find any exterior faces");
1808 }
1809 HashDestroy (facehash, NULL);
1810 }
1811
1812 /*-------------------------------------------------------------------*/
1813
polyhedral_faces(Zone * z,Regn * r,cgsize_t * conn,cgsize_t * conn_offsets)1814 static void polyhedral_faces (Zone *z, Regn *r, cgsize_t *conn, cgsize_t *conn_offsets)
1815 {
1816 int j, nf, flag;
1817 cgsize_t ne, nn, istart, nelems, nfaces, id;
1818 cgsize_t rind0, rind1;
1819 HASH *facehash;
1820 Face *face;
1821 PolyFace poly, *pf, **polylist;
1822 static char *funcname = "polyhedral_faces";
1823
1824 istart = r->data[1];
1825 nelems = r->data[2] - istart + 1;
1826 rind0 = r->data[3];
1827 rind1 = nelems - r->data[4];
1828
1829 facehash = HashCreate (nelems > 1024 ? (size_t)nelems / 3 : 127,
1830 compare_poly, hash_poly);
1831 if (NULL == facehash)
1832 FATAL ("polyhedral_faces:face hash table creation failed");
1833
1834 nfaces = 0;
1835 for (nn = 0, ne = 0; ne < nelems; ne++) {
1836 flag = (ne < rind0 || ne >= rind1) ? 1 : 0;
1837 for (j = conn_offsets[ne]; j < conn_offsets[ne+1]; j++) {
1838 id = conn[j];
1839 face = find_face (z, abs(id));
1840 poly.face = face;
1841 pf = (PolyFace *) HashFind (facehash, &poly);
1842 if (NULL == pf) {
1843 pf = (PolyFace *)MALLOC(funcname, sizeof(PolyFace));
1844 pf->face = face;
1845 pf->num = ++nfaces;
1846 pf->flags = flag;
1847 (void) HashAdd (facehash, pf);
1848 }
1849 else {
1850 if ((pf->flags & 1) != flag)
1851 pf->flags &= ~1;
1852 pf->flags |= 2;
1853 }
1854 conn[j] = id < 0 ? -(pf->num) : pf->num;
1855 }
1856 }
1857
1858 nfaces = (cgsize_t) HashSize (facehash);
1859 polylist = (PolyFace **) MALLOC (funcname, nfaces * sizeof(PolyFace *));
1860 nPolyFaces = 0;
1861 HashList (facehash, poly_faces, polylist);
1862 HashDestroy (facehash, NULL);
1863
1864 qsort(polylist, nfaces, sizeof(PolyFace *), poly_sort);
1865
1866 for (nn = 0, ne = 0; ne < nfaces; ne++) {
1867 if ((polylist[ne]->flags & 2) == 0) nn++;
1868 }
1869 r->nfaces = nn;
1870 r->faces = (Face **) MALLOC (funcname, nn * sizeof(Face *));
1871 for (nn = 0, ne = 0; ne < nfaces; ne++) {
1872 if ((polylist[ne]->flags & 2) == 0) {
1873 r->faces[nn] = copy_face(funcname, polylist[ne]->face);
1874 r->faces[nn]->id = 0;
1875 r->faces[nn]->flags = (polylist[ne]->flags & 1);
1876 nn++;
1877 }
1878 }
1879
1880 r->npoly = nfaces;
1881 r->poly = (Face **) MALLOC (funcname, nfaces * sizeof(Face *));
1882 for (nn = 0; nn < nfaces; nn++) {
1883 face = polylist[nn]->face;
1884 face->flags |= (polylist[nn]->flags & 2);
1885 r->poly[nn] = face;
1886 free(polylist[nn]);
1887 }
1888 free(polylist);
1889
1890 r->elemtype = CGNS_ENUMV(NFACE_n);
1891 r->nelems = nelems;
1892 r->elems = conn;
1893 r->elem_offsets = conn_offsets;
1894 }
1895
1896 /*-------------------------------------------------------------------*/
1897
sort_elemsets(const void * v1,const void * v2)1898 static int sort_elemsets(const void *v1, const void *v2)
1899 {
1900 return (((Regn *)v2)->dim - ((Regn *)v1)->dim);
1901 }
1902
1903 /*-------------------------------------------------------------------*/
1904
unstructured_region(int nregs,Regn * regs,Regn * r,CGNS_ENUMT (PointSetType_t)ptype,cgsize_t nlist,cgsize_t * list)1905 static cgsize_t unstructured_region (int nregs, Regn *regs, Regn *r,
1906 CGNS_ENUMT(PointSetType_t) ptype,
1907 cgsize_t nlist, cgsize_t *list)
1908 {
1909 int nr, nn;
1910 cgsize_t nf, nfaces, maxfaces;
1911 Face **faces, *f;
1912 static char *funcname = "unstructured_region";
1913
1914 if (ptype == CGNS_ENUMV(PointList) ||
1915 ptype == CGNS_ENUMV(ElementList)) {
1916 for (nf = 1; nf < nlist; nf++) {
1917 if (list[nf] < list[nf-1]) {
1918 qsort (list, (size_t)nlist, sizeof(cgsize_t), compare_ints);
1919 break;
1920 }
1921 }
1922 maxfaces = nlist;
1923 }
1924 else if (ptype == CGNS_ENUMV(PointRange) ||
1925 ptype == CGNS_ENUMV(ElementRange)) {
1926 if (list[0] > list[1]) {
1927 nf = list[0];
1928 list[0] = list[1];
1929 list[1] = nf;
1930 }
1931 maxfaces = list[1] - list[0] + 1;
1932 }
1933 else {
1934 strcpy (r->errmsg, "invalid point set type");
1935 maxfaces = 0;
1936 }
1937
1938 if (maxfaces < 1) return 0;
1939 faces = (Face **)MALLOC(funcname, (size_t)maxfaces * sizeof(Face *));
1940
1941 nfaces = 0;
1942 for (nr = 0; nr < nregs; nr++) {
1943 if (!regs[nr].nfaces) continue;
1944 for (nf = 0; nf < regs[nr].nfaces; nf++) {
1945 f = regs[nr].faces[nf];
1946 switch (ptype) {
1947 case CGNS_ENUMV(PointList):
1948 if (f->id) continue;
1949 for (nn = 0; nn < f->nnodes; nn++) {
1950 if (!find_int (f->nodes[nn]+1, nlist, list))
1951 break;
1952 }
1953 if (nn == f->nnodes) break;
1954 continue;
1955 case CGNS_ENUMV(PointRange):
1956 if (f->id) continue;
1957 for (nn = 0; nn < f->nnodes; nn++) {
1958 if (f->nodes[nn]+1 < list[0] ||
1959 f->nodes[nn]+1 > list[1])
1960 break;
1961 }
1962 if (nn == f->nnodes) break;
1963 continue;
1964 case CGNS_ENUMV(ElementList):
1965 if (f->id && find_int (f->id, nlist, list)) break;
1966 continue;
1967 case CGNS_ENUMV(ElementRange):
1968 if (f->id >= list[0] && f->id <= list[1]) break;
1969 continue;
1970 default:
1971 continue;
1972 }
1973 if (nfaces == maxfaces) {
1974 maxfaces += 100;
1975 faces = (Face **) REALLOC (funcname,
1976 (size_t)maxfaces * sizeof(Face *), faces);
1977 }
1978 faces[nfaces++] = f;
1979 }
1980 }
1981 if (nfaces) {
1982 r->nfaces = nfaces;
1983 r->faces = (Face **)MALLOC(funcname, (size_t)nfaces * sizeof(Face *));
1984 for (nf = 0; nf < nfaces; nf++)
1985 r->faces[nf] = copy_face(funcname, faces[nf]);
1986 }
1987 else
1988 strcpy (r->errmsg, "couldn't find any exterior faces");
1989 free (faces);
1990 return nfaces;
1991 }
1992
1993 /*-------------------------------------------------------------------*/
1994
unstructured_zone(Tcl_Interp * interp)1995 static int unstructured_zone (Tcl_Interp *interp)
1996 {
1997 int i, ns, nb, ip, nr, haspoly, nsets;
1998 int nsect, nints, nconns, nholes, nbocos, nrmlindex[3];
1999 int transform[3], rind[2];
2000 cgsize_t is, ie, np, n, ne, nf;
2001 cgsize_t nelem, elemsize, *conn, *conn_offsets;
2002 cgsize_t range[6], d_range[6];
2003 CGNS_ENUMT(GridLocation_t) location;
2004 CGNS_ENUMT(GridConnectivityType_t) type;
2005 CGNS_ENUMT(PointSetType_t) ptype, d_ptype;
2006 CGNS_ENUMT(ZoneType_t) d_ztype;
2007 CGNS_ENUMT(DataType_t) datatype;
2008 CGNS_ENUMT(BCType_t) bctype;
2009 CGNS_ENUMT(ElementType_t) elemtype;
2010 char name[33], d_name[33];
2011 Zone *z = &zones[cgnszone-1];
2012 static char *funcname = "unstructured_zone";
2013 static char *dspmsg = "finding exterior faces for";
2014
2015 if (cg_nsections (cgnsfn, cgnsbase, cgnszone, &nsect) ||
2016 cg_n1to1 (cgnsfn, cgnsbase, cgnszone, &nints) ||
2017 cg_nconns (cgnsfn, cgnsbase, cgnszone, &nconns) ||
2018 cg_nholes (cgnsfn, cgnsbase, cgnszone, &nholes) ||
2019 cg_nbocos (cgnsfn, cgnsbase, cgnszone, &nbocos)) {
2020 Tcl_SetResult (interp, (char *)cg_get_error(), TCL_STATIC);
2021 return 1;
2022 }
2023 z->nregs = nsect + nints + nconns + nholes + nbocos;
2024 z->regs = (Regn *) MALLOC (funcname, z->nregs * sizeof(Regn));
2025
2026 /* element sets */
2027
2028 haspoly = 0;
2029 for (nr = 0, ns = 1; ns <= nsect; ns++, nr++) {
2030 if (cg_section_read (cgnsfn, cgnsbase, cgnszone, ns,
2031 name, &elemtype, &is, &ie, &nb, &ip) ||
2032 cg_ElementDataSize (cgnsfn, cgnsbase, cgnszone, ns, &elemsize)) {
2033 Tcl_SetResult (interp, (char *)cg_get_error(), TCL_STATIC);
2034 return 1;
2035 }
2036 zone_message (dspmsg, name);
2037 strcpy (z->regs[nr].name, name);
2038 z->regs[nr].type = REG_ELEM;
2039 z->regs[nr].data[0] = elemtype;
2040 z->regs[nr].data[1] = is;
2041 z->regs[nr].data[2] = ie;
2042
2043 if (cg_goto (cgnsfn, cgnsbase, "Zone_t", cgnszone,
2044 "Elements_t", ns, "end") ||
2045 cg_rind_read (rind)) {
2046 rind[0] = rind[1] = 0;
2047 }
2048 z->regs[nr].data[3] = rind[0];
2049 z->regs[nr].data[4] = rind[1];
2050
2051 if (elemtype < CGNS_ENUMV(BAR_2) ||
2052 elemtype >= NofValidElementTypes) {
2053 strcpy (z->regs[nr].errmsg, "invalid element type");
2054 continue;
2055 }
2056
2057 /* do this after reading all the sections */
2058
2059 if (elemtype == CGNS_ENUMV(NFACE_n)) {
2060 z->regs[nr].dim = 3;
2061 haspoly++;
2062 continue;
2063 }
2064
2065 nelem = ie - is + 1;
2066 conn = (cgsize_t *) MALLOC (funcname, (size_t)elemsize * sizeof(cgsize_t));
2067 conn_offsets = NULL;
2068 if (elemtype == CGNS_ENUMV(MIXED) ||
2069 elemtype == CGNS_ENUMV(NGON_n) ) {
2070 conn_offsets = (cgsize_t *) MALLOC (funcname, (size_t) (nelem+1)*sizeof(cgsize_t));
2071 if (cg_poly_elements_read (cgnsfn, cgnsbase, cgnszone, ns, conn, conn_offsets, 0)) {
2072 free (conn);
2073 if (conn_offsets) {
2074 free(conn_offsets);
2075 }
2076 Tcl_SetResult (interp, (char *)cg_get_error(), TCL_STATIC);
2077 return 1;
2078 }
2079 }
2080 else {
2081 if (cg_elements_read (cgnsfn, cgnsbase, cgnszone, ns, conn, 0)) {
2082 free (conn);
2083 Tcl_SetResult (interp, (char *)cg_get_error(), TCL_STATIC);
2084 return 1;
2085 }
2086 }
2087
2088 /* check element indices */
2089
2090 if (elemtype == CGNS_ENUMV(MIXED)) {
2091 int dim;
2092 CGNS_ENUMT(ElementType_t) type;
2093 z->regs[nr].dim = -1;
2094 for (n = 0, ne = 0; ne < nelem; ne++) {
2095 type = (CGNS_ENUMT(ElementType_t))conn[n++];
2096 if (cg_npe (type, &ip) || ip <= 0) {
2097 strcpy(z->regs[nr].errmsg,
2098 "unhandled element type found in MIXED");
2099 break;
2100 }
2101 for (i = 0; i < ip; i++) {
2102 if (conn[n] < 1 || conn[n] > z->nnodes) {
2103 strcpy(z->regs[nr].errmsg, "invalid element index");
2104 break;
2105 }
2106 n++;
2107 }
2108 if (i < ip) break;
2109 dim = element_dimension(type);
2110 if (z->regs[nr].dim < dim) z->regs[nr].dim = dim;
2111 }
2112 }
2113 else if (elemtype == CGNS_ENUMV(NGON_n)) {
2114 z->regs[nr].dim = 2;
2115 for (ne = 0; ne < nelem; ne++) {
2116 ip = (int) conn_offsets[ne+1];
2117 for (n = conn_offsets[ne]; n < conn_offsets[ne+1]; n++) {
2118 if (conn[n] < 1 || conn[n] > z->nnodes) {
2119 strcpy(z->regs[nr].errmsg, "invalid element index");
2120 break;
2121 }
2122 }
2123 if (n < ip) break;
2124 }
2125 }
2126 else {
2127 z->regs[nr].dim = element_dimension(elemtype);
2128 cg_npe (elemtype, &ip);
2129 for (n = 0, ne = 0; ne < nelem; ne++) {
2130 for (i = 0; i < ip; i++) {
2131 if (conn[n] < 1 || conn[n] > z->nnodes) {
2132 strcpy(z->regs[nr].errmsg, "invalid element index");
2133 break;
2134 }
2135 n++;
2136 }
2137 if (i < ip) break;
2138 }
2139 }
2140 if (ne == nelem && z->regs[nr].dim > 0) {
2141
2142 if (z->regs[nr].dim == 1)
2143 edge_elements(&z->regs[nr], conn);
2144 else if (z->regs[nr].dim == 2)
2145 face_elements (&z->regs[nr], conn, conn_offsets);
2146 else
2147 exterior_faces (z, &z->regs[nr], conn, conn_offsets);
2148
2149 #ifndef NO_CUTTING_PLANE
2150 if (z->regs[nr].dim > 1) {
2151 z->regs[nr].elemtype = elemtype;
2152 z->regs[nr].nelems = nelem;
2153 z->regs[nr].elems = conn;
2154 z->regs[nr].elem_offsets = conn_offsets;
2155
2156 /* fix element indexing */
2157 cg_npe (elemtype, &ip);
2158 for (n = 0, ne = 0; ne < nelem; ne++) {
2159 if (elemtype == CGNS_ENUMT(MIXED)) {
2160 nb = (int)conn[n++];
2161 cg_npe ((CGNS_ENUMT(ElementType_t))nb, &ip);
2162 }
2163 else if (elemtype == CGNS_ENUMT(NGON_n)) {
2164 ip = (int)(conn_offsets[ne+1] - conn_offsets[ne]);
2165 }
2166 for (i = 0; i < ip; i++) {
2167 (conn[n])--;
2168 n++;
2169 }
2170 }
2171 }
2172 else
2173 #endif
2174 free (conn);
2175 }
2176 }
2177
2178 /* process NFACE_n sections */
2179
2180 if (haspoly) {
2181 for (ns = 0; ns < nsect; ns++) {
2182 if (z->regs[ns].data[0] != CGNS_ENUMV(NFACE_n)) continue;
2183 zone_message (dspmsg, z->regs[ns].name);
2184 nelem = z->regs[ns].data[2] - z->regs[ns].data[1] + 1;
2185 cg_ElementDataSize (cgnsfn, cgnsbase, cgnszone, ns+1, &elemsize);
2186 conn = (cgsize_t *) MALLOC (funcname, (size_t)elemsize * sizeof(cgsize_t));
2187 conn_offsets = (cgsize_t *) MALLOC (funcname, (size_t) (nelem+1)*sizeof(cgsize_t));
2188 if (cg_poly_elements_read (cgnsfn, cgnsbase, cgnszone, ns+1, conn, conn_offsets, 0)) {
2189 free (conn);
2190 free(conn_offsets);
2191 Tcl_SetResult (interp, (char *)cg_get_error(), TCL_STATIC);
2192 return 1;
2193 }
2194
2195 /* check element indices */
2196
2197 for (n = 0, ne = 0; ne < nelem; ne++) {
2198 ip = (int) conn_offsets[ne+1];
2199 for (n = conn_offsets[ne]; n < conn_offsets[ne+1]; n++) {
2200 if (NULL == find_face(z, abs(conn[n]))) {
2201 strcpy(z->regs[nr].errmsg, "invalid face index");
2202 break;
2203 }
2204 }
2205 if (n < ip) break;
2206 }
2207
2208 if (ne == nelem) {
2209 #ifndef NO_CUTTING_PLANE
2210 polyhedral_faces (z, &z->regs[ns], conn, conn_offsets);
2211 #else
2212 exterior_faces (z, &z->regs[ns], conn, conn_offsets);
2213 free (conn);
2214 #endif
2215 }
2216 }
2217 }
2218
2219 qsort(z->regs, nr, sizeof(Regn), sort_elemsets);
2220
2221 /* 1to1 connectivities */
2222
2223 for (ns = 1; ns <= nints; ns++) {
2224 if (cg_1to1_read (cgnsfn, cgnsbase, cgnszone, ns,
2225 name, d_name, range, d_range, transform)) {
2226 Tcl_SetResult (interp, (char *)cg_get_error(), TCL_STATIC);
2227 return 1;
2228 }
2229 zone_message (dspmsg, name);
2230 strcpy (z->regs[nr].name, name);
2231 z->regs[nr].type = REG_1TO1;
2232 z->regs[nr].data[0] = 2;
2233 z->regs[nr].data[1] = range[0];
2234 z->regs[nr].data[2] = range[1];
2235 strcpy (z->regs[nr].d_name, d_name);
2236 if (unstructured_region (nsect, z->regs, &z->regs[nr],
2237 CGNS_ENUMV(PointRange), 2, range)) z->regs[nr].dim = 2;
2238 nr++;
2239 }
2240
2241 /* general connectivities */
2242
2243 for (ns = 1; ns <= nconns; ns++) {
2244 if (cg_conn_info (cgnsfn, cgnsbase, cgnszone, ns, name,
2245 &location, &type, &ptype, &np, d_name, &d_ztype,
2246 &d_ptype, &datatype, &ie)) {
2247 Tcl_SetResult (interp, (char *)cg_get_error(), TCL_STATIC);
2248 return 1;
2249 }
2250 zone_message (dspmsg, name);
2251 conn = (cgsize_t *) MALLOC (funcname, (size_t)np * sizeof(cgsize_t));
2252 if (cg_conn_read_short (cgnsfn, cgnsbase, cgnszone, ns, conn)) {
2253 free (conn);
2254 Tcl_SetResult (interp, (char *)cg_get_error(), TCL_STATIC);
2255 return 1;
2256 }
2257 strcpy (z->regs[nr].name, name);
2258 z->regs[nr].type = REG_CONN;
2259 z->regs[nr].data[0] = type;
2260 z->regs[nr].data[1] = location;
2261 z->regs[nr].data[2] = ptype;
2262 z->regs[nr].data[3] = np;
2263 if (ptype == CGNS_ENUMT(PointRange)) {
2264 z->regs[nr].data[4] = conn[0];
2265 z->regs[nr].data[5] = conn[1];
2266 }
2267 strcpy (z->regs[nr].d_name, d_name);
2268
2269 if (type != CGNS_ENUMV(Abutting) &&
2270 type != CGNS_ENUMV(Abutting1to1)) {
2271 strcpy(z->regs[nr].errmsg,
2272 "can only handle Abutting or Abutting1to1 currently");
2273 }
2274 else if (ptype != CGNS_ENUMV(PointList) &&
2275 ptype != CGNS_ENUMV(PointRange)) {
2276 strcpy(z->regs[nr].errmsg,
2277 "point set type not PointList or PointRange");
2278 }
2279 else if (location != CGNS_ENUMV(Vertex) &&
2280 location != CGNS_ENUMV(CellCenter) &&
2281 location != CGNS_ENUMV(FaceCenter)) {
2282 strcpy(z->regs[nr].errmsg,
2283 "location not Vertex, CellCenter or FaceCenter");
2284 }
2285 else {
2286 if (location != CGNS_ENUMV(Vertex)) {
2287 ptype = (ptype == CGNS_ENUMV(PointRange) ?
2288 CGNS_ENUMV(ElementRange) : CGNS_ENUMV(ElementList));
2289 }
2290 if (unstructured_region (nsect, z->regs, &z->regs[nr],
2291 ptype, np, conn)) z->regs[nr].dim = 2;
2292 }
2293 free (conn);
2294 nr++;
2295 }
2296
2297 /* holes */
2298
2299 for (ns = 1; ns <= nholes; ns++) {
2300 if (cg_hole_info (cgnsfn, cgnsbase, cgnszone, ns, name,
2301 &location, &ptype, &nsets, &np)) {
2302 Tcl_SetResult (interp, (char *)cg_get_error(), TCL_STATIC);
2303 return 1;
2304 }
2305 conn = (cgsize_t *) MALLOC (funcname, (size_t)(3 * np * nsets) * sizeof(cgsize_t));
2306 if (cg_hole_read (cgnsfn, cgnsbase, cgnszone, ns, conn)) {
2307 free (conn);
2308 Tcl_SetResult (interp, (char *)cg_get_error(), TCL_STATIC);
2309 return 1;
2310 }
2311 strcpy (z->regs[nr].name, name);
2312 z->regs[nr].type = REG_HOLE;
2313 z->regs[nr].data[0] = nsets;
2314 z->regs[nr].data[1] = location;
2315 z->regs[nr].data[2] = ptype;
2316 z->regs[nr].data[3] = np;
2317
2318 if (ptype == CGNS_ENUMV(PointRange)) {
2319 z->regs[nr].data[4] = conn[0];
2320 z->regs[nr].data[5] = conn[1];
2321 }
2322 else if (ptype != CGNS_ENUMV(PointList)) {
2323 strcpy(z->regs[nr].errmsg,
2324 "point set type not PointList or PointRange");
2325 }
2326 else {
2327 if (location == CGNS_ENUMV(Vertex)) {
2328 d_ptype = ptype;
2329 }
2330 else {
2331 d_ptype = (ptype == CGNS_ENUMV(PointRange) ?
2332 CGNS_ENUMV(ElementRange) : CGNS_ENUMV(ElementList));
2333 }
2334 if (unstructured_region (nsect, z->regs, &z->regs[nr],
2335 d_ptype, np, conn)) z->regs[nr].dim = 2;
2336 if (z->regs[nr].dim && nsets > 1 &&
2337 ptype == CGNS_ENUMV(PointRange)) {
2338 Edge *edges = z->regs[nr].edges;
2339 Face **faces = z->regs[nr].faces;
2340 ne = z->regs[nr].nedges;
2341 nf = z->regs[nr].nfaces;
2342 for (ip = 1; ip < nsets; ip++) {
2343 z->regs[nr].nedges = z->regs[nr].nfaces = 0;
2344 is = unstructured_region (nsect, z->regs, &z->regs[nr],
2345 d_ptype, np, &conn[ip*2]);
2346 if (is && z->regs[nr].nedges) {
2347 edges = (Edge *) REALLOC (funcname,
2348 (ne + z->regs[nr].nedges) * sizeof(Edge), edges);
2349 for (i = 0; i < z->regs[nr].nedges; i++) {
2350 edges[ne].nodes[0] = z->regs[nr].edges[i].nodes[0];
2351 edges[ne].nodes[1] = z->regs[nr].edges[i].nodes[1];
2352 ne++;
2353 }
2354 free(z->regs[nr].edges);
2355 }
2356 if (is && z->regs[nr].nfaces) {
2357 faces = (Face **) REALLOC (funcname,
2358 (nf + z->regs[nr].nfaces) * sizeof(Face *), faces);
2359 for (i = 0; i < z->regs[nr].nfaces; i++)
2360 faces[nf++] = z->regs[nr].faces[i];
2361 free(z->regs[nr].faces);
2362 }
2363 }
2364 z->regs[nr].nedges = ne;
2365 z->regs[nr].edges = edges;
2366 z->regs[nr].nfaces = nf;
2367 z->regs[nr].faces = faces;
2368 }
2369 }
2370 free (conn);
2371 nr++;
2372 }
2373
2374 /* boundary conditions */
2375
2376 for (ns = 1; ns <= nbocos; ns++) {
2377 if (cg_boco_info (cgnsfn, cgnsbase, cgnszone, ns, name,
2378 &bctype, &ptype, &np, nrmlindex, &is, &datatype, &nb) ||
2379 cg_boco_gridlocation_read (cgnsfn, cgnsbase, cgnszone, ns,
2380 &location)) {
2381 Tcl_SetResult (interp, (char *)cg_get_error(), TCL_STATIC);
2382 return 1;
2383 }
2384 zone_message (dspmsg, name);
2385 conn = (cgsize_t *) MALLOC (funcname, (size_t)np * sizeof(cgsize_t));
2386 if (cg_boco_read (cgnsfn, cgnsbase, cgnszone, ns, conn, 0)) {
2387 free (conn);
2388 Tcl_SetResult (interp, (char *)cg_get_error(), TCL_STATIC);
2389 return 1;
2390 }
2391 strcpy (z->regs[nr].name, name);
2392 z->regs[nr].type = REG_BOCO;
2393 z->regs[nr].data[0] = bctype;
2394 z->regs[nr].data[1] = location;
2395 z->regs[nr].data[2] = ptype;
2396 z->regs[nr].data[3] = np;
2397 if (ptype == CGNS_ENUMV(PointRange) ||
2398 ptype == CGNS_ENUMV(ElementRange)) {
2399 z->regs[nr].data[4] = conn[0];
2400 z->regs[nr].data[5] = conn[1];
2401 }
2402
2403 if ((ptype == CGNS_ENUMV(PointRange) ||
2404 ptype == CGNS_ENUMV(PointList)) &&
2405 (location == CGNS_ENUMV(FaceCenter) ||
2406 location == CGNS_ENUMV(CellCenter))) {
2407 ptype = (ptype == CGNS_ENUMV(PointRange) ?
2408 CGNS_ENUMV(ElementRange) : CGNS_ENUMV(ElementList));
2409 }
2410 if (unstructured_region (nsect, z->regs, &z->regs[nr],
2411 ptype, np, conn)) z->regs[nr].dim = 2;
2412 free (conn);
2413 nr++;
2414 }
2415
2416 z->nregs = nr;
2417 return 0;
2418 }
2419
2420 /*====================================================================
2421 * find region edges
2422 *====================================================================*/
2423
2424 /*-------------------------------------------------------------------*/
2425
compare_edges(void * v1,void * v2)2426 static int compare_edges (void *v1, void *v2)
2427 {
2428 Edge *e1 = (Edge *)v1;
2429 Edge *e2 = (Edge *)v2;
2430
2431 if (e1->nodes[0] < e1->nodes[1]) {
2432 if (e2->nodes[0] < e2->nodes[1])
2433 return (int)(e1->nodes[0] - e2->nodes[0]);
2434 return (int)(e1->nodes[0] - e2->nodes[1]);
2435 }
2436 if (e2->nodes[0] < e2->nodes[1])
2437 return (int)(e1->nodes[1] - e2->nodes[0]);
2438 return (int)(e1->nodes[1] - e2->nodes[1]);
2439 }
2440
2441 /*-------------------------------------------------------------------*/
2442
hash_edge(void * v)2443 static size_t hash_edge (void *v)
2444 {
2445 Edge *e = (Edge *)v;
2446
2447 return ((size_t)(e->nodes[0] + e->nodes[1]));
2448 }
2449
2450 /*-------------------------------------------------------------------*/
2451
get_edges(void * ve,void * vr)2452 static size_t get_edges (void *ve, void *vr)
2453 {
2454 Edge *e = (Edge *)ve;
2455 Regn *r = (Regn *)vr;
2456
2457 r->edges[r->nedges].nodes[0] = e->nodes[0];
2458 r->edges[r->nedges].nodes[1] = e->nodes[1];
2459 (r->nedges)++;
2460 return 1;
2461 }
2462
2463 /*-------------------------------------------------------------------*/
2464
extract_edges(Regn * r)2465 static void extract_edges (Regn *r)
2466 {
2467 int i, k;
2468 cgsize_t j, n;
2469 size_t ne;
2470 Face *f;
2471 Edge edge, *ep;
2472 HASH *edgehash;
2473 float dot;
2474 static char *funcname = "extract_edges";
2475
2476 if (!r->nfaces) return;
2477 edgehash = HashCreate ((size_t)r->nfaces, compare_edges, hash_edge);
2478 if (NULL == edgehash)
2479 FATAL ("edge hash table creation failed");
2480 for (j = 0; j < r->nfaces; j++) {
2481 f = r->faces[j];
2482 if (f->flags) continue;
2483 for (i = 0, k = f->nnodes-1; i < f->nnodes; k = i++) {
2484 if (f->nodes[i] == f->nodes[k]) continue;
2485 if (f->nodes[i] < f->nodes[k]) {
2486 edge.nodes[0] = f->nodes[i];
2487 edge.nodes[1] = f->nodes[k];
2488 }
2489 else {
2490 edge.nodes[0] = f->nodes[k];
2491 edge.nodes[1] = f->nodes[i];
2492 }
2493 ep = (Edge *) HashFind (edgehash, &edge);
2494 if (NULL == ep) {
2495 ep = (Edge *) MALLOC (funcname, sizeof(Edge));
2496 ep->nodes[0] = edge.nodes[0];
2497 ep->nodes[1] = edge.nodes[1];
2498 ep->id = j;
2499 (void) HashAdd (edgehash, ep);
2500 }
2501 else {
2502 n = ep->id;
2503 dot = r->faces[n]->normal[0] * f->normal[0] +
2504 r->faces[n]->normal[1] * f->normal[1] +
2505 r->faces[n]->normal[2] * f->normal[2];
2506 if (dot > EDGE_ANGLE) {
2507 HashDelete (edgehash, ep);
2508 free (ep);
2509 }
2510 }
2511 }
2512 }
2513
2514 ne = HashSize (edgehash);
2515 if (ne) {
2516 r->nedges = 0;
2517 r->edges = (Edge *) MALLOC (funcname, ne * sizeof(Edge));
2518 HashList (edgehash, get_edges, r);
2519 }
2520 HashDestroy (edgehash, NULL);
2521 }
2522
2523 /*===================================================================
2524 * region manipulation
2525 *===================================================================*/
2526
2527 /*-------------------------------------------------------------------*/
2528
compute_normal(Node n0,Node n1,Node n2,Node n3)2529 static float *compute_normal (Node n0, Node n1, Node n2, Node n3)
2530 {
2531 int j;
2532 double xn, yn, zn, sn;
2533 double d1[3], d2[3];
2534 static float norm[3];
2535
2536 /* triangle */
2537
2538 if (NULL == n3) {
2539 for (j = 0; j < 3; j++) {
2540 d1[j] = n1[j] - n0[j];
2541 d2[j] = n2[j] - n0[j];
2542 }
2543 sn = 0.5;
2544 }
2545
2546 /* quadrilateral */
2547
2548 else {
2549 for (j = 0; j < 3; j++) {
2550 d1[j] = n2[j] - n0[j];
2551 d2[j] = n3[j] - n1[j];
2552 }
2553 sn = 1.0;
2554 }
2555 xn = sn * (d1[1] * d2[2] - d2[1] * d1[2]);
2556 yn = sn * (d1[2] * d2[0] - d2[2] * d1[0]);
2557 zn = sn * (d1[0] * d2[1] - d2[0] * d1[1]);
2558 sn = sqrt (xn * xn + yn * yn + zn * zn);
2559 if (sn == 0.0) sn = 1.0;
2560 norm[0] = (float)(xn / sn);
2561 norm[1] = (float)(yn / sn);
2562 norm[2] = (float)(zn / sn);
2563 return norm;
2564 }
2565
2566 /*-------------------------------------------------------------------*/
2567
face_normal(Zone * z,int nnodes,cgsize_t * nodes)2568 static float *face_normal (Zone *z, int nnodes, cgsize_t *nodes)
2569 {
2570 int i, n;
2571 float *n0, *n1, *n2, *n3;
2572 float *norm, sn;
2573 static float sum[3];
2574
2575 if (nnodes < 3) {
2576 for (i = 0; i < 3; i++)
2577 sum[i] = 0.0;
2578 return sum;
2579 }
2580 if (nnodes <= 4) {
2581 n0 = z->nodes[nodes[0]];
2582 n1 = z->nodes[nodes[1]];
2583 n2 = z->nodes[nodes[2]];
2584 n3 = nnodes == 4 ? z->nodes[nodes[3]] : NULL;
2585 return compute_normal(n0, n1, n2, n3);
2586 }
2587
2588 for (i = 0; i < 3; i++)
2589 sum[i] = 0.0;
2590 n0 = z->nodes[nodes[0]];
2591 n1 = z->nodes[nodes[1]];
2592 for (n = 2; n < nnodes; n++) {
2593 n2 = z->nodes[nodes[n]];
2594 norm = compute_normal(n0, n1, n2, NULL);
2595 for (i = 0; i < 3; i++)
2596 sum[i] += norm[i];
2597 n1 = n2;
2598 }
2599 sn = (float)sqrt(sum[0]*sum[0] + sum[1]*sum[1] + sum[2]*sum[2]);
2600 if (sn == 0.0) sn = 1.0;
2601 for (i = 0; i < 3; i++)
2602 sum[i] /= sn;
2603 return sum;
2604 }
2605
2606 /*-------------------------------------------------------------------*/
2607
region_normals(Zone * z,Regn * r)2608 static void region_normals (Zone *z, Regn *r)
2609 {
2610 int i, n;
2611 Face *f;
2612 float *norm;
2613
2614 for (n = 0; n < r->nfaces; n++) {
2615 f = r->faces[n];
2616 norm = face_normal(z, f->nnodes, f->nodes);
2617 for (i = 0; i < 3; i++)
2618 f->normal[i] = norm[i];
2619 }
2620 }
2621
2622 /*-------------------------------------------------------------------*/
2623
bounding_box(Zone * z,Regn * r)2624 static void bounding_box (Zone *z, Regn *r)
2625 {
2626 int i, j, n;
2627
2628 if (r->nfaces) {
2629 Face *f = r->faces[0];
2630 for (j = 0; j < 3; j++)
2631 r->bbox[j][0] = r->bbox[j][1] = z->nodes[f->nodes[0]][j];
2632 for (n = 0; n < r->nfaces; n++) {
2633 f = r->faces[n];
2634 for (i = 0; i < f->nnodes; i++) {
2635 for (j = 0; j < 3; j++) {
2636 if (r->bbox[j][0] > z->nodes[f->nodes[i]][j])
2637 r->bbox[j][0] = z->nodes[f->nodes[i]][j];
2638 if (r->bbox[j][1] < z->nodes[f->nodes[i]][j])
2639 r->bbox[j][1] = z->nodes[f->nodes[i]][j];
2640 }
2641 }
2642 }
2643 }
2644 else if (r->nedges) {
2645 Edge *e = r->edges;
2646 for (j = 0; j < 3; j++)
2647 r->bbox[j][0] = r->bbox[j][1] = z->nodes[e->nodes[0]][j];
2648 for (n = 0; n < r->nedges; n++, e++) {
2649 for (i = 0; i < 2; i++) {
2650 for (j = 0; j < 3; j++) {
2651 if (r->bbox[j][0] > z->nodes[e->nodes[i]][j])
2652 r->bbox[j][0] = z->nodes[e->nodes[i]][j];
2653 if (r->bbox[j][1] < z->nodes[e->nodes[i]][j])
2654 r->bbox[j][1] = z->nodes[e->nodes[i]][j];
2655 }
2656 }
2657 }
2658 }
2659 else {
2660 for (j = 0; j < 3; j++)
2661 r->bbox[j][0] = r->bbox[j][1] = 0.0;
2662 }
2663 }
2664
2665 /*-------------------------------------------------------------------*/
2666
get_bounds(int all,float bbox[3][2])2667 static void get_bounds (int all, float bbox[3][2])
2668 {
2669 int nz, nr, n, first = 1;
2670
2671 for (nz = 0; nz < nzones; nz++) {
2672 for (nr = 0; nr < zones[nz].nregs; nr++) {
2673 if (zones[nz].regs[nr].nfaces &&
2674 (all || zones[nz].regs[nr].mode)) {
2675 if (first) {
2676 for (n = 0; n < 3; n++) {
2677 bbox[n][0] = zones[nz].regs[nr].bbox[n][0];
2678 bbox[n][1] = zones[nz].regs[nr].bbox[n][1];
2679 }
2680 first = 0;
2681 }
2682 else {
2683 for (n = 0; n < 3; n++) {
2684 if (bbox[n][0] > zones[nz].regs[nr].bbox[n][0])
2685 bbox[n][0] = zones[nz].regs[nr].bbox[n][0];
2686 if (bbox[n][1] < zones[nz].regs[nr].bbox[n][1])
2687 bbox[n][1] = zones[nz].regs[nr].bbox[n][1];
2688 }
2689 }
2690 }
2691 }
2692 }
2693 if (first) {
2694 for (n = 0; n < 3; n++) {
2695 bbox[n][0] = 0.0;
2696 bbox[n][1] = 1.0;
2697 }
2698 }
2699 }
2700
2701 /*-------------------------------------------------------------------*/
2702
draw_outlines(Zone * z,Regn * r)2703 static void draw_outlines (Zone *z, Regn *r)
2704 {
2705 glDisable (GL_LIGHTING);
2706 glShadeModel (GL_FLAT);
2707 glBegin (GL_LINES);
2708 if (r->nedges) {
2709 int ne;
2710 for (ne = 0; ne < r->nedges; ne++) {
2711 glVertex3fv (z->nodes[r->edges[ne].nodes[0]]);
2712 glVertex3fv (z->nodes[r->edges[ne].nodes[1]]);
2713 }
2714 }
2715 else {
2716 glVertex3f (r->bbox[0][0], r->bbox[1][0], r->bbox[2][0]);
2717 glVertex3f (r->bbox[0][1], r->bbox[1][0], r->bbox[2][0]);
2718 glVertex3f (r->bbox[0][1], r->bbox[1][0], r->bbox[2][0]);
2719 glVertex3f (r->bbox[0][1], r->bbox[1][1], r->bbox[2][0]);
2720 glVertex3f (r->bbox[0][1], r->bbox[1][1], r->bbox[2][0]);
2721 glVertex3f (r->bbox[0][0], r->bbox[1][1], r->bbox[2][0]);
2722 glVertex3f (r->bbox[0][0], r->bbox[1][1], r->bbox[2][0]);
2723 glVertex3f (r->bbox[0][0], r->bbox[1][0], r->bbox[2][0]);
2724 glVertex3f (r->bbox[0][0], r->bbox[1][0], r->bbox[2][1]);
2725 glVertex3f (r->bbox[0][1], r->bbox[1][0], r->bbox[2][1]);
2726 glVertex3f (r->bbox[0][1], r->bbox[1][0], r->bbox[2][1]);
2727 glVertex3f (r->bbox[0][1], r->bbox[1][1], r->bbox[2][1]);
2728 glVertex3f (r->bbox[0][1], r->bbox[1][1], r->bbox[2][1]);
2729 glVertex3f (r->bbox[0][0], r->bbox[1][1], r->bbox[2][1]);
2730 glVertex3f (r->bbox[0][0], r->bbox[1][1], r->bbox[2][1]);
2731 glVertex3f (r->bbox[0][0], r->bbox[1][0], r->bbox[2][1]);
2732 glVertex3f (r->bbox[0][0], r->bbox[1][0], r->bbox[2][0]);
2733 glVertex3f (r->bbox[0][0], r->bbox[1][0], r->bbox[2][1]);
2734 glVertex3f (r->bbox[0][1], r->bbox[1][0], r->bbox[2][0]);
2735 glVertex3f (r->bbox[0][1], r->bbox[1][0], r->bbox[2][1]);
2736 glVertex3f (r->bbox[0][1], r->bbox[1][1], r->bbox[2][0]);
2737 glVertex3f (r->bbox[0][1], r->bbox[1][1], r->bbox[2][1]);
2738 glVertex3f (r->bbox[0][0], r->bbox[1][1], r->bbox[2][0]);
2739 glVertex3f (r->bbox[0][0], r->bbox[1][1], r->bbox[2][1]);
2740 }
2741 glEnd ();
2742 }
2743
2744 /*-------------------------------------------------------------------*/
2745
draw_mesh(Zone * z,Regn * r)2746 static void draw_mesh (Zone *z, Regn *r)
2747 {
2748 int nf, nn;
2749 Face *f;
2750
2751 glEnable (GL_LIGHTING);
2752 glShadeModel (GL_FLAT);
2753 glPolygonMode (GL_FRONT_AND_BACK, GL_LINE);
2754 for (nf = 0; nf < r->nfaces; nf++) {
2755 f = r->faces[nf];
2756 if (f->flags || f->nnodes < 2) continue;
2757 if (f->nnodes == 2)
2758 glBegin (GL_LINES);
2759 else if (f->nnodes == 3)
2760 glBegin (GL_TRIANGLES);
2761 else if (f->nnodes == 4)
2762 glBegin (GL_QUADS);
2763 else
2764 glBegin (GL_POLYGON);
2765 glNormal3fv (f->normal);
2766 for (nn = 0; nn < f->nnodes; nn++)
2767 glVertex3fv (z->nodes[f->nodes[nn]]);
2768 glEnd ();
2769 }
2770 }
2771
2772 /*-------------------------------------------------------------------*/
2773
draw_shaded(Zone * z,Regn * r)2774 static void draw_shaded (Zone *z, Regn *r)
2775 {
2776 int nf, nn;
2777 Face *f;
2778
2779 glEnable (GL_LIGHTING);
2780 glShadeModel (GL_FLAT);
2781 glPolygonMode (GL_FRONT_AND_BACK, GL_FILL);
2782 for (nf = 0; nf < r->nfaces; nf++) {
2783 f = r->faces[nf];
2784 if (f->flags || f->nnodes < 3) continue;
2785 if (f->nnodes == 3)
2786 glBegin (GL_TRIANGLES);
2787 else if (f->nnodes == 4)
2788 glBegin (GL_QUADS);
2789 else
2790 glBegin (GL_POLYGON);
2791 glNormal3fv (f->normal);
2792 for (nn = 0; nn < f->nnodes; nn++)
2793 glVertex3fv (z->nodes[f->nodes[nn]]);
2794 glEnd ();
2795 }
2796 }
2797
2798 /*===================================================================
2799 * tcl interface
2800 *===================================================================*/
2801
2802 /*---------- CGNSopen ----------------------------------------------
2803 * open a CGNS file - return bases
2804 *------------------------------------------------------------------*/
2805
CGNSopen(ClientData data,Tcl_Interp * interp,int argc,char ** argv)2806 static int CGNSopen (ClientData data, Tcl_Interp *interp, int argc, char **argv)
2807 {
2808 int fn, nb, idum;
2809 static char buff[33];
2810
2811 if (argc != 2) {
2812 Tcl_SetResult (interp, "usage: CGNSopen filename", TCL_STATIC);
2813 return TCL_ERROR;
2814 }
2815 Tcl_ResetResult (interp);
2816
2817 if (cgnsfn) {
2818 cg_close (cgnsfn);
2819 cgnsfn = 0;
2820 }
2821 if (cg_open (argv[1], CG_MODE_READ, &fn)) {
2822 Tcl_SetResult (interp, (char *)cg_get_error(), TCL_STATIC);
2823 return TCL_ERROR;
2824 }
2825 if (cg_nbases (fn, &nb)) {
2826 Tcl_SetResult (interp, (char *)cg_get_error(), TCL_STATIC);
2827 cg_close (fn);
2828 return TCL_ERROR;
2829 }
2830 if (nb < 1) {
2831 Tcl_SetResult (interp, "no bases defined", TCL_STATIC);
2832 cg_close (fn);
2833 return TCL_ERROR;
2834 }
2835 nbases = nb;
2836 for (nb = 1; nb <= nbases; nb++) {
2837 if (cg_base_read (fn, nb, buff, &idum, &idum)) {
2838 Tcl_SetResult (interp, (char *)cg_get_error(), TCL_STATIC);
2839 cg_close (fn);
2840 return TCL_ERROR;
2841 }
2842 Tcl_AppendElement (interp, buff);
2843 }
2844 cgnsfn = fn;
2845 return TCL_OK;
2846 }
2847
2848 /*---------- CGNSclose ---------------------------------------------
2849 * close the open CGNS file
2850 *------------------------------------------------------------------*/
2851
CGNSclose(ClientData data,Tcl_Interp * interp,int argc,char ** argv)2852 static int CGNSclose (ClientData data, Tcl_Interp *interp, int argc, char **argv)
2853 {
2854 if (cgnsfn && cg_close (cgnsfn)) {
2855 Tcl_SetResult (interp, (char *)cg_get_error(), TCL_STATIC);
2856 return TCL_ERROR;
2857 }
2858 cgnsfn = 0;
2859 free_all ();
2860 Tcl_ResetResult (interp);
2861 return TCL_OK;
2862 }
2863
2864 /*---------- CGNSbase ----------------------------------------------
2865 * set the CGNS base - return zones
2866 *------------------------------------------------------------------*/
2867
CGNSbase(ClientData data,Tcl_Interp * interp,int argc,char ** argv)2868 static int CGNSbase (ClientData data, Tcl_Interp *interp, int argc, char **argv)
2869 {
2870 int base, nz;
2871 cgsize_t sizes[9];
2872 char buff[33];
2873
2874 if (!cgnsfn) {
2875 Tcl_SetResult (interp, "CGNS file not open", TCL_STATIC);
2876 return TCL_ERROR;
2877 }
2878 if (argc != 2) {
2879 Tcl_SetResult (interp, "usage: CGNSbase basenum", TCL_STATIC);
2880 return TCL_ERROR;
2881 }
2882 base = atoi (argv[1]) + 1;
2883 if (base < 1 || base > nbases) {
2884 Tcl_SetResult (interp, "base number out of range", TCL_STATIC);
2885 return TCL_ERROR;
2886 }
2887 Tcl_ResetResult (interp);
2888 cgnsbase = base;
2889 if (cg_base_read (cgnsfn, cgnsbase, BaseName, &CellDim, &PhyDim) ||
2890 cg_nzones (cgnsfn, cgnsbase, &nz)) {
2891 Tcl_SetResult (interp, (char *)cg_get_error(), TCL_STATIC);
2892 return TCL_ERROR;
2893 }
2894 free_all ();
2895 if (CellDim < 2 || CellDim > 3 || PhyDim < 2 || PhyDim > 3) {
2896 Tcl_SetResult (interp, "CellDim and Phydim not 2 or 3", TCL_STATIC);
2897 return TCL_ERROR;
2898 }
2899 nzones = nz;
2900 zones = (Zone *) MALLOC ("CGNSbase", nzones * sizeof(Zone));
2901 for (nz = 1; nz <= nzones; nz++) {
2902 if (cg_zone_read (cgnsfn, cgnsbase, nz, buff, sizes)) {
2903 Tcl_SetResult (interp, (char *)cg_get_error(), TCL_STATIC);
2904 return TCL_ERROR;
2905 }
2906 strcpy (zones[nz-1].name, buff);
2907 Tcl_AppendElement (interp, buff);
2908 }
2909 return TCL_OK;
2910 }
2911
2912 /*---------- CGNSzone ----------------------------------------------
2913 * set the CGNS zone - return regions
2914 *------------------------------------------------------------------*/
2915
CGNSzone(ClientData data,Tcl_Interp * interp,int argc,char ** argv)2916 static int CGNSzone (ClientData data, Tcl_Interp *interp, int argc, char **argv)
2917 {
2918 int i, j, zone, nr, ncoords;
2919 cgsize_t sizes[9], rng[2][3], n, nnodes;
2920 int rind[6];
2921 CGNS_ENUMT(DataType_t) datatype;
2922 CGNS_ENUMT(ZoneType_t) zonetype;
2923 Node *nodes;
2924 float *xyz;
2925 double rad, theta, phi;
2926 char buff[65], coordtype[4];
2927 Zone *z;
2928 #ifdef NO_CUTTING_PLANE
2929 int *tag, nf, nn;
2930 #endif
2931 static char *funcname = "CGNSzone";
2932
2933 if (!cgnsfn) {
2934 Tcl_SetResult (interp, "CGNS file not open", TCL_STATIC);
2935 return TCL_ERROR;
2936 }
2937 if (argc != 2) {
2938 Tcl_SetResult (interp, "usage: CGNSzone zonenum", TCL_STATIC);
2939 return TCL_ERROR;
2940 }
2941 zone = atoi (argv[1]) + 1;
2942 if (zone < 1 || zone > nzones) {
2943 Tcl_SetResult (interp, "zone number out of range", TCL_STATIC);
2944 return TCL_ERROR;
2945 }
2946
2947 if (cg_zone_read (cgnsfn, cgnsbase, zone, buff, sizes) ||
2948 cg_zone_type (cgnsfn, cgnsbase, zone, &zonetype)) {
2949 Tcl_SetResult (interp, (char *)cg_get_error(), TCL_STATIC);
2950 return TCL_ERROR;
2951 }
2952 if (zonetype == CGNS_ENUMV(Structured)) {
2953 for (i = 0; i < CellDim; i++) {
2954 if (sizes[i] < 2) {
2955 Tcl_SetResult (interp, "zone dimension < 2", TCL_STATIC);
2956 return TCL_ERROR;
2957 }
2958 }
2959 }
2960 else if (zonetype != CGNS_ENUMV(Unstructured)) {
2961 Tcl_SetResult (interp, "invalid zone type", TCL_STATIC);
2962 return TCL_ERROR;
2963 }
2964 cgnszone = zone;
2965 z = &zones[zone-1];
2966
2967 /* get number of coordinates */
2968
2969 if (cg_ncoords (cgnsfn, cgnsbase, cgnszone, &ncoords)) {
2970 Tcl_SetResult (interp, (char *)cg_get_error(), TCL_STATIC);
2971 return TCL_ERROR;
2972 }
2973 if (ncoords < PhyDim) {
2974 Tcl_SetResult (interp, "less than PhyDim coordinates", TCL_STATIC);
2975 return TCL_ERROR;
2976 }
2977
2978 /* check for rind */
2979
2980 if (cg_goto (cgnsfn, cgnsbase, "Zone_t", zone,
2981 "GridCoordinates_t", 1, "end")) {
2982 Tcl_SetResult (interp, (char *)cg_get_error(), TCL_STATIC);
2983 return TCL_ERROR;
2984 }
2985 if ((i = cg_rind_read (rind)) != CG_OK) {
2986 if (i != CG_NODE_NOT_FOUND) {
2987 Tcl_SetResult (interp, (char *)cg_get_error(), TCL_STATIC);
2988 return TCL_ERROR;
2989 }
2990 for (i = 0; i < 6; i++)
2991 rind[i] = 0;
2992 }
2993
2994 /* get grid coordinate range */
2995
2996 if (zonetype == CGNS_ENUMV(Structured)) {
2997 for (i = 0; i < 3; i++) {
2998 rng[0][i] = 1;
2999 rng[1][i] = 1;
3000 }
3001 nnodes = 1;
3002 for (i = 0; i < CellDim; i++) {
3003 rng[0][i] = rind[2*i] + 1;
3004 rng[1][i] = rind[2*i] + sizes[i];
3005 nnodes *= sizes[i];
3006 }
3007 }
3008 else {
3009 nnodes = sizes[0] + rind[0] + rind[1];
3010 rng[0][0] = 1;
3011 rng[1][0] = nnodes;
3012 }
3013
3014 /* read the nodes */
3015
3016 strcpy (coordtype, " ");
3017 zone_message ("reading coordinates", NULL);
3018 xyz = (float *) MALLOC (funcname, (size_t)nnodes * sizeof(float));
3019 nodes = (Node *) MALLOC (funcname, (size_t)nnodes * sizeof(Node));
3020 for (i = 1; i <= ncoords; i++) {
3021 if (cg_coord_info (cgnsfn, cgnsbase, cgnszone, i, &datatype, buff) ||
3022 cg_coord_read (cgnsfn, cgnsbase, cgnszone, buff,
3023 CGNS_ENUMV(RealSingle), rng[0], rng[1], xyz)) {
3024 free (xyz);
3025 free (nodes);
3026 Tcl_SetResult (interp, (char *)cg_get_error(), TCL_STATIC);
3027 return TCL_ERROR;
3028 }
3029 if (0 == strcmp (buff, "CoordinateX") ||
3030 0 == strcmp (buff, "CoordinateR"))
3031 j = 0;
3032 else if (0 == strcmp (buff, "CoordinateY") ||
3033 0 == strcmp (buff, "CoordinateTheta"))
3034 j = 1;
3035 else if (0 == strcmp (buff, "CoordinateZ") ||
3036 0 == strcmp (buff, "CoordinatePhi"))
3037 j = 2;
3038 else
3039 continue;
3040 if (coordtype[j] == ' ' || strchr ("XYZ", buff[10]) != NULL)
3041 coordtype[j] = buff[10];
3042 for (n = 0; n < nnodes; n++)
3043 nodes[n][j] = xyz[n];
3044 }
3045 free (xyz);
3046 if (0 == strncmp (coordtype, "RTZ", PhyDim)) {
3047 for (n = 0; n < nnodes; n++) {
3048 rad = nodes[n][0];
3049 theta = nodes[n][1];
3050 nodes[n][0] = (float)(rad * cos (theta));
3051 nodes[n][1] = (float)(rad * sin (theta));
3052 }
3053 }
3054 else if (0 == strcmp (coordtype, "RTP")) {
3055 for (n = 0; n < nnodes; n++) {
3056 rad = nodes[n][0];
3057 theta = nodes[n][1];
3058 phi = nodes[n][2];
3059 nodes[n][0] = (float)(rad * sin (theta) * cos (phi));
3060 nodes[n][1] = (float)(rad * sin (theta) * sin (phi));
3061 nodes[n][2] = (float)(rad * cos (theta));
3062 }
3063 }
3064 else if (strncmp (coordtype, "XYZ", PhyDim)) {
3065 free (nodes);
3066 Tcl_SetResult (interp, "unknown coordinate system", TCL_STATIC);
3067 return TCL_ERROR;
3068 }
3069
3070 z->nnodes = nnodes;
3071 z->nodes = nodes;
3072
3073 /* build regions */
3074
3075 if (zonetype == CGNS_ENUMV(Structured)) {
3076 if (structured_zone (interp, sizes))
3077 return TCL_ERROR;
3078 }
3079 else {
3080 if (unstructured_zone (interp))
3081 return TCL_ERROR;
3082 }
3083
3084 #ifdef NO_CUTTING_PLANE
3085
3086 tag = (int *) MALLOC (funcname, nnodes * sizeof(int));
3087 for (n = 0; n < nnodes; n++)
3088 tag[n] = -1;
3089
3090 /* tag nodes which are actually used */
3091
3092 for (nn = 0, nr = 0; nr < z->nregs; nr++) {
3093 for (nf = 0; nf < z->regs[nr].nfaces; nf++) {
3094 for (n = 0; n < z->regs[nr].faces[nf]->nnodes; n++) {
3095 i = z->regs[nr].faces[nf]->nodes[n];
3096 if (tag[i] < 0)
3097 tag[i] = nn++;
3098 }
3099 }
3100 }
3101
3102 nodes = (Node *) MALLOC (funcname, nn * sizeof(Node));
3103 for (n = 0; n < nnodes; n++) {
3104 if (tag[n] >= 0) {
3105 j = tag[n];
3106 for (i = 0; i < 3; i++)
3107 nodes[j][i] = z->nodes[n][i];
3108 }
3109 }
3110
3111 free(z->nodes);
3112 z->nodes = nodes;
3113 z->nnodes = nn;
3114
3115 /* re-index region faces */
3116
3117 for (nr = 0; nr < z->nregs; nr++) {
3118 for (nf = 0; nf < z->regs[nr].nfaces; nf++) {
3119 for (n = 0; n < z->regs[nr].faces[nf]->nnodes; n++) {
3120 i = z->regs[nr].faces[nf]->nodes[n];
3121 z->regs[nr].faces[nf]->nodes[n] = tag[i];
3122 }
3123 }
3124 }
3125
3126 free(tag);
3127
3128 #endif
3129
3130 /* find region bounding box, edges and normals */
3131
3132 zone_message ("finding normals and edges", NULL);
3133 for (nr = 0; nr < z->nregs; nr++) {
3134 if (z->regs[nr].nfaces) {
3135 bounding_box (z, &z->regs[nr]);
3136 region_normals (z, &z->regs[nr]);
3137 extract_edges (&z->regs[nr]);
3138 }
3139 }
3140
3141 Tcl_ResetResult (interp);
3142 for (nr = 0; nr < z->nregs; nr++) {
3143 switch (z->regs[nr].type) {
3144 case REG_MESH:
3145 strcpy(buff, z->regs[nr].name);
3146 break;
3147 case REG_ELEM:
3148 sprintf(buff, "<Element Sections>/%s", z->regs[nr].name);
3149 break;
3150 case REG_1TO1:
3151 sprintf(buff, "<1to1 Connections>/%s", z->regs[nr].name);
3152 break;
3153 case REG_CONN:
3154 sprintf(buff, "<General Connections>/%s", z->regs[nr].name);
3155 break;
3156 case REG_HOLE:
3157 sprintf(buff, "<Overset Holes>/%s", z->regs[nr].name);
3158 break;
3159 case REG_BOCO:
3160 sprintf(buff, "<Boundary Conditions>/%s", z->regs[nr].name);
3161 break;
3162 case REG_BNDS:
3163 sprintf(buff, "<Mesh Boundaries>/%s", z->regs[nr].name);
3164 break;
3165 }
3166 Tcl_AppendElement (interp, buff);
3167 }
3168 return TCL_OK;
3169 }
3170
3171 /*---------- CGNSsummary -------------------------------------------
3172 * return info summary string
3173 *------------------------------------------------------------------*/
3174
CGNSsummary(ClientData data,Tcl_Interp * interp,int argc,char ** argv)3175 static int CGNSsummary (ClientData data, Tcl_Interp *interp, int argc, char **argv)
3176 {
3177 int n, nz;
3178 char *p, buff[128];
3179 Regn *r;
3180
3181 if (!cgnsfn) {
3182 Tcl_SetResult (interp, "CGNS file not open", TCL_STATIC);
3183 return TCL_ERROR;
3184 }
3185 if (argc < 1 || argc > 3) {
3186 Tcl_SetResult (interp, "usage: CGNSsummary [zone [reg]]", TCL_STATIC);
3187 return TCL_ERROR;
3188 }
3189 Tcl_ResetResult (interp);
3190 if (argc == 1) {
3191 sprintf (buff, "Physical Dim = %d, Cell Dim = %d", PhyDim, CellDim);
3192 Tcl_AppendResult (interp, buff, NULL);
3193 return TCL_OK;
3194 }
3195
3196 nz = atoi (argv[1]);
3197 if (nz < 0 || nz >= nzones) {
3198 Tcl_SetResult (interp, "zone number out of range", TCL_STATIC);
3199 return TCL_ERROR;
3200 }
3201
3202 if (argc == 2) {
3203 cgsize_t sizes[9];
3204 CGNS_ENUMT(ZoneType_t) zonetype;
3205 if (cg_zone_read (cgnsfn, cgnsbase, nz+1, buff, sizes) ||
3206 cg_zone_type (cgnsfn, cgnsbase, nz+1, &zonetype)) {
3207 Tcl_SetResult (interp, (char *)cg_get_error(), TCL_STATIC);
3208 return TCL_ERROR;
3209 }
3210 Tcl_AppendResult (interp, cg_ZoneTypeName(zonetype),
3211 " Zone : ", NULL);
3212 if (zonetype == CGNS_ENUMV(Unstructured)) {
3213 sprintf (buff, "%ld vertices, %ld elements",
3214 (long)sizes[0], (long)sizes[1]);
3215 Tcl_AppendResult (interp, buff, NULL);
3216 }
3217 else {
3218 sprintf (buff, "%ld", (long)sizes[0]);
3219 for (n = 1; n < CellDim; n++) {
3220 p = buff + strlen(buff);
3221 sprintf (p, " x %ld", (long)sizes[n]);
3222 }
3223 Tcl_AppendResult (interp, buff, " vertices", NULL);
3224 }
3225 return TCL_OK;
3226 }
3227
3228 n = atoi (argv[2]);
3229 if (n < 0 || n >= zones[nz].nregs) {
3230 Tcl_SetResult (interp, "region number out of range", TCL_STATIC);
3231 return TCL_ERROR;
3232 }
3233 r = &zones[nz].regs[n];
3234
3235 switch (r->type) {
3236 case REG_MESH:
3237 if (CellDim == 2) {
3238 sprintf (buff, "%ld x %ld",
3239 (long)r->data[0], (long)r->data[1]);
3240 }
3241 else {
3242 sprintf (buff, "%ld x %ld x %ld",
3243 (long)r->data[0], (long)r->data[1], (long)r->data[2]);
3244 }
3245 Tcl_AppendResult (interp, "Structured Mesh : ",
3246 buff, " vertices", NULL);
3247 break;
3248 case REG_ELEM:
3249 sprintf (buff, "%ld", (long)(r->data[2] - r->data[1] + 1));
3250 Tcl_AppendResult (interp, cg_ElementTypeName(r->data[0]),
3251 " Element Set : ", buff, " elements", NULL);
3252 break;
3253 case REG_1TO1:
3254 if (r->data[0] == 2)
3255 sprintf (buff, "%ld", (long)(r->data[2] - r->data[1] + 1));
3256 else if (CellDim == 2) {
3257 sprintf (buff, "%ld x %ld",
3258 (long)(r->data[3] - r->data[1] + 1),
3259 (long)(r->data[4] - r->data[2] + 1));
3260 }
3261 else {
3262 sprintf (buff, "%ld x %ld x %ld",
3263 (long)(r->data[4] - r->data[1] + 1),
3264 (long)(r->data[5] - r->data[2] + 1),
3265 (long)(r->data[6] - r->data[3] + 1));
3266 }
3267 Tcl_AppendResult (interp, "1to1 Connection : PointRange ",
3268 buff, " -> ", r->d_name, NULL);
3269 break;
3270 case REG_CONN:
3271 if (r->data[2] == CGNS_ENUMV(PointList) ||
3272 r->data[2] == CGNS_ENUMV(ElementList))
3273 sprintf (buff, " %ld", (long)r->data[3]);
3274 else if (r->data[3] == 2)
3275 sprintf (buff, " %ld", (long)(r->data[5] - r->data[4] + 1));
3276 else if (CellDim == 2) {
3277 sprintf (buff, " %ld x %ld",
3278 (long)(r->data[6] - r->data[4] + 1),
3279 (long)(r->data[7] - r->data[5] + 1));
3280 }
3281 else {
3282 sprintf (buff, " %ld x %ld x %ld",
3283 (long)(r->data[7] - r->data[4] + 1),
3284 (long)(r->data[8] - r->data[5] + 1),
3285 (long)(r->data[9] - r->data[6] + 1));
3286 }
3287 Tcl_AppendResult (interp,
3288 cg_GridConnectivityTypeName(r->data[0]),
3289 " Connection : ", cg_PointSetTypeName(r->data[2]),
3290 buff, " -> ", r->d_name, NULL);
3291 break;
3292 case REG_HOLE:
3293 if (r->data[2] == CGNS_ENUMV(PointList) ||
3294 r->data[2] == CGNS_ENUMV(ElementList))
3295 sprintf (buff, " %ld", (long)r->data[3]);
3296 else if (r->data[3] == 2)
3297 sprintf (buff, " %ld", (long)(r->data[5] - r->data[4] + 1));
3298 else if (CellDim == 2) {
3299 sprintf (buff, " %ld x %ld",
3300 (long)(r->data[6] - r->data[4] + 1),
3301 (long)(r->data[7] - r->data[5] + 1));
3302 }
3303 else {
3304 sprintf (buff, " %ld x %ld x %ld",
3305 (long)(r->data[7] - r->data[4] + 1),
3306 (long)(r->data[8] - r->data[5] + 1),
3307 (long)(r->data[9] - r->data[6] + 1));
3308 }
3309 Tcl_AppendResult (interp, "Overset Hole : ",
3310 cg_PointSetTypeName(r->data[2]), buff, NULL);
3311 break;
3312 case REG_BOCO:
3313 if (r->data[2] == CGNS_ENUMV(PointList) ||
3314 r->data[2] == CGNS_ENUMV(ElementList))
3315 sprintf (buff, " %ld", (long)r->data[3]);
3316 else if (r->data[3] == 2)
3317 sprintf (buff, " %ld", (long)(r->data[5] - r->data[4] + 1));
3318 else if (CellDim == 2) {
3319 sprintf (buff, " %ld x %ld",
3320 (long)(r->data[6] - r->data[4] + 1),
3321 (long)(r->data[7] - r->data[5] + 1));
3322 }
3323 else {
3324 sprintf (buff, " %ld x %ld x %ld",
3325 (long)(r->data[7] - r->data[4] + 1),
3326 (long)(r->data[8] - r->data[5] + 1),
3327 (long)(r->data[9] - r->data[6] + 1));
3328 }
3329 Tcl_AppendResult (interp, cg_BCTypeName(r->data[0]),
3330 " Boundary Condition : ", cg_PointSetTypeName(r->data[2]),
3331 buff, NULL);
3332 break;
3333 case REG_BNDS:
3334 if (CellDim == 2) {
3335 sprintf (buff, "%ld x %ld",
3336 (long)(r->data[2] - r->data[0] + 1),
3337 (long)(r->data[3] - r->data[1] + 1));
3338 }
3339 else {
3340 sprintf (buff, "%ld x %ld x %ld",
3341 (long)(r->data[3] - r->data[0] + 1),
3342 (long)(r->data[4] - r->data[1] + 1),
3343 (long)(r->data[5] - r->data[2] + 1));
3344 }
3345 Tcl_AppendResult (interp, "Mesh Boundary : ", buff,
3346 " vertices", NULL);
3347 break;
3348 }
3349 return TCL_OK;
3350 }
3351
3352 /*---------- CGNSgetbase -------------------------------------------
3353 * get base properties
3354 *------------------------------------------------------------------*/
3355
CGNSgetbase(ClientData data,Tcl_Interp * interp,int argc,char ** argv)3356 static int CGNSgetbase (ClientData data, Tcl_Interp *interp, int argc, char **argv)
3357 {
3358 char cd[16], pd[16], nz[16];
3359
3360 if (!cgnsfn) {
3361 Tcl_SetResult (interp, "CGNS file not open", TCL_STATIC);
3362 return TCL_ERROR;
3363 }
3364 if (argc != 1) {
3365 Tcl_SetResult (interp, "usage: CGNSgetbase", TCL_STATIC);
3366 return TCL_ERROR;
3367 }
3368 Tcl_ResetResult (interp);
3369
3370 sprintf (pd, "%d", PhyDim);
3371 sprintf (cd, "%d", CellDim);
3372 sprintf (nz, "%d", nzones);
3373 Tcl_AppendResult (interp,
3374 "Base Name : ", BaseName,
3375 "\nPhysical Dim: ", pd,
3376 "\nCell Dim : ", cd,
3377 "\nNumber Zones: ", nz, NULL);
3378 return TCL_OK;
3379 }
3380
3381 /*---------- CGNSgetzone -------------------------------------------
3382 * get zone properties
3383 *------------------------------------------------------------------*/
3384
CGNSgetzone(ClientData data,Tcl_Interp * interp,int argc,char ** argv)3385 static int CGNSgetzone (ClientData data, Tcl_Interp *interp, int argc, char **argv)
3386 {
3387 int n, ndim, zone, cnts[4];
3388 cgsize_t sizes[9];
3389 CGNS_ENUMT(ZoneType_t) zonetype;
3390 char *p, buff[65];
3391 static char *cntname[] = {
3392 "Element Sections",
3393 "1to1 Connections",
3394 "General Connections",
3395 "Boundary Conditions"
3396 };
3397
3398 if (!cgnsfn) {
3399 Tcl_SetResult (interp, "CGNS file not open", TCL_STATIC);
3400 return TCL_ERROR;
3401 }
3402 if (argc != 2) {
3403 Tcl_SetResult (interp, "usage: CGNSgetzone zonenum", TCL_STATIC);
3404 return TCL_ERROR;
3405 }
3406 zone = atoi (argv[1]) + 1;
3407 if (zone < 1 || zone > nzones) {
3408 Tcl_SetResult (interp, "zone number out of range", TCL_STATIC);
3409 return TCL_ERROR;
3410 }
3411
3412 if (cg_zone_read (cgnsfn, cgnsbase, zone, buff, sizes) ||
3413 cg_zone_type (cgnsfn, cgnsbase, zone, &zonetype) ||
3414 cg_nsections (cgnsfn, cgnsbase, zone, &cnts[0]) ||
3415 cg_n1to1 (cgnsfn, cgnsbase, zone, &cnts[1]) ||
3416 cg_nconns (cgnsfn, cgnsbase, zone, &cnts[2]) ||
3417 cg_nbocos (cgnsfn, cgnsbase, zone, &cnts[3])) {
3418 Tcl_SetResult (interp, (char *)cg_get_error(), TCL_STATIC);
3419 return TCL_ERROR;
3420 }
3421 Tcl_ResetResult (interp);
3422
3423 Tcl_AppendResult (interp, "Zone Name : ", buff,
3424 "\nType of Zone : ", cg_ZoneTypeName(zonetype),
3425 "\nVertex Dimensions : ", NULL);
3426
3427 ndim = zonetype == CGNS_ENUMV(Unstructured) ? 1 : CellDim;
3428 sprintf (buff, "%ld", (long)sizes[0]);
3429 for (n = 1; n < ndim; n++) {
3430 p = buff + strlen(buff);
3431 sprintf (p, " x %ld", (long)sizes[n]);
3432 }
3433 Tcl_AppendResult (interp, buff, "\nCell Dimensions : ", NULL);
3434
3435 sprintf (buff, "%ld", (long)sizes[ndim]);
3436 for (n = 1; n < ndim; n++) {
3437 p = buff + strlen(buff);
3438 sprintf (p, " x %ld", (long)sizes[n+CellDim]);
3439 }
3440 Tcl_AppendResult (interp, buff, NULL);
3441
3442 for (n = 0; n < 4; n++) {
3443 sprintf (buff, "\n%-19s: %d", cntname[n], cnts[n]);
3444 Tcl_AppendResult (interp, buff, NULL);
3445 }
3446 return TCL_OK;
3447 }
3448
3449 /*---------- CGNSgetregion -----------------------------------------
3450 * get region properties
3451 *------------------------------------------------------------------*/
3452
CGNSgetregion(ClientData data,Tcl_Interp * interp,int argc,char ** argv)3453 static int CGNSgetregion (ClientData data, Tcl_Interp *interp, int argc, char **argv)
3454 {
3455 int n;
3456 char buff[128];
3457 Zone *z;
3458 Regn *r;
3459
3460 if (!cgnsfn) {
3461 Tcl_SetResult (interp, "CGNS file not open", TCL_STATIC);
3462 return TCL_ERROR;
3463 }
3464 if (argc != 3) {
3465 Tcl_SetResult (interp, "usage: CGNSgetregion zone reg", TCL_STATIC);
3466 return TCL_ERROR;
3467 }
3468 n = atoi (argv[1]);
3469 if (n < 0 || n >= nzones) {
3470 Tcl_SetResult (interp, "zone number out of range", TCL_STATIC);
3471 return TCL_ERROR;
3472 }
3473 z = &zones[n];
3474 n = atoi (argv[2]);
3475 if (n < 0 || n >= z->nregs) {
3476 Tcl_SetResult (interp, "region number out of range", TCL_STATIC);
3477 return TCL_ERROR;
3478 }
3479 r = &z->regs[n];
3480
3481 Tcl_ResetResult (interp);
3482 switch (r->type) {
3483 case REG_MESH:
3484 if (CellDim == 2) {
3485 sprintf (buff, "%ld x %ld",
3486 (long)r->data[0], (long)r->data[1]);
3487 }
3488 else {
3489 sprintf (buff, "%ld x %ld x %ld",
3490 (long)r->data[0], (long)r->data[1], (long)r->data[2]);
3491 }
3492 Tcl_AppendResult (interp,
3493 "Region Name : ", r->name,
3494 "\nType of Region : Structured Mesh",
3495 "\nMesh Dimensions: ", buff, NULL);
3496 break;
3497 case REG_ELEM:
3498 sprintf (buff, "%ld -> %ld",
3499 (long)r->data[1], (long)r->data[2]);
3500 Tcl_AppendResult (interp,
3501 "Region Name : ", r->name,
3502 "\nType of Region : Element Set",
3503 "\nElement Set Type: ", cg_ElementTypeName(r->data[0]),
3504 "\nElement Range : ", buff, NULL);
3505 if (r->data[3] || r->data[4]) {
3506 sprintf (buff, "%ld %ld",
3507 (long)r->data[3], (long)r->data[4]);
3508 Tcl_AppendResult (interp,
3509 "\nRind Elements : ", buff, NULL);
3510 }
3511 break;
3512 case REG_1TO1:
3513 Tcl_AppendResult (interp,
3514 "Region Name : ", r->name,
3515 "\nType of Region: 1to1 Connectivity",
3516 "\nPoint Set Type: PointRange", NULL);
3517 if (r->data[0] == 2) {
3518 sprintf (buff, "%ld -> %ld",
3519 (long)r->data[1], (long)r->data[2]);
3520 Tcl_AppendResult (interp,
3521 "\nIndex Range : ", buff, NULL);
3522 }
3523 else {
3524 sprintf (buff, "%ld -> %ld",
3525 (long)r->data[1], (long)r->data[1+CellDim]);
3526 Tcl_AppendResult (interp,
3527 "\nI Index Range : ", buff, NULL);
3528 sprintf (buff, "%ld -> %ld",
3529 (long)r->data[2], (long)r->data[2+CellDim]);
3530 Tcl_AppendResult (interp,
3531 "\nJ Index Range : ", buff, NULL);
3532 if (CellDim == 3) {
3533 sprintf (buff, "%ld -> %ld",
3534 (long)r->data[3], (long)r->data[6]);
3535 Tcl_AppendResult (interp,
3536 "\nK Index Range : ", buff, NULL);
3537 }
3538 }
3539 Tcl_AppendResult (interp, "\nDonor Zone : ", r->d_name, NULL);
3540 break;
3541 case REG_CONN:
3542 Tcl_AppendResult (interp,
3543 "Region Name : ", r->name,
3544 "\nType of Region : General Connectivity",
3545 "\nConnectivity Type: ",
3546 cg_GridConnectivityTypeName(r->data[0]),
3547 "\nGrid Location : ", cg_GridLocationName(r->data[1]),
3548 "\nPoint Set Type : ", cg_PointSetTypeName(r->data[2]),
3549 NULL);
3550 if (r->data[2] == CGNS_ENUMV(PointList) || r->data[2] == CGNS_ENUMV(ElementList)) {
3551 sprintf (buff, "%ld", (long)r->data[3]);
3552 Tcl_AppendResult (interp, "\nNumber of Points : ", buff, NULL);
3553 }
3554 else if (r->data[3] == 2) {
3555 sprintf (buff, "%ld -> %ld",
3556 (long)r->data[4], (long)r->data[5]);
3557 Tcl_AppendResult (interp, "\nIndex Range : ", buff, NULL);
3558 }
3559 else {
3560 sprintf (buff, "%ld -> %ld",
3561 (long)r->data[4], (long)r->data[4+CellDim]);
3562 Tcl_AppendResult (interp,
3563 "\nI Index Range : ", buff, NULL);
3564 sprintf (buff, "%ld -> %ld",
3565 (long)r->data[5], (long)r->data[5+CellDim]);
3566 Tcl_AppendResult (interp,
3567 "\nJ Index Range : ", buff, NULL);
3568 if (CellDim == 3) {
3569 sprintf (buff, "%ld -> %ld",
3570 (long)r->data[6], (long)r->data[9]);
3571 Tcl_AppendResult (interp,
3572 "\nK Index Range : ", buff, NULL);
3573 }
3574 }
3575 Tcl_AppendResult (interp, "\nDonor Zone : ", r->d_name, NULL);
3576 break;
3577 case REG_HOLE:
3578 sprintf(buff, "%ld", (long)r->data[0]);
3579 Tcl_AppendResult (interp,
3580 "Region Name : ", r->name,
3581 "\nType of Region : Overset Hole",
3582 "\nNumber Sets : ", buff,
3583 "\nGrid Location : ", cg_GridLocationName(r->data[1]),
3584 "\nPoint Set Type : ", cg_PointSetTypeName(r->data[2]),
3585 NULL);
3586 if (r->data[2] == CGNS_ENUMV(PointList) || r->data[2] == CGNS_ENUMV(ElementList)) {
3587 sprintf (buff, "%ld", (long)r->data[3]);
3588 Tcl_AppendResult (interp, "\nNumber of Points : ", buff, NULL);
3589 }
3590 else if (r->data[3] == 2) {
3591 sprintf (buff, "%ld -> %ld",
3592 (long)r->data[4], (long)r->data[5]);
3593 Tcl_AppendResult (interp, "\nIndex Range : ", buff, NULL);
3594 }
3595 else {
3596 sprintf (buff, "%ld -> %ld",
3597 (long)r->data[4], (long)r->data[4+CellDim]);
3598 Tcl_AppendResult (interp,
3599 "\nI Index Range : ", buff, NULL);
3600 sprintf (buff, "%ld -> %ld",
3601 (long)r->data[5], (long)r->data[5+CellDim]);
3602 Tcl_AppendResult (interp,
3603 "\nJ Index Range : ", buff, NULL);
3604 if (CellDim == 3) {
3605 sprintf (buff, "%ld -> %ld",
3606 (long)r->data[6], (long)r->data[9]);
3607 Tcl_AppendResult (interp,
3608 "\nK Index Range : ", buff, NULL);
3609 }
3610 }
3611 break;
3612 case REG_BOCO:
3613 Tcl_AppendResult (interp,
3614 "Region Name : ", r->name,
3615 "\nType of Region : Boundary Condition",
3616 "\nType of BC : ", cg_BCTypeName(r->data[0]),
3617 "\nGrid Location : ", cg_GridLocationName(r->data[1]),
3618 "\nPoint Set Type : ", cg_PointSetTypeName(r->data[2]),
3619 NULL);
3620 if (r->data[2] == CGNS_ENUMV(PointList) ||
3621 r->data[2] == CGNS_ENUMV(ElementList)) {
3622 sprintf (buff, "%ld", (long)r->data[3]);
3623 Tcl_AppendResult (interp, "\nNumber of Points: ", buff, NULL);
3624 }
3625 else if (r->data[3] == 2) {
3626 sprintf (buff, "%ld -> %ld",
3627 (long)r->data[4], (long)r->data[5]);
3628 Tcl_AppendResult (interp, "\nIndex Range : ", buff, NULL);
3629 }
3630 else {
3631 sprintf (buff, "%ld -> %ld",
3632 (long)r->data[4], (long)r->data[4+CellDim]);
3633 Tcl_AppendResult (interp,
3634 "\nI Index Range : ", buff, NULL);
3635 sprintf (buff, "%ld -> %ld",
3636 (long)r->data[5], (long)r->data[5+CellDim]);
3637 Tcl_AppendResult (interp,
3638 "\nJ Index Range : ", buff, NULL);
3639 if (CellDim == 3) {
3640 sprintf (buff, "%ld -> %ld",
3641 (long)r->data[6], (long)r->data[9]);
3642 Tcl_AppendResult (interp,
3643 "\nK Index Range : ", buff, NULL);
3644 }
3645 }
3646 break;
3647 case REG_BNDS:
3648 strcpy (buff, r->name);
3649 Tcl_AppendResult (interp,
3650 "Region Name : ", r->name,
3651 "\nType of Region: Mesh Boundary", NULL);
3652 sprintf (buff, "%ld -> %ld",
3653 (long)r->data[0], (long)r->data[CellDim]);
3654 Tcl_AppendResult (interp,
3655 "\nI Index Range : ", buff, NULL);
3656 sprintf (buff, "%ld -> %ld",
3657 (long)r->data[1], (long)r->data[1+CellDim]);
3658 Tcl_AppendResult (interp,
3659 "\nJ Index Range : ", buff, NULL);
3660 if (CellDim == 3) {
3661 sprintf (buff, "%ld -> %ld",
3662 (long)r->data[2], (long)r->data[5]);
3663 Tcl_AppendResult (interp,
3664 "\nK Index Range : ", buff, NULL);
3665 }
3666 break;
3667 }
3668 return TCL_OK;
3669 }
3670
3671 /*---------- CGNSregiondim -----------------------------------------
3672 * return dimension of a region
3673 *------------------------------------------------------------------*/
3674
CGNSregiondim(ClientData data,Tcl_Interp * interp,int argc,char ** argv)3675 static int CGNSregiondim (ClientData data, Tcl_Interp *interp, int argc, char **argv)
3676 {
3677 int n;
3678 char buff[16];
3679 Zone *z;
3680
3681 if (!cgnsfn) {
3682 Tcl_SetResult (interp, "CGNS file not open", TCL_STATIC);
3683 return TCL_ERROR;
3684 }
3685 if (argc != 3) {
3686 Tcl_SetResult (interp, "usage: CGNSregtype zone reg", TCL_STATIC);
3687 return TCL_ERROR;
3688 }
3689 n = atoi (argv[1]);
3690 if (n < 0 || n >= nzones) {
3691 Tcl_SetResult (interp, "zone number out of range", TCL_STATIC);
3692 return TCL_ERROR;
3693 }
3694 z = &zones[n];
3695 n = atoi (argv[2]);
3696 if (n < 0 || n >= z->nregs) {
3697 Tcl_SetResult (interp, "region number out of range", TCL_STATIC);
3698 return TCL_ERROR;
3699 }
3700 Tcl_ResetResult (interp);
3701 if (!z->regs[n].dim) {
3702 if (*(z->regs[n].errmsg))
3703 Tcl_SetResult (interp, z->regs[n].errmsg, TCL_STATIC);
3704 return TCL_ERROR;
3705 }
3706 sprintf (buff, "%d", z->regs[n].dim);
3707 Tcl_AppendResult (interp, buff, NULL);
3708 return TCL_OK;
3709 }
3710
3711 /*---------- CGNSbounds --------------------------------------------
3712 * get bounding box
3713 *------------------------------------------------------------------*/
3714
transform_bounds(float m[16],float bb[3][2])3715 static void transform_bounds (float m[16], float bb[3][2])
3716 {
3717 int i, j;
3718 float x, y, z, bbox[3][2];
3719
3720 x = m[0] * bb[0][0] + m[4] * bb[1][0] + m[8] * bb[2][0] + m[12];
3721 y = m[1] * bb[0][0] + m[5] * bb[1][0] + m[9] * bb[2][0] + m[13];
3722 z = m[2] * bb[0][0] + m[6] * bb[1][0] + m[10] * bb[2][0] + m[14];
3723 bbox[0][0] = bbox[0][1] = x;
3724 bbox[1][0] = bbox[1][1] = y;
3725 bbox[2][0] = bbox[2][1] = z;
3726
3727 x = m[0] * bb[0][1] + m[4] * bb[1][0] + m[8] * bb[2][0] + m[12];
3728 y = m[1] * bb[0][1] + m[5] * bb[1][0] + m[9] * bb[2][0] + m[13];
3729 z = m[2] * bb[0][1] + m[6] * bb[1][0] + m[10] * bb[2][0] + m[14];
3730 if (bbox[0][0] > x) bbox[0][0] = x;
3731 if (bbox[0][1] < x) bbox[0][1] = x;
3732 if (bbox[1][0] > y) bbox[1][0] = y;
3733 if (bbox[1][1] < y) bbox[1][1] = y;
3734 if (bbox[2][0] > z) bbox[2][0] = z;
3735 if (bbox[2][1] < z) bbox[2][1] = z;
3736
3737 x = m[0] * bb[0][0] + m[4] * bb[1][1] + m[8] * bb[2][0] + m[12];
3738 y = m[1] * bb[0][0] + m[5] * bb[1][1] + m[9] * bb[2][0] + m[13];
3739 z = m[2] * bb[0][0] + m[6] * bb[1][1] + m[10] * bb[2][0] + m[14];
3740 if (bbox[0][0] > x) bbox[0][0] = x;
3741 if (bbox[0][1] < x) bbox[0][1] = x;
3742 if (bbox[1][0] > y) bbox[1][0] = y;
3743 if (bbox[1][1] < y) bbox[1][1] = y;
3744 if (bbox[2][0] > z) bbox[2][0] = z;
3745 if (bbox[2][1] < z) bbox[2][1] = z;
3746
3747 x = m[0] * bb[0][1] + m[4] * bb[1][1] + m[8] * bb[2][0] + m[12];
3748 y = m[1] * bb[0][1] + m[5] * bb[1][1] + m[9] * bb[2][0] + m[13];
3749 z = m[2] * bb[0][1] + m[6] * bb[1][1] + m[10] * bb[2][0] + m[14];
3750 if (bbox[0][0] > x) bbox[0][0] = x;
3751 if (bbox[0][1] < x) bbox[0][1] = x;
3752 if (bbox[1][0] > y) bbox[1][0] = y;
3753 if (bbox[1][1] < y) bbox[1][1] = y;
3754 if (bbox[2][0] > z) bbox[2][0] = z;
3755 if (bbox[2][1] < z) bbox[2][1] = z;
3756
3757 x = m[0] * bb[0][0] + m[4] * bb[1][0] + m[8] * bb[2][1] + m[12];
3758 y = m[1] * bb[0][0] + m[5] * bb[1][0] + m[9] * bb[2][1] + m[13];
3759 z = m[2] * bb[0][0] + m[6] * bb[1][0] + m[10] * bb[2][1] + m[14];
3760 if (bbox[0][0] > x) bbox[0][0] = x;
3761 if (bbox[0][1] < x) bbox[0][1] = x;
3762 if (bbox[1][0] > y) bbox[1][0] = y;
3763 if (bbox[1][1] < y) bbox[1][1] = y;
3764 if (bbox[2][0] > z) bbox[2][0] = z;
3765 if (bbox[2][1] < z) bbox[2][1] = z;
3766
3767 x = m[0] * bb[0][1] + m[4] * bb[1][0] + m[8] * bb[2][1] + m[12];
3768 y = m[1] * bb[0][1] + m[5] * bb[1][0] + m[9] * bb[2][1] + m[13];
3769 z = m[2] * bb[0][1] + m[6] * bb[1][0] + m[10] * bb[2][1] + m[14];
3770 if (bbox[0][0] > x) bbox[0][0] = x;
3771 if (bbox[0][1] < x) bbox[0][1] = x;
3772 if (bbox[1][0] > y) bbox[1][0] = y;
3773 if (bbox[1][1] < y) bbox[1][1] = y;
3774 if (bbox[2][0] > z) bbox[2][0] = z;
3775 if (bbox[2][1] < z) bbox[2][1] = z;
3776
3777 x = m[0] * bb[0][0] + m[4] * bb[1][1] + m[8] * bb[2][1] + m[12];
3778 y = m[1] * bb[0][0] + m[5] * bb[1][1] + m[9] * bb[2][1] + m[13];
3779 z = m[2] * bb[0][0] + m[6] * bb[1][1] + m[10] * bb[2][1] + m[14];
3780 if (bbox[0][0] > x) bbox[0][0] = x;
3781 if (bbox[0][1] < x) bbox[0][1] = x;
3782 if (bbox[1][0] > y) bbox[1][0] = y;
3783 if (bbox[1][1] < y) bbox[1][1] = y;
3784 if (bbox[2][0] > z) bbox[2][0] = z;
3785 if (bbox[2][1] < z) bbox[2][1] = z;
3786
3787 x = m[0] * bb[0][1] + m[4] * bb[1][1] + m[8] * bb[2][1] + m[12];
3788 y = m[1] * bb[0][1] + m[5] * bb[1][1] + m[9] * bb[2][1] + m[13];
3789 z = m[2] * bb[0][1] + m[6] * bb[1][1] + m[10] * bb[2][1] + m[14];
3790 if (bbox[0][0] > x) bbox[0][0] = x;
3791 if (bbox[0][1] < x) bbox[0][1] = x;
3792 if (bbox[1][0] > y) bbox[1][0] = y;
3793 if (bbox[1][1] < y) bbox[1][1] = y;
3794 if (bbox[2][0] > z) bbox[2][0] = z;
3795 if (bbox[2][1] < z) bbox[2][1] = z;
3796
3797 for (i = 0; i < 3; i++)
3798 for (j = 0; j < 2; j++)
3799 bb[i][j] = bbox[i][j];
3800 }
3801
3802 /*-------------------------------------------------------------------*/
3803
CGNSbounds(ClientData data,Tcl_Interp * interp,int argc,char ** argv)3804 static int CGNSbounds (ClientData data, Tcl_Interp *interp, int argc, char **argv)
3805 {
3806 float bbox[3][2], matrix[16];
3807 int n, all = 0;
3808 CONST char **args;
3809 char sbb[65];
3810
3811 if (argc > 1) all = atoi(argv[1]);
3812 get_bounds (all, bbox);
3813 if (argc > 2) {
3814 if (TCL_OK != Tcl_SplitList (interp, argv[2], &n, &args))
3815 return TCL_ERROR;
3816 for (n = 0; n < 16; n++)
3817 matrix[n] = (float) atof (args[n]);
3818 Tcl_Free ((char *)args);
3819 transform_bounds (matrix, bbox);
3820 }
3821 Tcl_ResetResult (interp);
3822 for (n = 0; n < 3; n++) {
3823 sprintf (sbb, "%f %f", bbox[n][0], bbox[n][1]);
3824 Tcl_AppendElement (interp, sbb);
3825 }
3826 return TCL_OK;
3827 }
3828
3829 /*---------- OGLregion ---------------------------------------------
3830 * create OGL display list for region
3831 *------------------------------------------------------------------*/
3832
OGLregion(ClientData data,Tcl_Interp * interp,int argc,char ** argv)3833 static int OGLregion (ClientData data, Tcl_Interp *interp, int argc, char **argv)
3834 {
3835 int zone, regn, nc;
3836 CONST char **args;
3837 Zone *z;
3838 Regn *r;
3839 static char slist[17];
3840
3841 if (argc != 5) {
3842 Tcl_SetResult (interp, "usage: OGLregion zone region mode color",
3843 TCL_STATIC);
3844 return TCL_ERROR;
3845 }
3846 zone = atoi (argv[1]);
3847 if (zone < 0 || zone >= nzones) {
3848 Tcl_SetResult (interp, "zone number out of range", TCL_STATIC);
3849 return TCL_ERROR;
3850 }
3851 z = &zones[zone];
3852 regn = atoi (argv[2]);
3853 if (regn < 0 || regn >= z->nregs) {
3854 Tcl_SetResult (interp, "region number out of range", TCL_STATIC);
3855 return TCL_ERROR;
3856 }
3857 r = &z->regs[regn];
3858
3859 if (r->nfaces || r->nedges) {
3860 r->mode = atoi (argv[3]);
3861 if (TCL_OK != Tcl_SplitList (interp, argv[4], &nc, &args))
3862 return TCL_ERROR;
3863 if (nc != 3) {
3864 Tcl_Free ((char *)args);
3865 Tcl_SetResult (interp, "invalid color", TCL_STATIC);
3866 return TCL_ERROR;
3867 }
3868 for (nc = 0; nc < 3; nc++)
3869 r->color[nc] = (float)atof (args[nc]);
3870 r->color[3] = 1.0;
3871 Tcl_Free ((char *)args);
3872
3873 if (!r->dlist) r->dlist = glGenLists (1);
3874 glNewList (r->dlist, GL_COMPILE);
3875 glColor3fv (r->color);
3876 glMaterialfv (GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, r->color);
3877 if (r->nfaces) {
3878 switch (r->mode) {
3879 case 1:
3880 draw_outlines (z, r);
3881 break;
3882 case 2:
3883 draw_mesh (z, r);
3884 break;
3885 case 3:
3886 draw_shaded (z, r);
3887 break;
3888 default:
3889 r->mode = 0;
3890 break;
3891 }
3892 }
3893 else if (r->mode < 1 || r->mode > 3) {
3894 r->mode = 0;
3895 }
3896 else {
3897 draw_outlines (z, r);
3898 }
3899 glEndList ();
3900 }
3901
3902 sprintf (slist, "%d", r->dlist);
3903 Tcl_SetResult (interp, slist, TCL_STATIC);
3904 return TCL_OK;
3905 }
3906
3907 /*---------- OGLaxis -----------------------------------------------
3908 * create OGL display list for axis
3909 *------------------------------------------------------------------*/
3910
3911 #define CHAR_W 8
3912 #define CHAR_H 13
3913
3914 static GLubyte x_raster[] = {
3915 0x00, 0x00, 0xc3, 0x66, 0x66, 0x3c, 0x3c, 0x18,
3916 0x3c, 0x3c, 0x66, 0x66, 0xc3
3917 };
3918 static GLubyte y_raster[] = {
3919 0x00, 0x00, 0x18, 0x18, 0x18, 0x18, 0x18, 0x18,
3920 0x3c, 0x3c, 0x66, 0x66, 0xc3
3921 };
3922 static GLubyte z_raster[] = {
3923 0x00, 0x00, 0xff, 0xc0, 0xc0, 0x60, 0x30, 0x7e,
3924 0x0c, 0x06, 0x03, 0x03, 0xff
3925 };
3926
3927 /*-------------------------------------------------------------------*/
3928
OGLaxis(ClientData data,Tcl_Interp * interp,int argc,char ** argv)3929 static int OGLaxis (ClientData data, Tcl_Interp *interp, int argc, char **argv)
3930 {
3931 int vis;
3932 float bbox[3][2];
3933 static char slist[17];
3934
3935 if (argc < 2 || argc > 3) {
3936 Tcl_SetResult (interp, "usage: OGLaxis visible [bounds]",
3937 TCL_STATIC);
3938 return TCL_ERROR;
3939 }
3940 vis = atoi (argv[1]);
3941 if (!AxisDL) AxisDL = glGenLists (1);
3942
3943 glNewList (AxisDL, GL_COMPILE);
3944 if (vis) {
3945 if (argc == 3) {
3946 int nb, n = 0;
3947 CONST char **args;
3948 if (TCL_OK != Tcl_SplitList (interp, argv[2], &nb, &args))
3949 return TCL_ERROR;
3950 if (nb == 3) {
3951 for (n = 0; n < nb; n++) {
3952 if (sscanf (args[n], "%f %f", &bbox[n][0], &bbox[n][1]) != 2)
3953 break;
3954 }
3955 }
3956 Tcl_Free ((char *)args);
3957 if (n != 3) {
3958 Tcl_SetResult (interp, "invalid bounding box", TCL_STATIC);
3959 return TCL_ERROR;
3960 }
3961 }
3962 else
3963 get_bounds (0, bbox);
3964 glLineWidth (3.0);
3965 glDisable (GL_LIGHTING);
3966 glShadeModel (GL_FLAT);
3967 glBegin (GL_LINES);
3968 glColor3f (1.0, 0.0, 0.0);
3969 glVertex3f (bbox[0][0], bbox[1][0], bbox[2][0]);
3970 glVertex3f (bbox[0][1], bbox[1][0], bbox[2][0]);
3971 glColor3f (1.0, 1.0, 0.0);
3972 glVertex3f (bbox[0][0], bbox[1][0], bbox[2][0]);
3973 glVertex3f (bbox[0][0], bbox[1][1], bbox[2][0]);
3974 glColor3f (0.0, 1.0, 0.0);
3975 glVertex3f (bbox[0][0], bbox[1][0], bbox[2][0]);
3976 glVertex3f (bbox[0][0], bbox[1][0], bbox[2][1]);
3977 glEnd ();
3978 glLineWidth (1.0);
3979
3980 glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
3981 glColor3f (1.0, 0.0, 0.0);
3982 glRasterPos3f (bbox[0][1], bbox[1][0], bbox[2][0]);
3983 glBitmap (CHAR_W, CHAR_H, -0.5 * (float)CHAR_W,
3984 0.5 * (float)CHAR_H, 0.0, 0.0, x_raster);
3985 glColor3f (1.0, 1.0, 0.0);
3986 glRasterPos3f (bbox[0][0], bbox[1][1], bbox[2][0]);
3987 glBitmap (CHAR_W, CHAR_H, -0.5 * (float)CHAR_W,
3988 0.5 * (float)CHAR_H, 0.0, 0.0, y_raster);
3989 glColor3f (0.0, 1.0, 0.0);
3990 glRasterPos3f (bbox[0][0], bbox[1][0], bbox[2][1]);
3991 glBitmap (CHAR_W, CHAR_H, -0.5 * (float)CHAR_W,
3992 0.5 * (float)CHAR_H, 0.0, 0.0, z_raster);
3993 }
3994 glEndList ();
3995
3996 sprintf (slist, "%d", AxisDL);
3997 Tcl_SetResult (interp, slist, TCL_STATIC);
3998 return TCL_OK;
3999 }
4000
4001 /*-------------------------------------------------------------------*/
4002
OGLcolor(ClientData data,Tcl_Interp * interp,int argc,char ** argv)4003 static int OGLcolor (ClientData data, Tcl_Interp *interp, int argc, char **argv)
4004 {
4005 int index, i, j;
4006 double r, g, b;
4007 double h, v, s;
4008 static char color[256];
4009 static double huemap[12] = {0,1,2,3,4,5,0.5,1.25,2.65,3.4,4.5,5.5};
4010
4011 if (argc != 2) {
4012 Tcl_SetResult (interp, "usage: OGLcolor index",
4013 TCL_STATIC);
4014 return TCL_ERROR;
4015 }
4016 index = abs(atoi(argv[1])) % 132;
4017 h = huemap[index % 12];
4018 i = (int)h;
4019 h -= (double)i;
4020 j = index / 12;
4021 if ((j % 2) == 0) {
4022 v = 1.0;
4023 s = 1.0 - sqrt((double)j / 22.0);
4024 }
4025 else {
4026 v = 1.0 - sqrt((double)j / 44.0);
4027 s = 1.0;
4028 }
4029 r = g = b = 0.0;
4030 switch (i) {
4031 case 6:
4032 h = 0.0;
4033 /* FALLTHRU */
4034 case 0:
4035 r = v;
4036 g = v * (1.0 - (s * (1.0 - h)));
4037 b = v * (1.0 - s);
4038 break;
4039 case 1:
4040 r = v * (1.0 - (s * h));
4041 g = v;
4042 b = v * (1.0 - s);
4043 break;
4044 case 2:
4045 r = v * (1.0 - s);
4046 g = v;
4047 b = v * (1.0 - (s * (1.0 - h)));
4048 break;
4049 case 3:
4050 r = v * (1.0 - s);
4051 g = v * (1.0 - (s * h));
4052 b = v;
4053 break;
4054 case 4:
4055 r = v * (1.0 - (s * (1.0 - h)));
4056 g = v * (1.0 - s);
4057 b = v;
4058 break;
4059 case 5:
4060 r = v;
4061 g = v * (1.0 - s);
4062 b = v * (1.0 - (s * h));
4063 break;
4064 }
4065
4066 if (r < 0.0) r = 0.0;
4067 if (r > 1.0) r = 1.0;
4068 if (g < 0.0) g = 0.0;
4069 if (g > 1.0) g = 1.0;
4070 if (b < 0.0) b = 0.0;
4071 if (b > 1.0) b = 1.0;
4072
4073 sprintf(color, "%g %g %g", r, g, b);
4074 Tcl_SetResult (interp, color, TCL_STATIC);
4075 return TCL_OK;
4076 }
4077
4078 /*==============================================================
4079 * cutting plane routines
4080 *==============================================================*/
4081
4082 #ifndef NO_CUTTING_PLANE
4083
4084 /*------------------------------------------------------------------*/
4085
init_cutplane(float plane[4])4086 static void init_cutplane (float plane[4])
4087 {
4088 int nz, nr;
4089
4090 for (nz = 0; nz < nzones; nz++) {
4091 for (nr = 0; nr < zones[nz].nregs; nr++) {
4092 if (zones[nz].regs[nr].cut.nelems) {
4093 free (zones[nz].regs[nr].cut.elems);
4094 zones[nz].regs[nr].cut.nelems = 0;
4095 }
4096 if (zones[nz].regs[nr].cut.nedges) {
4097 free (zones[nz].regs[nr].cut.edges);
4098 zones[nz].regs[nr].cut.nedges = 0;
4099 }
4100 }
4101 }
4102 if (cutplane.nnodes) free (cutplane.nodes);
4103 cutplane.nelems = 0;
4104 cutplane.nedges = 0;
4105 cutplane.nnodes = 0;
4106 for (nr = 0; nr < 4; nr++)
4107 cutplane.plane[nr] = plane[nr];
4108 }
4109
4110 /*------------------------------------------------------------------*/
4111
classify_element(Zone * z,int nnodes,cgsize_t * nodeid)4112 static int classify_element (Zone *z, int nnodes, cgsize_t *nodeid)
4113 {
4114 int n, index = 0;
4115 int mask = (1 << (nnodes - 1)) - 1;
4116 float *node;
4117 double s;
4118
4119 for (n = 0; n < nnodes; n++) {
4120 node = z->nodes[nodeid[n]];
4121 s = node[0] * cutplane.plane[0] + node[1] * cutplane.plane[1] +
4122 node[2] * cutplane.plane[2] - cutplane.plane[3];
4123 if (s >= 0.0) index |= (1 << n);
4124 }
4125 if (index > mask) index ^= ((mask << 1) | 1);
4126 return index;
4127 }
4128
4129 /*------------------------------------------------------------------*/
4130
classify_polygon(Zone * z,int nnodes,cgsize_t * nodeid)4131 static int classify_polygon (Zone *z, int nnodes, cgsize_t *nodeid)
4132 {
4133 int n, start, diff;
4134 float *node;
4135 double s;
4136
4137 node = z->nodes[nodeid[0]];
4138 s = node[0] * cutplane.plane[0] + node[1] * cutplane.plane[1] +
4139 node[2] * cutplane.plane[2] - cutplane.plane[3];
4140 start = s >= 0.0 ? 1 : -1;
4141
4142 for (n = 1; n < nnodes; n++) {
4143 node = z->nodes[nodeid[n]];
4144 s = node[0] * cutplane.plane[0] + node[1] * cutplane.plane[1] +
4145 node[2] * cutplane.plane[2] - cutplane.plane[3];
4146 diff = start * (s >= 0.0 ? 1 : -1);
4147 if (diff < 0) return 1;
4148 }
4149 return 0;
4150 }
4151
4152 /*------------------------------------------------------------------*/
4153
find_elements()4154 static cgsize_t find_elements ()
4155 {
4156 #define ELEM_INC 50
4157 int nz, nnodes, nn, nr, nf;
4158 cgsize_t n, ne, maxelems, nelems, *elems;
4159 CGNS_ENUMT(ElementType_t) type;
4160 Zone *z;
4161 Regn *r;
4162 Face *f;
4163
4164 for (nz = 0; nz < nzones; nz++) {
4165 z = &zones[nz];
4166 for (nr = 0; nr < z->nregs; nr++) {
4167 r = &z->regs[nr];
4168 if (r->dim < 2 || (r->mode == 0 && !ignorevis)) continue;
4169 type = r->elemtype;
4170 cg_npe(type, &nnodes);
4171 maxelems = nelems = 0;
4172 elems = NULL;
4173
4174 if (type == CGNS_ENUMV(NGON_n)) {
4175 for (n = 0, ne = 0; ne < r->nelems; ne++) {
4176 nn = r->elem_offsets[ne+1] - r->elem_offsets[ne];
4177 if (nn > 2 &&
4178 classify_polygon(z, nn, &r->elems[n])) {
4179 if (nelems >= maxelems) {
4180 maxelems += ELEM_INC;
4181 elems = (cgsize_t *) REALLOC ("find_elements",
4182 (size_t)maxelems * sizeof(cgsize_t), elems);
4183 }
4184 elems[nelems++] = ne;
4185 }
4186 n += nn;
4187 }
4188 }
4189 else if (type == CGNS_ENUMV(NFACE_n)) {
4190 for (n = 0, ne = 0; ne < r->nelems; ne++) {
4191 nf = r->elem_offsets[ne+1] - r->elem_offsets[ne];
4192 for (nn = 0; nn < nf; nn++) {
4193 f = r->poly[abs(r->elems[n+nn])-1];
4194 if (f->nnodes > 2 &&
4195 classify_polygon(z, f->nnodes, f->nodes)) {
4196 if (nelems >= maxelems) {
4197 maxelems += ELEM_INC;
4198 elems = (cgsize_t *) REALLOC ("find_elements",
4199 (size_t)maxelems * sizeof(cgsize_t), elems);
4200 }
4201 elems[nelems++] = ne;
4202 break;
4203 }
4204 }
4205 n += nf;
4206 }
4207 }
4208 else {
4209 for (n = 0, ne = 0; ne < r->nelems; ne++) {
4210 if (r->elemtype == CGNS_ENUMV(MIXED)) {
4211 type = (CGNS_ENUMT(ElementType_t))r->elems[n++];
4212 cg_npe(type, &nnodes);
4213 }
4214 switch (type) {
4215 case CGNS_ENUMV(TRI_3):
4216 case CGNS_ENUMV(TRI_6):
4217 case CGNS_ENUMV(TRI_9):
4218 case CGNS_ENUMV(TRI_10):
4219 case CGNS_ENUMV(TRI_12):
4220 case CGNS_ENUMV(TRI_15):
4221 nn = 3;
4222 break;
4223 case CGNS_ENUMV(QUAD_4):
4224 case CGNS_ENUMV(QUAD_8):
4225 case CGNS_ENUMV(QUAD_9):
4226 case CGNS_ENUMV(QUAD_12):
4227 case CGNS_ENUMV(QUAD_16):
4228 case CGNS_ENUMV(QUAD_P4_16):
4229 case CGNS_ENUMV(QUAD_25):
4230 case CGNS_ENUMV(TETRA_4):
4231 case CGNS_ENUMV(TETRA_10):
4232 case CGNS_ENUMV(TETRA_16):
4233 case CGNS_ENUMV(TETRA_20):
4234 case CGNS_ENUMV(TETRA_22):
4235 case CGNS_ENUMV(TETRA_34):
4236 case CGNS_ENUMV(TETRA_35):
4237 nn = 4;
4238 break;
4239 case CGNS_ENUMV(PYRA_5):
4240 case CGNS_ENUMV(PYRA_13):
4241 case CGNS_ENUMV(PYRA_14):
4242 case CGNS_ENUMV(PYRA_21):
4243 case CGNS_ENUMV(PYRA_29):
4244 case CGNS_ENUMV(PYRA_30):
4245 case CGNS_ENUMV(PYRA_P4_29):
4246 case CGNS_ENUMV(PYRA_50):
4247 case CGNS_ENUMV(PYRA_55):
4248 nn = 5;
4249 break;
4250 case CGNS_ENUMV(PENTA_6):
4251 case CGNS_ENUMV(PENTA_15):
4252 case CGNS_ENUMV(PENTA_18):
4253 case CGNS_ENUMV(PENTA_24):
4254 case CGNS_ENUMV(PENTA_38):
4255 case CGNS_ENUMV(PENTA_40):
4256 case CGNS_ENUMV(PENTA_33):
4257 case CGNS_ENUMV(PENTA_66):
4258 case CGNS_ENUMV(PENTA_75):
4259 nn = 6;
4260 break;
4261 case CGNS_ENUMV(HEXA_8):
4262 case CGNS_ENUMV(HEXA_20):
4263 case CGNS_ENUMV(HEXA_27):
4264 case CGNS_ENUMV(HEXA_32):
4265 case CGNS_ENUMV(HEXA_56):
4266 case CGNS_ENUMV(HEXA_64):
4267 case CGNS_ENUMV(HEXA_44):
4268 case CGNS_ENUMV(HEXA_98):
4269 case CGNS_ENUMV(HEXA_125):
4270 nn = 8;
4271 break;
4272 default:
4273 nn = 0;
4274 break;
4275 }
4276 if (nn && classify_element(z, nn, &r->elems[n])) {
4277 if (nelems >= maxelems) {
4278 maxelems += ELEM_INC;
4279 elems = (cgsize_t *) REALLOC ("find_elements",
4280 (size_t)maxelems * sizeof(cgsize_t), elems);
4281 }
4282 if (r->elemtype == CGNS_ENUMV(MIXED))
4283 elems[nelems] = n - 1;
4284 else
4285 elems[nelems] = n;
4286 nelems++;
4287 }
4288 n += nnodes;
4289 }
4290 }
4291 r->cut.nelems = nelems;
4292 r->cut.elems = elems;
4293 cutplane.nelems += nelems;
4294 }
4295 }
4296
4297 return cutplane.nelems;
4298 }
4299
4300 /*-------------------------------------------------------------------*/
4301
compare_cut_node(void * v1,void * v2)4302 static int compare_cut_node (void *v1, void *v2)
4303 {
4304 int i;
4305 CutNode *c1 = (CutNode *)v1;
4306 CutNode *c2 = (CutNode *)v2;
4307
4308 for (i = 0; i < 3; i++) {
4309 if (c1->nodes[i] != c2->nodes[i])
4310 return (int)(c1->nodes[i] - c2->nodes[i]);
4311 }
4312 return 0;
4313 }
4314
4315 /*-------------------------------------------------------------------*/
4316
hash_cut_node(void * v)4317 static size_t hash_cut_node (void *v)
4318 {
4319 CutNode *c = (CutNode *)v;
4320
4321 return ((size_t)(c->nodes[0] + c->nodes[1] + c->nodes[2]));
4322 }
4323
4324 /*-------------------------------------------------------------------*/
4325
get_cut_node(void * v)4326 static void get_cut_node (void *v)
4327 {
4328 int i;
4329 CutNode *c = (CutNode *)v;
4330 Zone *z = &zones[c->nodes[0]];
4331 float *n, *n1, *n2;
4332
4333 n = cutplane.nodes[c->id];
4334 n1 = z->nodes[c->nodes[1]];
4335 n2 = z->nodes[c->nodes[2]];
4336
4337 for (i = 0; i < 3; i++)
4338 n[i] = n1[i] + c->ratio * (n2[i] - n1[i]);
4339 free (c);
4340 }
4341
4342 /*-------------------------------------------------------------------*/
4343
get_cut_edge(void * ve,void * vc)4344 static size_t get_cut_edge (void *ve, void *vc)
4345 {
4346 Edge *e = (Edge *)ve;
4347 CutData *c = (CutData *)vc;
4348
4349 c->edges[c->nedges].nodes[0] = e->nodes[0];
4350 c->edges[c->nedges].nodes[1] = e->nodes[1];
4351 (c->nedges)++;
4352 return 1;
4353 }
4354
4355 /*----- tri elements -----*/
4356
4357 #define TRI_SIZE 3
4358 #define TRI_EDGES 3
4359
4360 static int triCuts[TRI_SIZE+1][4] = {
4361 {0},
4362 {2,0,2, 0},
4363 {2,0,1, 0},
4364 {2,1,2, 0}
4365 };
4366 static int triEdges[TRI_EDGES][2] = {
4367 {0,1},
4368 {1,2},
4369 {2,0}
4370 };
4371
4372 /*----- quad elements -----*/
4373
4374 #define QUAD_SIZE 7
4375 #define QUAD_EDGES 4
4376
4377 static int quadCuts[QUAD_SIZE+1][4] = {
4378 {0},
4379 {2,0,3, 0},
4380 {2,0,1, 0},
4381 {2,1,3, 0},
4382 {2,1,2, 0},
4383 {2,0,3, 0},
4384 {2,0,2, 0},
4385 {2,2,3, 0}
4386 };
4387 static int quadEdges[QUAD_EDGES][2] = {
4388 {0,1},
4389 {1,2},
4390 {2,3},
4391 {3,0}
4392 };
4393
4394 /*----- tet elements -----*/
4395
4396 #define TET_SIZE 7
4397 #define TET_EDGES 6
4398
4399 static int tetCuts[TET_SIZE+1][6] = {
4400 {0},
4401 {3,0,3,2, 0},
4402 {3,0,1,4, 0},
4403 {4,1,4,3,2, 0},
4404 {3,1,2,5, 0},
4405 {4,0,3,5,1, 0},
4406 {4,0,2,5,4, 0},
4407 {3,3,5,4, 0}
4408 };
4409 static int tetEdges[TET_EDGES][2] = {
4410 {0,1},
4411 {1,2},
4412 {2,0},
4413 {0,3},
4414 {1,3},
4415 {2,3}
4416 };
4417
4418 /*----- pyramid elements -----*/
4419
4420 #define PYR_SIZE 15
4421 #define PYR_EDGES 8
4422
4423 static int pyrCuts[PYR_SIZE+1][9] = {
4424 {0},
4425 {3,0,4,3, 0},
4426 {3,0,1,5, 0},
4427 {4,1,5,4,3, 0},
4428 {3,1,2,6, 0},
4429 {3,0,4,3, 3,1,2,6, 0},
4430 {4,0,2,6,5, 0},
4431 {5,3,2,6,5,4, 0},
4432 {3,2,3,7, 0},
4433 {4,0,4,7,2, 0},
4434 {3,2,3,7, 3,0,1,5, 0},
4435 {5,1,5,4,7,2, 0},
4436 {4,1,3,7,6, 0},
4437 {5,0,4,7,6,1, 0},
4438 {5,0,3,7,6,5, 0},
4439 {4,4,7,6,5, 0},
4440 };
4441 static int pyrEdges[PYR_EDGES][2] = {
4442 {0,1},
4443 {1,2},
4444 {2,3},
4445 {3,0},
4446 {0,4},
4447 {1,4},
4448 {2,4},
4449 {3,4}
4450 };
4451
4452 /*----- wedge elements -----*/
4453
4454 #define WDG_SIZE 31
4455 #define WDG_EDGES 9
4456
4457 static int wdgCuts[WDG_SIZE+1][10] = {
4458 {0},
4459 {3,0,3,2, 0},
4460 {3,0,1,4, 0},
4461 {4,1,4,3,2, 0},
4462 {3,1,2,5, 0},
4463 {4,0,3,5,1, 0},
4464 {4,0,2,5,4, 0},
4465 {3,3,5,4, 0},
4466 {3,3,6,8, 0},
4467 {4,0,6,8,2, 0},
4468 {3,3,6,8, 3,0,1,4, 0},
4469 {5,1,4,6,8,2, 0},
4470 {3,3,6,8, 3,1,2,5, 0},
4471 {5,0,6,8,5,1, 0},
4472 {3,3,6,8, 4,0,2,5,4, 0},
4473 {4,4,6,8,5, 0},
4474 {3,4,7,6, 0},
4475 {3,4,7,6, 3,0,3,2, 0},
4476 {4,0,1,7,6, 0},
4477 {5,1,7,6,3,2, 0},
4478 {3,4,7,6, 3,1,2,5, 0},
4479 {3,4,7,6, 4,0,3,5,1, 0},
4480 {5,0,2,5,7,6, 0},
4481 {4,3,5,7,6, 0},
4482 {4,3,4,7,8, 0},
4483 {5,0,4,7,8,2, 0},
4484 {5,0,1,7,8,3, 0},
4485 {4,1,7,8,2, 0},
4486 {4,3,4,7,8, 3,1,2,5, 0},
4487 {3,0,4,1, 3,5,7,8, 0},
4488 {3,0,2,3, 3,5,7,8, 0},
4489 {3,5,7,8, 0}
4490 };
4491 static int wdgEdges[WDG_EDGES][2] = {
4492 {0,1},
4493 {1,2},
4494 {2,0},
4495 {0,3},
4496 {1,4},
4497 {2,5},
4498 {3,4},
4499 {4,5},
4500 {5,3}
4501 };
4502
4503 /*----- hex elements -----*/
4504
4505 #define HEX_SIZE 127
4506 #define HEX_EDGES 12
4507
4508 static int hexCuts[HEX_SIZE+1][17] = {
4509 {0},
4510 {3,3,0,4, 0},
4511 {3,0,1,5, 0},
4512 {4,3,1,5,4, 0},
4513 {3,1,2,6, 0},
4514 {3,3,0,4, 3,1,2,6, 0},
4515 {4,0,2,6,5, 0},
4516 {5,3,2,6,5,4, 0},
4517 {3,3,7,2, 0},
4518 {4,0,4,7,2, 0},
4519 {3,0,1,5, 3,3,7,2, 0},
4520 {5,4,7,2,1,5, 0},
4521 {4,3,7,6,1, 0},
4522 {5,0,4,7,6,1, 0},
4523 {5,0,3,7,6,5, 0},
4524 {4,7,6,5,4, 0},
4525 {3,4,11,8, 0},
4526 {4,0,8,11,3, 0},
4527 {3,0,1,5, 3,4,11,8, 0},
4528 {5,3,11,8,5,1, 0},
4529 {3,1,2,6, 3,4,11,8, 0},
4530 {4,0,8,11,3, 3,1,2,6, 0},
4531 {4,0,2,6,5, 3,4,11,8, 0},
4532 {6,3,2,6,5,8,11, 0},
4533 {3,3,7,2, 3,4,11,8, 0},
4534 {5,0,2,7,11,8, 0},
4535 {3,0,1,5, 3,3,7,2, 3,4,11,8, 0},
4536 {6,2,1,5,8,11,7, 0},
4537 {4,3,7,6,1, 3,4,11,8, 0},
4538 {6,0,8,11,7,6,1, 0},
4539 {5,0,3,7,6,5, 3,4,11,8, 0},
4540 {5,11,7,6,5,8, 0},
4541 {3,8,9,5, 0},
4542 {3,8,9,5, 3,3,0,4, 0},
4543 {4,0,8,9,1, 0},
4544 {5,4,3,1,9,8, 0},
4545 {3,8,9,5, 3,1,2,6, 0},
4546 {3,1,9,5, 3,3,0,4, 3,1,2,6, 0},
4547 {5,0,2,6,9,8, 0},
4548 {6,3,2,6,9,8,4, 0},
4549 {3,8,9,5, 3,3,7,2, 0},
4550 {4,0,4,7,2, 3,8,9,5, 0},
4551 {4,0,8,9,1, 3,3,7,2, 0},
4552 {6,4,7,2,1,9,8, 0},
4553 {4,3,7,6,1, 3,8,9,5, 0},
4554 {5,7,6,1,0,4, 3,8,9,5, 0},
4555 {6,0,3,7,6,9,8, 0},
4556 {5,4,7,6,9,8, 0},
4557 {4,4,11,9,5, 0},
4558 {5,0,3,11,9,5, 0},
4559 {5,0,4,11,9,1, 0},
4560 {4,3,1,9,11, 0},
4561 {4,4,11,9,5, 3,1,2,6, 0},
4562 {6,0,3,11,9,6,2, 0},
4563 {6,0,2,6,9,11,4, 0},
4564 {5,3,2,6,9,11, 0},
4565 {4,4,11,9,5, 3,3,7,2, 0},
4566 {6,0,2,7,11,9,5, 0},
4567 {5,11,9,1,0,4, 3,3,7,2, 0},
4568 {5,11,7,2,1,9, 0},
4569 {4,3,7,6,1, 4,4,11,9,5, 0},
4570 {4,11,7,6,9, 3,0,1,5, 0},
4571 {4,11,7,6,9, 3,3,0,4, 0},
4572 {4,11,7,6,9, 0},
4573 {3,9,10,6, 0},
4574 {3,9,10,6, 3,3,0,4, 0},
4575 {3,9,10,6, 3,0,1,5, 0},
4576 {4,4,3,1,5, 3,9,10,6, 0},
4577 {4,2,1,9,10, 0},
4578 {4,2,1,9,10, 3,3,0,4, 0},
4579 {5,0,2,10,9,5, 0},
4580 {6,3,2,10,9,5,4, 0},
4581 {3,3,7,2, 3,9,10,6, 0},
4582 {4,7,2,0,4, 3,9,10,6, 0},
4583 {3,0,1,5, 3,2,3,7, 3,9,10,6, 0},
4584 {4,4,7,10,9, 4,1,2,6,5, 0},
4585 {5,3,7,10,9,1, 0},
4586 {6,4,7,10,9,1,0, 0},
4587 {6,3,7,10,9,5,0, 0},
4588 {5,4,7,10,9,5, 0},
4589 {3,4,11,8, 3,9,10,6, 0},
4590 {4,0,8,11,3, 3,9,10,6, 0},
4591 {3,0,1,5, 3,4,11,8, 3,9,10,6, 0},
4592 {5,3,11,10,6,1, 3,1,9,5, 0},
4593 {4,1,2,10,9, 3,4,11,8, 0},
4594 {4,3,11,10,2, 4,0,8,9,1, 0},
4595 {5,4,11,10,2,0, 3,8,9,5, 0},
4596 {4,3,11,10,2, 3,8,9,5, 0},
4597 {3,3,7,2, 3,9,10,6, 3,4,11,8, 0},
4598 {5,0,2,6,9,8, 3,11,7,10, 0},
4599 {3,0,4,3, 3,11,7,10, 3,1,2,6, 3,8,9,5, 0},
4600 {3,1,9,5, 3,11,7,10, 3,1,2,6, 0},
4601 {5,3,1,9,8,4, 3,11,7,10, 0},
4602 {4,0,8,9,1, 3,11,7,10, 0},
4603 {3,11,7,10, 3,8,9,5, 3,3,0,4, 0},
4604 {3,11,7,10, 3,8,9,5, 0},
4605 {4,8,10,6,5, 0},
4606 {4,8,10,6,5, 3,3,0,4, 0},
4607 {5,0,8,10,6,1, 0},
4608 {6,3,4,8,10,6,1, 0},
4609 {5,10,2,1,5,8, 0},
4610 {5,3,2,10,8,4, 3,0,1,5, 0},
4611 {4,0,8,10,2, 0},
4612 {5,3,2,10,8,4, 0},
4613 {4,8,10,6,5, 3,3,7,2, 0},
4614 {4,4,7,10,8, 4,0,2,6,5, 0},
4615 {5,3,7,10,8,0, 3,1,2,6, 0},
4616 {4,4,7,10,8, 3,1,2,6, 0},
4617 {6,3,7,10,8,5,1, 0},
4618 {4,4,7,10,8, 3,0,1,5, 0},
4619 {5,3,7,10,8,0, 0},
4620 {4,4,7,10,8, 0},
4621 {5,4,11,10,6,5, 0},
4622 {6,0,3,11,10,6,5, 0},
4623 {6,4,11,10,6,1,0, 0},
4624 {5,3,11,10,6,1, 0},
4625 {6,4,11,10,2,1,5, 0},
4626 {4,3,11,10,2, 3,0,1,5, 0},
4627 {5,4,11,10,2,0, 0},
4628 {4,11,10,2,3, 0},
4629 {5,3,2,6,5,4, 3,11,7,10, 0},
4630 {4,0,2,6,5, 3,11,7,10, 0},
4631 {3,3,0,4, 3,1,2,6, 3,11,7,10, 0},
4632 {3,1,2,6, 3,11,7,10, 0},
4633 {4,4,3,1,5, 3,11,7,10, 0},
4634 {3,0,1,5, 3,11,7,10, 0},
4635 {3,11,7,10, 3,3,0,4, 0},
4636 {3,11,7,10, 0}
4637 };
4638 static int hexEdges[HEX_EDGES][2] = {
4639 {0,1},
4640 {1,2},
4641 {2,3},
4642 {3,0},
4643 {0,4},
4644 {1,5},
4645 {2,6},
4646 {3,7},
4647 {4,5},
4648 {5,6},
4649 {6,7},
4650 {7,4}
4651 };
4652
4653 static int n_cut_nodes;
4654 static HASH *cut_hash;
4655
4656 /*------------------------------------------------------------------*/
4657
intersect_polygon(int zonenum,int nnodes,cgsize_t * nodeid,HASH * edgehash)4658 static void intersect_polygon (int zonenum, int nnodes, cgsize_t *nodeid,
4659 HASH *edgehash)
4660 {
4661 int i, n, nn, i1, i2;
4662 cgsize_t id = -1;
4663 float *node;
4664 double s1, s2;
4665 CutNode cnode, *cn;
4666 Edge edge, *ep;
4667 Zone *z = &zones[zonenum];
4668 static char *funcname = "intersect_polygon";
4669
4670 if (nnodes < 3) return;
4671 node = z->nodes[nodeid[0]];
4672 for (s1 = 0.0, i = 0; i < 3; i++)
4673 s1 += (node[i] * cutplane.plane[i]);
4674 i1 = (s1 - cutplane.plane[3]) >= 0.0 ? 1 : -1;
4675
4676 for (n = 1; n <= nnodes; n++) {
4677 nn = n % nnodes;
4678 node = z->nodes[nodeid[nn]];
4679 for (s2 = 0.0, i = 0; i < 3; i++)
4680 s2 += (node[i] * cutplane.plane[i]);
4681 i2 = (s2 - cutplane.plane[3]) >= 0.0 ? 1 : -1;
4682 if (i1 * i2 < 0) {
4683 cnode.nodes[0] = zonenum;
4684 if (nodeid[n-1] < nodeid[nn]) {
4685 cnode.nodes[1] = nodeid[n-1];
4686 cnode.nodes[2] = nodeid[nn];
4687 if (s1 == s2)
4688 cnode.ratio = 0.0;
4689 else
4690 cnode.ratio = (float)((cutplane.plane[3] - s1) / (s2 - s1));
4691 }
4692 else {
4693 cnode.nodes[1] = nodeid[nn];
4694 cnode.nodes[2] = nodeid[n-1];
4695 if (s1 == s2)
4696 cnode.ratio = 0.0;
4697 else
4698 cnode.ratio = (float)((cutplane.plane[3] - s2) / (s1 - s2));
4699 }
4700 cn = (CutNode *) HashFind (cut_hash, &cnode);
4701 if (cn == NULL) {
4702 cn = (CutNode *) MALLOC (funcname, sizeof(CutNode));
4703 for (i = 0; i < 3; i++)
4704 cn->nodes[i] = cnode.nodes[i];
4705 cn->id = n_cut_nodes++;
4706 cn->ratio = cnode.ratio;
4707 (void) HashAdd (cut_hash, cn);
4708 }
4709 if (id >= 0) {
4710 edge.nodes[0] = id;
4711 edge.nodes[1] = cn->id;
4712 ep = (Edge *) HashFind (edgehash, &edge);
4713 if (NULL == ep) {
4714 ep = (Edge *) MALLOC (funcname, sizeof(Edge));
4715 ep->nodes[0] = edge.nodes[0];
4716 ep->nodes[1] = edge.nodes[1];
4717 (void) HashAdd (edgehash, ep);
4718 }
4719 }
4720 id = cn->id;
4721 }
4722 s1 = s2;
4723 i1 = i2;
4724 }
4725 }
4726
4727 /*------------------------------------------------------------------*/
4728
intersect_element(int zonenum,CGNS_ENUMT (ElementType_t)elemtype,cgsize_t * nodeid,HASH * edgehash)4729 static void intersect_element (int zonenum, CGNS_ENUMT(ElementType_t) elemtype,
4730 cgsize_t *nodeid, HASH *edgehash)
4731 {
4732 int i, n, index, nn, nc, *cuts;
4733 int edgemap[HEX_EDGES][2], ids[7];
4734 cgsize_t swap;
4735 double s1, s2;
4736 float *node;
4737 CutNode cnode, *cn;
4738 Edge edge, *ep;
4739 Zone *z = &zones[zonenum];
4740 static char *funcname = "intersect_element";
4741
4742 /* get intersection lookup table */
4743
4744 switch (elemtype) {
4745 case CGNS_ENUMV(TRI_3):
4746 case CGNS_ENUMV(TRI_6):
4747 case CGNS_ENUMV(TRI_9):
4748 case CGNS_ENUMV(TRI_10):
4749 case CGNS_ENUMV(TRI_12):
4750 case CGNS_ENUMV(TRI_15):
4751 index = classify_element(z, 3, nodeid);
4752 if (index < 1 || index > TRI_SIZE) return;
4753 cuts = triCuts[index];
4754 for (n = 0; n < TRI_EDGES; n++) {
4755 edgemap[n][0] = triEdges[n][0];
4756 edgemap[n][1] = triEdges[n][1];
4757 }
4758 break;
4759 case CGNS_ENUMV(QUAD_4):
4760 case CGNS_ENUMV(QUAD_8):
4761 case CGNS_ENUMV(QUAD_9):
4762 case CGNS_ENUMV(QUAD_12):
4763 case CGNS_ENUMV(QUAD_16):
4764 case CGNS_ENUMV(QUAD_P4_16):
4765 case CGNS_ENUMV(QUAD_25):
4766 index = classify_element(z, 4, nodeid);
4767 if (index < 1 || index > QUAD_SIZE) return;
4768 cuts = quadCuts[index];
4769 for (n = 0; n < QUAD_EDGES; n++) {
4770 edgemap[n][0] = quadEdges[n][0];
4771 edgemap[n][1] = quadEdges[n][1];
4772 }
4773 break;
4774 case CGNS_ENUMV(TETRA_4):
4775 case CGNS_ENUMV(TETRA_10):
4776 case CGNS_ENUMV(TETRA_16):
4777 case CGNS_ENUMV(TETRA_20):
4778 case CGNS_ENUMV(TETRA_22):
4779 case CGNS_ENUMV(TETRA_34):
4780 case CGNS_ENUMV(TETRA_35):
4781 index = classify_element(z, 4, nodeid);
4782 if (index < 1 || index > TET_SIZE) return;
4783 cuts = tetCuts[index];
4784 for (n = 0; n < TET_EDGES; n++) {
4785 edgemap[n][0] = tetEdges[n][0];
4786 edgemap[n][1] = tetEdges[n][1];
4787 }
4788 break;
4789 case CGNS_ENUMV(PYRA_5):
4790 case CGNS_ENUMV(PYRA_13):
4791 case CGNS_ENUMV(PYRA_14):
4792 case CGNS_ENUMV(PYRA_21):
4793 case CGNS_ENUMV(PYRA_29):
4794 case CGNS_ENUMV(PYRA_30):
4795 case CGNS_ENUMV(PYRA_P4_29):
4796 case CGNS_ENUMV(PYRA_50):
4797 case CGNS_ENUMV(PYRA_55):
4798 index = classify_element(z, 5, nodeid);
4799 if (index < 1 || index > PYR_SIZE) return;
4800 cuts = pyrCuts[index];
4801 for (n = 0; n < PYR_EDGES; n++) {
4802 edgemap[n][0] = pyrEdges[n][0];
4803 edgemap[n][1] = pyrEdges[n][1];
4804 }
4805 break;
4806 case CGNS_ENUMV(PENTA_6):
4807 case CGNS_ENUMV(PENTA_15):
4808 case CGNS_ENUMV(PENTA_18):
4809 case CGNS_ENUMV(PENTA_24):
4810 case CGNS_ENUMV(PENTA_38):
4811 case CGNS_ENUMV(PENTA_40):
4812 case CGNS_ENUMV(PENTA_33):
4813 case CGNS_ENUMV(PENTA_66):
4814 case CGNS_ENUMV(PENTA_75):
4815 index = classify_element(z, 6, nodeid);
4816 if (index < 1 || index > WDG_SIZE) return;
4817 cuts = wdgCuts[index];
4818 for (n = 0; n < WDG_EDGES; n++) {
4819 edgemap[n][0] = wdgEdges[n][0];
4820 edgemap[n][1] = wdgEdges[n][1];
4821 }
4822 break;
4823 case CGNS_ENUMV(HEXA_8):
4824 case CGNS_ENUMV(HEXA_20):
4825 case CGNS_ENUMV(HEXA_27):
4826 case CGNS_ENUMV(HEXA_32):
4827 case CGNS_ENUMV(HEXA_56):
4828 case CGNS_ENUMV(HEXA_64):
4829 case CGNS_ENUMV(HEXA_44):
4830 case CGNS_ENUMV(HEXA_98):
4831 case CGNS_ENUMV(HEXA_125):
4832 index = classify_element(z, 8, nodeid);
4833 if (index < 1 || index > HEX_SIZE) return;
4834 cuts = hexCuts[index];
4835 for (n = 0; n < HEX_EDGES; n++) {
4836 edgemap[n][0] = hexEdges[n][0];
4837 edgemap[n][1] = hexEdges[n][1];
4838 }
4839 break;
4840 default:
4841 return;
4842 }
4843
4844 /* get the edge intersections */
4845
4846 for (nc = 0; cuts[nc];) {
4847 nn = cuts[nc];
4848 for (n = 1; n <= nn; n++) {
4849
4850 /* get edge nodes */
4851
4852 cnode.nodes[0] = zonenum;
4853 cnode.nodes[1] = nodeid[edgemap[cuts[nc+n]][0]];
4854 cnode.nodes[2] = nodeid[edgemap[cuts[nc+n]][1]];
4855 if (cnode.nodes[2] < cnode.nodes[1]) {
4856 swap = cnode.nodes[1];
4857 cnode.nodes[1] = cnode.nodes[2];
4858 cnode.nodes[2] = swap;
4859 }
4860 cn = (CutNode *) HashFind (cut_hash, &cnode);
4861
4862 /* add node to hash table if not there */
4863
4864 if (NULL == cn) {
4865 cn = (CutNode *) MALLOC (funcname, sizeof(CutNode));
4866 cn->id = n_cut_nodes++;
4867 for (i = 0; i < 3; i++)
4868 cn->nodes[i] = cnode.nodes[i];
4869 node = z->nodes[cn->nodes[1]];
4870 for (s1 = 0.0, i = 0; i < 3; i++)
4871 s1 += (node[i] * cutplane.plane[i]);
4872 node = z->nodes[cn->nodes[2]];
4873 for (s2 = 0.0, i = 0; i < 3; i++)
4874 s2 += (node[i] * cutplane.plane[i]);
4875 if (s1 == s2)
4876 cn->ratio = 0.0;
4877 else
4878 cn->ratio = (float)((cutplane.plane[3] - s1) / (s2 - s1));
4879 (void) HashAdd (cut_hash, cn);
4880 }
4881 ids[n-1] = cn->id;
4882 }
4883 ids[nn] = ids[0];
4884
4885 /* add cutplane edge */
4886
4887 for (n = 0; n < nn; n++) {
4888 edge.nodes[0] = ids[n];
4889 edge.nodes[1] = ids[n+1];
4890 ep = (Edge *) HashFind (edgehash, &edge);
4891 if (NULL == ep) {
4892 ep = (Edge *) MALLOC (funcname, sizeof(Edge));
4893 ep->nodes[0] = edge.nodes[0];
4894 ep->nodes[1] = edge.nodes[1];
4895 (void) HashAdd (edgehash, ep);
4896 }
4897 }
4898
4899 /* next cut */
4900
4901 nc += (nn + 1);
4902 }
4903 }
4904
4905 /*------------------------------------------------------------------*/
4906
find_intersects()4907 static cgsize_t find_intersects ()
4908 {
4909 int nz, nr, nf, nfaces, nnodes;
4910 cgsize_t n, ne;
4911 size_t nn;
4912 CGNS_ENUMT(ElementType_t) type;
4913 Face *f;
4914 Regn *r;
4915 HASH *edgehash;
4916 static char *funcname = "find_intersects";
4917
4918 /* create hash table to store nodes at edge intersections */
4919
4920 n_cut_nodes = 0;
4921 nn = (size_t)cutplane.nelems;
4922 cut_hash = HashCreate (nn > 1024 ? nn / 3 : 127,
4923 compare_cut_node, hash_cut_node);
4924 cutplane.nedges = 0;
4925
4926 for (nz = 0; nz < nzones; nz++) {
4927 for (nr = 0; nr < zones[nz].nregs; nr++) {
4928 r = &zones[nz].regs[nr];
4929 if (r->cut.nelems == 0) continue;
4930 type = r->elemtype;
4931
4932 nn = (size_t)r->cut.nelems;
4933 edgehash = HashCreate (nn > 1024 ? nn / 3: 127,
4934 compare_edges, hash_edge);
4935
4936 for (n = 0, ne = 0; ne < r->cut.nelems; ne++) {
4937 n = r->cut.elems[ne];
4938 if (r->elemtype == CGNS_ENUMV(NGON_n)) {
4939 nnodes = r->elem_offsets[n+1]-r->elem_offsets[n];
4940 intersect_polygon(nz, nnodes, &r->elems[r->elem_offsets[n]], edgehash);
4941 }
4942 else if (r->elemtype == CGNS_ENUMV(NFACE_n)) {
4943 nfaces = r->elem_offsets[n+1]-r->elem_offsets[n];
4944 for (nf = 0; nf < nfaces; nf++) {
4945 f = r->poly[abs(r->elems[r->elem_offsets[n]+nf])-1];
4946 intersect_polygon(nz, f->nnodes, f->nodes, edgehash);
4947 }
4948 }
4949 else {
4950 cg_npe(type, &nnodes);
4951 if (r->elemtype == CGNS_ENUMV(MIXED)) {
4952 type = (CGNS_ENUMT(ElementType_t))r->elems[n++];
4953 cg_npe(type, &nnodes);
4954 }
4955 intersect_element (nz, type, &r->elems[n], edgehash);
4956 }
4957 }
4958
4959 r->cut.nedges = 0;
4960 nn = HashSize (edgehash);
4961 if (nn) {
4962 r->cut.edges = (Edge *) MALLOC (funcname, nn * sizeof(Edge));
4963 HashList (edgehash, get_cut_edge, &r->cut);
4964 }
4965 HashDestroy (edgehash, NULL);
4966 cutplane.nedges += r->cut.nedges;
4967 }
4968 }
4969
4970 nn = HashSize (cut_hash);
4971 cutplane.nnodes = (cgsize_t)nn;
4972 cutplane.nodes = (Node *) MALLOC (funcname, nn * sizeof(Node));
4973 HashDestroy (cut_hash, get_cut_node);
4974
4975 return cutplane.nelems;
4976 }
4977
4978 /*------------------------------------------------------------------*/
4979
draw_edges()4980 static void draw_edges ()
4981 {
4982 int nz, nr;
4983 cgsize_t ne, nn;
4984 Zone *z;
4985 Regn *r;
4986
4987 glDisable(GL_LIGHTING);
4988 glShadeModel(GL_FLAT);
4989 glBegin(GL_LINES);
4990 glColor3fv(cutcolor);
4991
4992 for (nz = 0; nz < nzones; nz++) {
4993 z = &zones[nz];
4994 for (nr = 0; nr < z->nregs; nr++) {
4995 r = &z->regs[nr];
4996 if (r->cut.nedges == 0) continue;
4997 if (!usecutclr) glColor3fv(r->color);
4998 for (ne = 0; ne < r->cut.nedges; ne++) {
4999 nn = r->cut.edges[ne].nodes[0];
5000 glVertex3fv(cutplane.nodes[nn]);
5001 nn = r->cut.edges[ne].nodes[1];
5002 glVertex3fv(cutplane.nodes[nn]);
5003 }
5004 }
5005 }
5006
5007 glEnd();
5008 }
5009
5010 /*------------------------------------------------------------------*/
5011
draw_elements(int mode)5012 static void draw_elements (int mode)
5013 {
5014 int nz, nr, i, j;
5015 int ip, nf, nnodes;
5016 cgsize_t n, ne, nn;
5017 float *nodes[4], norm[3];
5018 CGNS_ENUMT(ElementType_t) type;
5019 Zone *z;
5020 Regn *r;
5021 Face *f;
5022
5023 glEnable(GL_LIGHTING);
5024 glShadeModel(GL_FLAT);
5025 glPolygonMode(GL_FRONT_AND_BACK, mode);
5026 glColor3fv(cutcolor);
5027 glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, cutcolor);
5028
5029 for (nz = 0; nz < nzones; nz++) {
5030 z = &zones[nz];
5031 for (nr = 0; nr < z->nregs; nr++) {
5032 r = &z->regs[nr];
5033 if (r->cut.nelems == 0) continue;
5034 type = r->elemtype;
5035
5036 if (!usecutclr) {
5037 glColor3fv(r->color);
5038 glMaterialfv(GL_FRONT_AND_BACK,
5039 GL_AMBIENT_AND_DIFFUSE, r->color);
5040 }
5041
5042 if (type == CGNS_ENUMV(NGON_n)) {
5043 for (ne = 0; ne < r->cut.nelems; ne++) {
5044 nn = r->cut.elems[ne];
5045 nnodes = r->elem_offsets[nn+1]- r->elem_offsets[nn];
5046 if (nnodes < 3) continue;
5047 if (nnodes == 3)
5048 glBegin(GL_TRIANGLES);
5049 else if (nnodes == 4)
5050 glBegin(GL_QUADS);
5051 else
5052 glBegin(GL_POLYGON);
5053 glNormal3fv(face_normal(z, nnodes, &r->elems[r->elem_offsets[nn]]));
5054 for (i = r->elem_offsets[nn]; i < r->elem_offsets[nn+1]; i++)
5055 glVertex3fv(z->nodes[r->elems[i]]);
5056 glEnd();
5057 }
5058 }
5059 else if (type == CGNS_ENUMV(NFACE_n)) {
5060 for (ne = 0; ne < r->cut.nelems; ne++) {
5061 nn = r->cut.elems[ne];
5062 nf = r->elem_offsets[nn+1]- r->elem_offsets[nn];
5063 for (j = 0; j < nf; j++) {
5064 f = r->poly[abs(r->elems[r->elem_offsets[nn]+j])-1];
5065 if (f->nnodes < 3) continue;
5066 if (f->nnodes == 3)
5067 glBegin(GL_TRIANGLES);
5068 else if (f->nnodes == 4)
5069 glBegin(GL_QUADS);
5070 else
5071 glBegin(GL_POLYGON);
5072 if (r->elems[r->elem_offsets[nn]+j] < 0) {
5073 for (i = 0; i < 3; i++)
5074 norm[i] = -f->normal[i];
5075 }
5076 else {
5077 for (i = 0; i < 3; i++)
5078 norm[i] = f->normal[i];
5079 }
5080 glNormal3fv(norm);
5081 for (i = 0; i < f->nnodes; i++)
5082 glVertex3fv(z->nodes[f->nodes[i]]);
5083 glEnd();
5084 }
5085 }
5086 }
5087 else {
5088 for (ne = 0; ne < r->cut.nelems; ne++) {
5089 n = r->cut.elems[ne];
5090 if (r->elemtype == CGNS_ENUMV(MIXED))
5091 type = (CGNS_ENUMT(ElementType_t))r->elems[n++];
5092 switch (type) {
5093 case CGNS_ENUMV(TRI_3):
5094 case CGNS_ENUMV(TRI_6):
5095 case CGNS_ENUMV(TRI_9):
5096 case CGNS_ENUMV(TRI_10):
5097 case CGNS_ENUMV(TRI_12):
5098 case CGNS_ENUMV(TRI_15):
5099 ip = 0;
5100 nf = 1;
5101 break;
5102 case CGNS_ENUMV(QUAD_4):
5103 case CGNS_ENUMV(QUAD_8):
5104 case CGNS_ENUMV(QUAD_9):
5105 case CGNS_ENUMV(QUAD_12):
5106 case CGNS_ENUMV(QUAD_16):
5107 case CGNS_ENUMV(QUAD_P4_16):
5108 case CGNS_ENUMV(QUAD_25):
5109 ip = 1;
5110 nf = 1;
5111 break;
5112 case CGNS_ENUMV(TETRA_4):
5113 case CGNS_ENUMV(TETRA_10):
5114 case CGNS_ENUMV(TETRA_16):
5115 case CGNS_ENUMV(TETRA_20):
5116 case CGNS_ENUMV(TETRA_22):
5117 case CGNS_ENUMV(TETRA_34):
5118 case CGNS_ENUMV(TETRA_35):
5119 ip = 2;
5120 nf = 4;
5121 break;
5122 case CGNS_ENUMV(PYRA_5):
5123 case CGNS_ENUMV(PYRA_13):
5124 case CGNS_ENUMV(PYRA_14):
5125 case CGNS_ENUMV(PYRA_21):
5126 case CGNS_ENUMV(PYRA_29):
5127 case CGNS_ENUMV(PYRA_30):
5128 case CGNS_ENUMV(PYRA_P4_29):
5129 case CGNS_ENUMV(PYRA_50):
5130 case CGNS_ENUMV(PYRA_55):
5131 ip = 6;
5132 nf = 5;
5133 break;
5134 case CGNS_ENUMV(PENTA_6):
5135 case CGNS_ENUMV(PENTA_15):
5136 case CGNS_ENUMV(PENTA_18):
5137 case CGNS_ENUMV(PENTA_24):
5138 case CGNS_ENUMV(PENTA_38):
5139 case CGNS_ENUMV(PENTA_40):
5140 case CGNS_ENUMV(PENTA_33):
5141 case CGNS_ENUMV(PENTA_66):
5142 case CGNS_ENUMV(PENTA_75):
5143 ip = 11;
5144 nf = 5;
5145 break;
5146 case CGNS_ENUMV(HEXA_8):
5147 case CGNS_ENUMV(HEXA_20):
5148 case CGNS_ENUMV(HEXA_27):
5149 case CGNS_ENUMV(HEXA_32):
5150 case CGNS_ENUMV(HEXA_56):
5151 case CGNS_ENUMV(HEXA_64):
5152 case CGNS_ENUMV(HEXA_44):
5153 case CGNS_ENUMV(HEXA_98):
5154 case CGNS_ENUMV(HEXA_125):
5155 ip = 16;
5156 nf = 6;
5157 break;
5158 default:
5159 ip = 0;
5160 nf = 0;
5161 break;
5162 }
5163 for (j = 0; j < nf; j++) {
5164 nnodes = facenodes[ip+j][0];
5165 for (i = 0; i < nnodes; i++) {
5166 nn = r->elems[n+facenodes[ip+j][i+1]];
5167 nodes[i] = z->nodes[nn];
5168 }
5169 if (nnodes == 4) {
5170 glBegin(GL_QUADS);
5171 glNormal3fv(compute_normal(nodes[0], nodes[1],
5172 nodes[2], nodes[3]));
5173 }
5174 else {
5175 glBegin(GL_TRIANGLES);
5176 glNormal3fv(compute_normal(nodes[0], nodes[1],
5177 nodes[2], NULL));
5178 }
5179 for (i = 0; i < nnodes; i++)
5180 glVertex3fv(nodes[i]);
5181 glEnd();
5182 }
5183 }
5184 }
5185 }
5186 }
5187 }
5188
5189 /*---------- OGLcutplane -------------------------------------------
5190 * create OGL display list for a cutting plane
5191 *------------------------------------------------------------------*/
5192
OGLcutplane(ClientData data,Tcl_Interp * interp,int argc,char ** argv)5193 static int OGLcutplane (ClientData data, Tcl_Interp *interp, int argc, char **argv)
5194 {
5195 int mode;
5196 float plane[4];
5197 static char slist[33];
5198
5199 if (argc < 1 || argc > 3) {
5200 Tcl_SetResult (interp, "usage: OGLcutplane [mode] [plane]",
5201 TCL_STATIC);
5202 return TCL_ERROR;
5203 }
5204
5205 /* create and return displaylist flag */
5206
5207 if (argc == 1) {
5208 if (!CutDL) CutDL = glGenLists (1);
5209 sprintf (slist, "%d", CutDL);
5210 Tcl_SetResult (interp, slist, TCL_STATIC);
5211 return TCL_OK;
5212 }
5213
5214 mode = atoi(argv[1]);
5215
5216 if (argc == 3) {
5217 int np;
5218 CONST char **args;
5219 if (TCL_OK != Tcl_SplitList (interp, argv[2], &np, &args))
5220 return TCL_ERROR;
5221 if (np != 4) {
5222 Tcl_Free ((char *)args);
5223 Tcl_SetResult (interp, "invalid plane", TCL_STATIC);
5224 return TCL_ERROR;
5225 }
5226 for (np = 0; np < 4; np++)
5227 plane[np] = (float) atof (args[np]);
5228 Tcl_Free ((char *)args);
5229 init_cutplane(plane);
5230 find_elements();
5231 }
5232
5233 if (!CutDL) CutDL = glGenLists (1);
5234
5235 glNewList (CutDL, GL_COMPILE);
5236 if (mode && cutplane.nelems) {
5237 if (mode == 1) {
5238 if (cutplane.nedges == 0)
5239 find_intersects();
5240 draw_edges ();
5241 }
5242 else {
5243 draw_elements (mode > 2 ? GL_FILL : GL_LINE);
5244 }
5245 }
5246 glEndList ();
5247
5248 sprintf (slist, "%ld %ld", (long)cutplane.nelems, (long)cutplane.nedges);
5249 Tcl_SetResult (interp, slist, TCL_STATIC);
5250 return TCL_OK;
5251 }
5252
5253 /*---------- OGLdrawplane ------------------------------------------
5254 * draw the cutting plane
5255 *------------------------------------------------------------------*/
5256
OGLdrawplane(ClientData data,Tcl_Interp * interp,int argc,char ** argv)5257 static int OGLdrawplane (ClientData data, Tcl_Interp *interp, int argc, char **argv)
5258 {
5259 int n, np, i, j, k, index, n0, n1;
5260 CONST char **args;
5261 float plane[4], bbox[3][2], s[8], ds;
5262 float node[8][3], pnode[6][3];
5263 static char slist[17];
5264
5265 if (argc < 1 || argc > 2) {
5266 Tcl_SetResult (interp, "usage: OGLdrawplane [plane]",
5267 TCL_STATIC);
5268 return TCL_ERROR;
5269 }
5270
5271 if (!PlaneDL) PlaneDL = glGenLists (1);
5272
5273 if (argc == 1) {
5274 glNewList (PlaneDL, GL_COMPILE);
5275 glEndList ();
5276 sprintf (slist, "%d", PlaneDL);
5277 Tcl_SetResult (interp, slist, TCL_STATIC);
5278 return TCL_OK;
5279 }
5280
5281 if (TCL_OK != Tcl_SplitList (interp, argv[1], &np, &args))
5282 return TCL_ERROR;
5283 if (np != 4) {
5284 Tcl_Free ((char *)args);
5285 Tcl_SetResult (interp, "invalid plane", TCL_STATIC);
5286 return TCL_ERROR;
5287 }
5288 for (n = 0; n < np; n++)
5289 plane[n] = (float) atof (args[n]);
5290 Tcl_Free ((char *)args);
5291
5292 get_bounds (ignorevis, bbox);
5293 index = n = 0;
5294 for (k = 0; k < 2; k++) {
5295 for (j = 0; j < 2; j++) {
5296 for (i = 0; i < 2; i++) {
5297 node[n][0] = bbox[0][(i+j)%2];
5298 node[n][1] = bbox[1][j];
5299 node[n][2] = bbox[2][k];
5300 s[n] = node[n][0] * plane[0] + node[n][1] * plane[1] +
5301 node[n][2] * plane[2];
5302 if (s[n] >= plane[3]) index |= (1 << n);
5303 n++;
5304 }
5305 }
5306 }
5307 if (index > 0x7f) index ^= 0xff;
5308 if (index < 1 || index > HEX_SIZE) {
5309 Tcl_SetResult (interp, "plane doesn't intersect", TCL_STATIC);
5310 return TCL_ERROR;
5311 }
5312
5313 np = hexCuts[index][0];
5314 for (n = 0; n < np; n++) {
5315 j = hexCuts[index][n+1];
5316 n0 = hexEdges[j][0];
5317 n1 = hexEdges[j][1];
5318 ds = s[n1] - s[n0];
5319 if (ds != 0.0)
5320 ds = (plane[3] - s[n0]) / ds;
5321 for (i = 0; i < 3; i++)
5322 pnode[n][i] = node[n0][i] + ds * (node[n1][i] - node[n0][i]);
5323 }
5324
5325 glNewList (PlaneDL, GL_COMPILE);
5326 glEnable(GL_BLEND);
5327 glDepthMask(GL_FALSE);
5328 glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
5329 glDisable(GL_LIGHTING);
5330 glShadeModel(GL_FLAT);
5331 glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
5332 glColor4fv(cutcolor);
5333 glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, cutcolor);
5334
5335 glBegin(GL_TRIANGLE_FAN);
5336 glNormal3fv(plane);
5337 for (n = 0; n < np; n++)
5338 glVertex3fv(pnode[n]);
5339 glEnd();
5340
5341 glDepthMask(GL_TRUE);
5342 glDisable(GL_BLEND);
5343 glEnable(GL_LIGHTING);
5344 glEndList ();
5345
5346 sprintf (slist, "%d", PlaneDL);
5347 Tcl_SetResult (interp, slist, TCL_STATIC);
5348 return TCL_OK;
5349 }
5350
5351 /*---------- OGLcutconfig ------------------------------------------
5352 * set the cutting plane color and operation
5353 *------------------------------------------------------------------*/
5354
OGLcutconfig(ClientData data,Tcl_Interp * interp,int argc,char ** argv)5355 static int OGLcutconfig (ClientData data, Tcl_Interp *interp, int argc, char **argv)
5356 {
5357 int n, np;
5358 CONST char **args;
5359
5360 if (argc < 2 || argc > 4) {
5361 Tcl_SetResult (interp, "usage: OGLcutconfig color [usecutclr] [ignorevis]",
5362 TCL_STATIC);
5363 return TCL_ERROR;
5364 }
5365
5366 if (TCL_OK != Tcl_SplitList (interp, argv[1], &np, &args))
5367 return TCL_ERROR;
5368 if (np < 3 || np > 4) {
5369 Tcl_Free ((char *)args);
5370 Tcl_SetResult (interp, "invalid color", TCL_STATIC);
5371 return TCL_ERROR;
5372 }
5373 for (n = 0; n < np; n++)
5374 cutcolor[n] = (float) atof (args[n]);
5375 Tcl_Free ((char *)args);
5376
5377 if (argc > 2) {
5378 usecutclr = atoi (argv[2]);
5379 if (argc > 3)
5380 ignorevis = atoi (argv[3]);
5381 }
5382 return TCL_OK;
5383 }
5384
5385 #endif
5386
5387 /*---------- Cgnstcl_Init --------------------------------------
5388 * Initialize and create the commands
5389 *--------------------------------------------------------------*/
5390
5391 #if defined(_WIN32) && defined(BUILD_DLL)
5392 __declspec(dllexport)
5393 #endif
Cgnstcl_Init(Tcl_Interp * interp)5394 int Cgnstcl_Init(Tcl_Interp *interp)
5395 {
5396 global_interp = interp;
5397 Tcl_CreateCommand (interp, "CGNSopen", (Tcl_CmdProc *)CGNSopen,
5398 (ClientData)0, (Tcl_CmdDeleteProc *)0);
5399 Tcl_CreateCommand (interp, "CGNSclose", (Tcl_CmdProc *)CGNSclose,
5400 (ClientData)0, (Tcl_CmdDeleteProc *)0);
5401 Tcl_CreateCommand (interp, "CGNSbase", (Tcl_CmdProc *)CGNSbase,
5402 (ClientData)0, (Tcl_CmdDeleteProc *)0);
5403 Tcl_CreateCommand (interp, "CGNSzone", (Tcl_CmdProc *)CGNSzone,
5404 (ClientData)0, (Tcl_CmdDeleteProc *)0);
5405 Tcl_CreateCommand (interp, "CGNSsummary", (Tcl_CmdProc *)CGNSsummary,
5406 (ClientData)0, (Tcl_CmdDeleteProc *)0);
5407 Tcl_CreateCommand (interp, "CGNSgetbase", (Tcl_CmdProc *)CGNSgetbase,
5408 (ClientData)0, (Tcl_CmdDeleteProc *)0);
5409 Tcl_CreateCommand (interp, "CGNSgetzone", (Tcl_CmdProc *)CGNSgetzone,
5410 (ClientData)0, (Tcl_CmdDeleteProc *)0);
5411 Tcl_CreateCommand (interp, "CGNSgetregion", (Tcl_CmdProc *)CGNSgetregion,
5412 (ClientData)0, (Tcl_CmdDeleteProc *)0);
5413 Tcl_CreateCommand (interp, "CGNSregiondim", (Tcl_CmdProc *)CGNSregiondim,
5414 (ClientData)0, (Tcl_CmdDeleteProc *)0);
5415 Tcl_CreateCommand (interp, "CGNSbounds", (Tcl_CmdProc *)CGNSbounds,
5416 (ClientData)0, (Tcl_CmdDeleteProc *)0);
5417 Tcl_CreateCommand (interp, "OGLregion", (Tcl_CmdProc *)OGLregion,
5418 (ClientData)0, (Tcl_CmdDeleteProc *)0);
5419 Tcl_CreateCommand (interp, "OGLaxis", (Tcl_CmdProc *)OGLaxis,
5420 (ClientData)0, (Tcl_CmdDeleteProc *)0);
5421 Tcl_CreateCommand (interp, "OGLcolor", (Tcl_CmdProc *)OGLcolor,
5422 (ClientData)0, (Tcl_CmdDeleteProc *)0);
5423 #ifndef NO_CUTTING_PLANE
5424 Tcl_CreateCommand (interp, "OGLcutplane", (Tcl_CmdProc *)OGLcutplane,
5425 (ClientData)0, (Tcl_CmdDeleteProc *)0);
5426 Tcl_CreateCommand (interp, "OGLdrawplane", (Tcl_CmdProc *)OGLdrawplane,
5427 (ClientData)0, (Tcl_CmdDeleteProc *)0);
5428 Tcl_CreateCommand (interp, "OGLcutconfig", (Tcl_CmdProc *)OGLcutconfig,
5429 (ClientData)0, (Tcl_CmdDeleteProc *)0);
5430 #endif
5431 return TCL_OK;
5432 }
5433
5434