1 /*
2 * import routines for CGNS
3 */
4
5 #include <stdio.h>
6 #include <stdlib.h>
7 #include <string.h>
8 #include <ctype.h>
9 #include <math.h>
10
11 #include "cgnsImport.h"
12 #include "hash.h"
13
14 #ifndef CG_MODE_READ
15 # define CG_MODE_READ MODE_READ
16 # define CG_MODE_WRITE MODE_WRITE
17 # define CG_MODE_MODIFY MODE_MODIFY
18 #endif
19
20 #ifndef DOUBLE /* data type for coordinates */
21 #define DOUBLE double
22 #endif
23
24 #ifndef TOLERANCE /* tolerance for duplicate checking */
25 #define TOLERANCE 1.e-03
26 #endif
27
28 #ifndef REGION_BASE /* base name for unnamed regions */
29 #define REGION_BASE "Region"
30 #endif
31
32 /*--- variables ---*/
33
34 static int num_vars = 0;
35 static char **var_names = 0;
36
37 /*--- node bit flags ---*/
38
39 #define USED_BIT 1
40 #define REGN_BIT 2
41 #define REFS_BIT 4
42
43 /*--- node structure ---*/
44
45 typedef struct {
46 cgsize_t id; /* node ID */
47 int flags; /* references to node */
48 DOUBLE x, y, z; /* coordinates */
49 DOUBLE dist; /* distance from origin */
50 DOUBLE *vars; /* variables */
51 } cgnsNODE;
52
53 /*--- node mapping data ---*/
54
55 #define NODEMAP_INC 50 /* realloc this many at a time */
56
57 typedef struct {
58 cgsize_t nodeid; /* node id */
59 cgsize_t mapped; /* set when mapped */
60 cgnsNODE *node; /* pointer to node data */
61 } NODEMAP;
62
63 static cgsize_t num_nodes = 0; /* number of nodes */
64 static cgsize_t max_nodes = 0; /* number of nodes malloced */
65 static NODEMAP *nodemap; /* list of nodes */
66
67 /*--- duplicate node checking ---*/
68
69 #define ENTRY_INC 50
70
71 static int no_check = 0;
72 static double def_tol = TOLERANCE;
73 static double tolerance = TOLERANCE;
74 DOUBLE xmin, xmax, ymin, ymax, zmin, zmax;
75
76 static cgsize_t num_entries = 0;
77 static cgsize_t max_entries = 0;
78 static cgnsNODE **node_list;
79
80 /*--- element data ---*/
81
82 #define ELEMENT_INC 50 /* realloc this many at a time */
83
84 typedef struct {
85 cgsize_t elemid; /* element ID */
86 int elemtype; /* element type (number of nodes) */
87 cgsize_t *nodeid; /* node ID's for element */
88 char facemap[6]; /* remapping of faces */
89 } cgnsELEM;
90
91 static cgsize_t num_elements = 0;/* number of elements */
92 static cgsize_t max_elements = 0;/* number of elements malloced */
93 static cgnsELEM *elemlist; /* list of elements */
94
95 /*--- region data ---*/
96
97 #define REGION_INC 50 /* step increment for realloc */
98
99 static char region_name[33]; /* region name */
100 static int region_type; /* type of region */
101 static cgsize_t region_id = 0; /* region ID */
102 static cgsize_t region_max = 0; /* malloced size of region_list */
103 static cgsize_t region_nobjs = 0; /* number objects in region */
104 static cgsize_t *region_list; /* list of nodes in region */
105
106 typedef struct {
107 char name[33]; /* region name */
108 int type; /* region type */
109 cgsize_t nobjs; /* number of objects */
110 cgsize_t *objid; /* object ID's */
111 } cgnsREGN;
112
113 static int num_regions = 0; /* number of regions */
114 static cgnsREGN *reglist; /* list of regions */
115
116 /*--- external faces */
117
118 typedef struct {
119 cgsize_t faceid;
120 int nnodes;
121 cgsize_t nodeid[4];
122 int flags;
123 } cgnsFACE;
124
125 static cgsize_t num_faces = 0;
126 static cgnsFACE **facelist;
127
128 /*--- CGNS data ---*/
129
130 int cgnsFile = 0;
131 int cgnsBase = 0;
132 int cgnsZone = 0;
133 char cgnsZoneName[33] = "";
134
135 /*--- error handling callback ---*/
136
137 static void (*errcallback)( /* callback pointer to user routine */
138 char *errmsg /* error message */
139 ) = NULL;
140
141 /*======================================================================
142 * Node routines
143 *======================================================================*/
144
145 /*---------- NewNode ---------------------------------------------------
146 * create a new node
147 *----------------------------------------------------------------------*/
148
NewNode(cgsize_t id,double x,double y,double z)149 static cgnsNODE *NewNode (cgsize_t id, double x, double y, double z)
150 {
151 cgnsNODE *node = (cgnsNODE *) malloc (sizeof(cgnsNODE));
152
153 if (NULL != node) {
154 node->id = id;
155 node->flags = 0;
156 node->x = x;
157 node->y = y;
158 node->z = z;
159 node->dist = 0.0;
160 node->vars = 0;
161 if (num_vars)
162 node->vars = (DOUBLE *) calloc (num_vars, sizeof(DOUBLE));
163 }
164 return (node);
165 }
166
167 /*---------- GetNode ---------------------------------------------------
168 * return the node for a given node id
169 *----------------------------------------------------------------------*/
170
GetNode(cgsize_t nodeid,cgsize_t * pos)171 static cgnsNODE *GetNode (cgsize_t nodeid, cgsize_t *pos)
172 {
173 cgsize_t lo = 0, hi = num_nodes - 1;
174
175 *pos = 0;
176 if (!num_nodes || nodeid < nodemap[0].nodeid)
177 return (NULL);
178 if (nodeid == nodemap[0].nodeid)
179 return (nodemap[0].node);
180 if (!hi || nodeid > nodemap[hi].nodeid) {
181 *pos = num_nodes;
182 return (NULL);
183 }
184 if (nodeid == nodemap[hi].nodeid) {
185 *pos = hi;
186 return (nodemap[hi].node);
187 }
188
189 while (1) {
190 *pos = (lo + hi) >> 1;
191 if (nodeid == nodemap[*pos].nodeid)
192 return (nodemap[*pos].node);
193 if (hi - lo <= 1)
194 break;
195 if (nodeid < nodemap[*pos].nodeid)
196 hi = *pos;
197 else
198 lo = *pos;
199 }
200 *pos = hi;
201 return (NULL);
202 }
203
204 /*---------- CompareNodes ----------------------------------------------
205 * compare two nodes, returns 0 if nodes are the same within
206 * the specified tolerance, else 1
207 *----------------------------------------------------------------------*/
208
CompareNodes(cgnsNODE * node1,cgnsNODE * node2)209 static int CompareNodes (cgnsNODE *node1, cgnsNODE *node2)
210 {
211 double dist = (node2->x - node1->x) * (node2->x - node1->x) +
212 (node2->y - node1->y) * (node2->y - node1->y) +
213 (node2->z - node1->z) * (node2->z - node1->z);
214 return (dist < (tolerance * tolerance) ? 0 : 1);
215 }
216
217 /*======================================================================
218 duplicate node checking routines
219
220 The nodes are stored in a sorted list based on radius from the origin.
221 Once the position in the list is determined, then a scan backwards and
222 forwards in the list is done for a matching node.
223 ========================================================================*/
224
225 /*---------- FindPosition ----------------------------------------------
226 * bisection search to locate position for node
227 *----------------------------------------------------------------------*/
228
FindPosition(cgnsNODE * node)229 static cgsize_t FindPosition (cgnsNODE *node)
230 {
231 cgsize_t mid, lo = 0, hi = num_entries - 1;
232
233 if (!num_entries || node->dist <= node_list[0]->dist)
234 return (0);
235
236 if (!hi || node->dist > node_list[hi]->dist)
237 return (num_entries);
238 if (node->dist == node_list[hi]->dist)
239 return (hi);
240
241 while ((hi - lo) > 1) {
242 mid = (lo + hi) >> 1;
243 if (node->dist == node_list[mid]->dist)
244 return (mid);
245 if (node->dist < node_list[mid]->dist)
246 hi = mid;
247 else
248 lo = mid;
249 }
250 return (hi);
251 }
252
253 /*---------- FindNode --------------------------------------------------
254 * search for matching node in Node List
255 *----------------------------------------------------------------------*/
256
FindNode(cgnsNODE * node,cgsize_t * pos)257 static cgnsNODE *FindNode (cgnsNODE *node, cgsize_t *pos)
258 {
259 cgsize_t n;
260
261 *pos = FindPosition (node);
262
263 for (n = *pos - 1; n >= 0; n--) {
264 if (fabs (node->dist - node_list[n]->dist) >= tolerance)
265 break;
266 if (!CompareNodes (node, node_list[n]))
267 return (node_list[n]);
268 }
269 for (n = *pos; n < num_entries; n++) {
270 if (fabs (node->dist - node_list[n]->dist) >= tolerance)
271 break;
272 if (!CompareNodes (node, node_list[n]))
273 return (node_list[n]);
274 }
275 return (NULL);
276 }
277
278 /*---------- AddNode ---------------------------------------------------
279 * add a new node to the duplicate node checking list
280 *----------------------------------------------------------------------*/
281
AddNode(cgnsNODE * node,cgsize_t pos)282 static void AddNode (cgnsNODE *node, cgsize_t pos)
283 {
284 cgsize_t n;
285
286 if (num_entries == max_entries) {
287 if (!max_entries)
288 node_list = (cgnsNODE **) malloc (ENTRY_INC * sizeof(cgnsNODE *));
289 else
290 node_list = (cgnsNODE **) realloc (node_list,
291 (size_t)(max_entries + ENTRY_INC) * sizeof(cgnsNODE *));
292 if (NULL == node_list)
293 cgnsImportFatal (
294 "AddNode:malloc failed for new node entry in duplicate node list");
295 max_entries += ENTRY_INC;
296 }
297 for (n = num_entries; n > pos; n--)
298 node_list[n] = node_list[n-1];
299 node_list[pos] = node;
300 num_entries++;
301 }
302
303 /*======================================================================
304 * Element routines
305 *======================================================================*/
306
307 /*---------- GetElement ------------------------------------------------
308 * return the element for a given element id
309 *----------------------------------------------------------------------*/
310
GetElement(cgsize_t elemid,cgsize_t * pos)311 static cgnsELEM *GetElement (cgsize_t elemid, cgsize_t *pos)
312 {
313 cgsize_t lo = 0, hi = num_elements - 1;
314
315 *pos = 0;
316 if (!num_elements || elemid < elemlist[0].elemid)
317 return (NULL);
318 if (elemid == elemlist[0].elemid)
319 return (&elemlist[0]);
320 if (!hi || elemid > elemlist[hi].elemid) {
321 *pos = num_elements;
322 return (NULL);
323 }
324 if (elemid == elemlist[hi].elemid) {
325 *pos = hi;
326 return (&elemlist[hi]);
327 }
328
329 while (1) {
330 *pos = (lo + hi) >> 1;
331 if (elemid == elemlist[*pos].elemid)
332 return (&elemlist[*pos]);
333 if (hi - lo <= 1)
334 break;
335 if (elemid < elemlist[*pos].elemid)
336 hi = *pos;
337 else
338 lo = *pos;
339 }
340 *pos = hi;
341 return (NULL);
342 }
343
344 /*---------- NewElement ------------------------------------------------
345 * add new element to list of elements
346 *----------------------------------------------------------------------*/
347
NewElement(cgsize_t pos)348 static cgnsELEM *NewElement (cgsize_t pos)
349 {
350 int i;
351 cgsize_t n;
352
353 /* malloc/realloc if needed */
354
355 if (num_elements == max_elements) {
356 if (!max_elements)
357 elemlist = (cgnsELEM *) malloc (ELEMENT_INC * sizeof(cgnsELEM));
358 else
359 elemlist = (cgnsELEM *) realloc (elemlist,
360 (size_t)(max_elements + ELEMENT_INC) * sizeof(cgnsELEM));
361 if (NULL == elemlist)
362 cgnsImportFatal ("AddElement:malloc failed for element list");
363 max_elements += ELEMENT_INC;
364 }
365
366 /* insert new element */
367
368 for (n = num_elements; n > pos; n--) {
369 elemlist[n].elemid = elemlist[n-1].elemid;
370 elemlist[n].elemtype = elemlist[n-1].elemtype;
371 elemlist[n].nodeid = elemlist[n-1].nodeid;
372 for (i = 0; i < 6; i++)
373 elemlist[n].facemap[i] = elemlist[n-1].facemap[i];
374 }
375 num_elements++;
376 return (&elemlist[pos]);
377 }
378
379 /*======================================================================
380 * Region routines
381 *======================================================================*/
382
383 /*---------- GetRegion -------------------------------------------------
384 * return a region for a given region name
385 *----------------------------------------------------------------------*/
386
GetRegion(char * name,int * pos)387 static cgnsREGN *GetRegion (char *name, int *pos)
388 {
389 int cmp, lo = 0, hi = num_regions - 1;
390
391 *pos = 0;
392 if (!num_regions || (cmp = strcmp (name, reglist[0].name)) < 0)
393 return (NULL);
394 if (0 == cmp)
395 return (®list[0]);
396 if (!hi || (cmp = strcmp (name, reglist[hi].name)) > 0) {
397 *pos = num_regions;
398 return (NULL);
399 }
400 if (0 == cmp) {
401 *pos = hi;
402 return (®list[hi]);
403 }
404
405 while (1) {
406 *pos = (lo + hi) >> 1;
407 if (0 == (cmp = strcmp (name, reglist[*pos].name)))
408 return (®list[*pos]);
409 if (hi - lo <= 1)
410 break;
411 if (cmp < 0)
412 hi = *pos;
413 else
414 lo = *pos;
415 }
416 *pos = hi;
417 return (NULL);
418 }
419
420 /*---------- NewRegion -------------------------------------------------
421 * add a new region to region list
422 *----------------------------------------------------------------------*/
423
NewRegion(char * name,int pos)424 static cgnsREGN *NewRegion (char *name, int pos)
425 {
426 int n;
427 static char *errmsg = "NewRegion:malloc failed for region list";
428
429 if (!num_regions)
430 reglist = (cgnsREGN *) malloc (sizeof(cgnsREGN));
431 else
432 reglist = (cgnsREGN *) realloc (reglist,
433 (num_regions + 1) * sizeof(cgnsREGN));
434 if (NULL == reglist)
435 cgnsImportFatal (errmsg);
436 for (n = num_regions; n > pos; n--) {
437 strcpy (reglist[n].name, reglist[n-1].name);
438 reglist[n].type = reglist[n-1].type;
439 reglist[n].nobjs = reglist[n-1].nobjs;
440 reglist[n].objid = reglist[n-1].objid;
441 }
442 strncpy (reglist[pos].name, name, 32);
443 reglist[pos].name[32] = 0;
444 reglist[pos].type = 0;
445 reglist[pos].nobjs = 0;
446 reglist[pos].objid = NULL;
447 num_regions++;
448 return (®list[pos]);
449 }
450
451 /*======================================================================
452 * external face regions
453 *======================================================================*/
454
455 /*---------- get_face_nodes -----------------------------------------
456 * get nodes for an element face
457 *-------------------------------------------------------------------*/
458
get_face_nodes(cgsize_t faceid,cgsize_t nodeid[4])459 static int get_face_nodes (cgsize_t faceid, cgsize_t nodeid[4])
460 {
461 cgnsELEM *elem;
462 cgsize_t elemid = faceid >> 3;
463 int facenum = (int)(faceid & 7);
464 int n, nfaces = 0, noff = 0, nnodes;
465 static int facenodes[20][5] = {
466 /* tet */
467 {3, 0, 2, 1, 0},
468 {3, 0, 1, 3, 0},
469 {3, 1, 2, 3, 0},
470 {3, 2, 0, 3, 0},
471 /* pyramid */
472 {4, 0, 3, 2, 1},
473 {3, 0, 1, 4, 0},
474 {3, 1, 2, 4, 0},
475 {3, 2, 3, 4, 0},
476 {3, 3, 0, 4, 0},
477 /* wedge */
478 {4, 0, 1, 4, 3},
479 {4, 1, 2, 5, 4},
480 {4, 2, 0, 3, 5},
481 {3, 0, 2, 1, 0},
482 {3, 3, 4, 5, 0},
483 /* hex */
484 {4, 0, 3, 2, 1},
485 {4, 0, 1, 5, 4},
486 {4, 1, 2, 6, 5},
487 {4, 2, 3, 7, 6},
488 {4, 0, 4, 7, 3},
489 {4, 4, 5, 6, 7}
490 };
491
492 if (elemid < 0 || elemid >= num_elements)
493 cgnsImportFatal ("get_face_nodes:invalid element number");
494 elem = &elemlist[elemid];
495 switch (elem->elemtype) {
496 case cgnsELEM_TET:
497 noff = 0;
498 nfaces = 4;
499 break;
500 case cgnsELEM_PYR:
501 noff = 4;
502 nfaces = 5;
503 break;
504 case cgnsELEM_WDG:
505 noff = 9;
506 nfaces = 5;
507 break;
508 case cgnsELEM_HEX:
509 noff = 14;
510 nfaces = 6;
511 break;
512 default:
513 cgnsImportFatal ("get_face_nodes:invalid element type");
514 }
515 if (facenum < 0 || facenum >= nfaces)
516 return (0);
517 noff += (int)elem->facemap[facenum];
518 nnodes = facenodes[noff][0];
519 for (n = 0; n < nnodes; n++)
520 nodeid[n] = elem->nodeid[facenodes[noff][n+1]];
521 return (nnodes);
522 }
523
524 /*---------- compare_faces -------------------------------------------
525 * face comparison routine
526 *--------------------------------------------------------------------*/
527
compare_faces(void * v1,void * v2)528 static int compare_faces (void *v1, void *v2)
529 {
530 cgnsFACE *f1 = (cgnsFACE *)v1;
531 cgnsFACE *f2 = (cgnsFACE *)v2;
532 int n;
533
534 if (f1->nnodes != f2->nnodes)
535 return (f1->nnodes - f2->nnodes);
536
537 /* the following assumes nodes have been sorted */
538
539 for (n = 0; n < f1->nnodes; n++) {
540 if (f1->nodeid[n] != f2->nodeid[n])
541 return (int)(f1->nodeid[n] - f2->nodeid[n]);
542 }
543 return (0);
544 }
545
546 /*---------- get_faces ----------------------------------------------
547 * get the exterior faces
548 *-------------------------------------------------------------------*/
549
get_faces(void * v)550 static void get_faces (void *v)
551 {
552 facelist[num_faces++] = (cgnsFACE *)v;
553 }
554
555 /*---------- hash_face -----------------------------------------------
556 * face hash function
557 *--------------------------------------------------------------------*/
558
hash_face(void * v)559 static size_t hash_face (void *v)
560 {
561 cgnsFACE *face = (cgnsFACE *)v;
562 int n;
563 size_t hash = 0;
564
565 for (n = 0; n < face->nnodes; n++)
566 hash += (size_t)face->nodeid[n];
567 return (hash);
568 }
569
570 /*---------- sortfaces -------------------------------------------------
571 * called by qsort to sort the list of faces
572 *----------------------------------------------------------------------*/
573
sortfaces(const void * f1,const void * f2)574 static int sortfaces (const void *f1, const void *f2)
575 {
576 return (int)((*((cgnsFACE **)f1))->faceid - (*((cgnsFACE **)f2))->faceid);
577 }
578
579 /*---------- exterior_faces -----------------------------------------
580 * find exterior faces
581 *-------------------------------------------------------------------*/
582
exterior_faces(void)583 static void exterior_faces (void)
584 {
585 int i, j, k, nfaces;
586 cgsize_t nodeid[4];
587 cgsize_t n, nn, id, faceid;
588 HASH FaceList;
589 cgnsFACE *pf, face;
590
591 FaceList = HashCreate (2047, compare_faces, hash_face);
592 if (NULL == FaceList)
593 cgnsImportFatal ("exterior_faces:malloc failed for face hash table");
594
595 for (n = 0; n < num_elements; n++) {
596 switch (elemlist[n].elemtype) {
597 case cgnsELEM_WDG:
598 nfaces = 5;
599 break;
600 case cgnsELEM_HEX:
601 nfaces = 6;
602 break;
603 default:
604 nfaces = elemlist[n].elemtype;
605 break;
606 }
607
608 /* loop over element faces */
609
610 for (j = 0; j < nfaces; j++) {
611
612 /* get face nodes and sort */
613
614 faceid = (n << 3) | j;
615 face.nnodes = get_face_nodes (faceid, nodeid);
616 for (i = 0; i < face.nnodes; i++) {
617 id = nodeid[i];
618 for (k = 0; k < i; k++) {
619 if (face.nodeid[k] > id) {
620 nn = face.nodeid[k];
621 face.nodeid[k] = id;
622 id = nn;
623 }
624 }
625 face.nodeid[i] = id;
626 }
627
628 if (NULL == (pf = (cgnsFACE *) HashFind (FaceList, &face))) {
629
630 /* create new face and add to list */
631
632 if (NULL == (pf = (cgnsFACE *) malloc (sizeof(cgnsFACE))))
633 cgnsImportFatal ("exterior_faces:malloc failed for new face");
634 pf->faceid = faceid;
635 pf->flags = 0;
636 pf->nnodes = face.nnodes;
637 for (i = 0; i < face.nnodes; i++)
638 pf->nodeid[i] = face.nodeid[i];
639 (void) HashAdd (FaceList, pf);
640 }
641
642 /* else already exists */
643
644 else {
645 HashDelete (FaceList, pf);
646 free (pf);
647 }
648 }
649 }
650
651 facelist = (cgnsFACE **) malloc (HashSize (FaceList) * sizeof(cgnsFACE *));
652 if (NULL == facelist)
653 cgnsImportFatal ("exterior_faces:malloc failed for exterior face list");
654 num_faces = 0;
655 HashDestroy (FaceList, get_faces);
656
657 /* check if faces need sorting */
658
659 for (n = 1; n < num_faces; n++) {
660 if (facelist[n]->faceid < facelist[n-1]->faceid)
661 break;
662 }
663 if (n < num_faces)
664 qsort (facelist, (size_t)num_faces, sizeof(cgnsFACE *), sortfaces);
665
666 /* get face nodes in the correct order */
667
668 for (n = 0; n < num_faces; n++) {
669 get_face_nodes (facelist[n]->faceid, nodeid);
670 for (i = 0; i < 4; i++)
671 facelist[n]->nodeid[i] = nodeid[i];
672 }
673 }
674
675 /*===================================================================
676 * write regions to cgns file
677 *===================================================================*/
678
679 /*---------- sortnodes -------------------------------------------------
680 * called by qsort to sort list of node ID mappings
681 *----------------------------------------------------------------------*/
682
sortnodes(const void * n1,const void * n2)683 static int sortnodes (const void *n1, const void *n2)
684 {
685 return (int)(*((cgsize_t *)n1) - *((cgsize_t *)n2));
686 }
687
688 /*---------- write_node_region --------------------------------------
689 * write region from node list
690 *-------------------------------------------------------------------*/
691
write_node_region(cgnsREGN * reg,cgsize_t offset)692 static cgsize_t write_node_region (cgnsREGN *reg, cgsize_t offset)
693 {
694 int nn, isect;
695 cgsize_t i, j, mid, lo, hi, pos;
696 cgsize_t nfaces, nc, *conns, *conns_offset;
697 CGNS_ENUMT(ElementType_t) elemtype = CGNS_ENUMV(ElementTypeNull);
698 cgnsNODE *node;
699
700 /* get exterior faces */
701
702 if (num_faces == 0) exterior_faces ();
703 for (j = 0; j < num_faces; j++)
704 facelist[j]->flags = 0;
705
706 /* sort region nodes */
707
708 for (i = 1; i < reg->nobjs; i++) {
709 if (reg->objid[i] < reg->objid[i-1])
710 break;
711 }
712 if (i < reg->nobjs)
713 qsort (reg->objid, (size_t)reg->nobjs, sizeof(cgsize_t), sortnodes);
714
715 /* scan list of exterior faces */
716
717 nfaces = nc = 0;
718 for (j = 0; j < num_faces; j++) {
719 if (facelist[j]->flags) continue;
720 for (nn = 0; nn < facelist[j]->nnodes; nn++) {
721 lo = 0;
722 hi = reg->nobjs - 1;
723 while (lo <= hi) {
724 mid = (lo + hi) >> 1;
725 if (facelist[j]->nodeid[nn] == reg->objid[mid])
726 break;
727 if (facelist[j]->nodeid[nn] < reg->objid[mid])
728 hi = mid - 1;
729 else
730 lo = mid + 1;
731 }
732 if (lo > hi)
733 break;
734 }
735 if (nn == facelist[j]->nnodes) {
736 nfaces++;
737 facelist[j]->flags = 1;
738 if (nc != nn) {
739 elemtype = nc ? CGNS_ENUMV(MIXED) :
740 (nn == 3 ? CGNS_ENUMV(TRI_3) : CGNS_ENUMV(QUAD_4));
741 nc = nn;
742 }
743 }
744 }
745 if (!nfaces) return 0;
746
747 conns = (cgsize_t *) malloc ((size_t)(5 * nfaces) * sizeof(cgsize_t));
748 if (NULL == conns)
749 cgnsImportFatal ("write_node_region:malloc failed for connectivity");
750
751 if (elemtype == CGNS_ENUMV(MIXED)){
752 conns_offset =(cgsize_t *) malloc ((size_t)(num_faces+1) * sizeof(cgsize_t));
753 if (NULL == conns_offset)
754 cgnsImportFatal ("write_node_region:malloc failed for connectivity offset");
755 conns_offset[0] = 0;
756 }
757
758 /* write face connectivities */
759
760 for (nc = 0, j = 0, i = 0; j < num_faces; j++) {
761 if (facelist[j]->flags) {
762 if (elemtype == CGNS_ENUMV(MIXED)){
763 conns[nc++] = facelist[j]->nnodes == 3 ?
764 CGNS_ENUMV(TRI_3) : CGNS_ENUMV(QUAD_4);
765 conns_offset[i+1] = conns_offset[i] + conns[nc-1] + 1;
766 i++;
767 }
768 for (nn = 0; nn < facelist[j]->nnodes; nn++) {
769 if (NULL == (node = GetNode (facelist[j]->nodeid[nn], &pos)))
770 cgnsImportFatal ("write_node_region:missing element node");
771 conns[nc++] = pos + 1;
772 }
773 }
774 }
775
776 if (elemtype == CGNS_ENUMV(MIXED)) {
777 if (cg_poly_section_write (cgnsFile, cgnsBase, cgnsZone, reg->name,
778 elemtype, offset, offset + nfaces - 1, 0, conns, conns_offset, &isect))
779 cgnsImportFatal ((char *)cg_get_error());
780 free(conns_offset);
781 }
782 else {
783 if (cg_section_write (cgnsFile, cgnsBase, cgnsZone, reg->name,
784 elemtype, offset, offset + nfaces - 1, 0, conns, &isect))
785 cgnsImportFatal ((char *)cg_get_error());
786 }
787
788 /* create parent cell mapping */
789
790 for (nc = 0, j = 0; j < num_faces; j++) {
791 if (facelist[j]->flags)
792 conns[nc++] = (facelist[j]->faceid >> 3) + 1;
793 }
794 for (j = 0; j < nfaces; j++)
795 conns[nc++] = 0;
796 for (j = 0; j < num_faces; j++) {
797 if (facelist[j]->flags)
798 conns[nc++] = (facelist[j]->faceid & 7) + 1;
799 }
800 for (j = 0; j < nfaces; j++)
801 conns[nc++] = 0;
802 if (cg_parent_data_write (cgnsFile, cgnsBase, cgnsZone, isect, conns))
803 cgnsImportFatal ((char *)cg_get_error());
804
805 free (conns);
806 return nfaces;
807 }
808
809 /*---------- write_face_region --------------------------------------
810 * write region from face list
811 *-------------------------------------------------------------------*/
812
write_face_region(cgnsREGN * reg,cgsize_t offset)813 static cgsize_t write_face_region (cgnsREGN *reg, cgsize_t offset)
814 {
815 int nn, facenum, i, isect;
816 cgsize_t elemid, nodeid[4];
817 cgsize_t n, nc, pos, *conns, *conns_offset;
818 CGNS_ENUMT(ElementType_t) elemtype = CGNS_ENUMV(ElementTypeNull);
819 cgnsELEM *elem;
820 cgnsNODE *node;
821
822 if (!reg->nobjs) return 0;
823 conns = (cgsize_t *) malloc ((size_t)(5 * reg->nobjs) * sizeof(cgsize_t));
824 if (NULL == conns)
825 cgnsImportFatal ("write_face_region:malloc failed for connectivity");
826
827 for (i = 0, n = 0; n < reg->nobjs; n++) {
828 elemid = reg->objid[n] >> 3;
829 facenum = (int)(reg->objid[n] & 7) - 1;
830 if (NULL == (elem = GetElement (elemid, &pos)))
831 cgnsImportFatal ("write_face_region:region element not found");
832 nn = get_face_nodes ((pos << 3) | facenum, nodeid);
833 if (i != nn) {
834 if (i) {
835 elemtype = CGNS_ENUMV(MIXED);
836 break;
837 }
838 i = nn;
839 elemtype = nn == 3 ? CGNS_ENUMV(TRI_3) : CGNS_ENUMV(QUAD_4);
840 }
841 }
842
843 if (elemtype == CGNS_ENUMV(MIXED))
844 {
845 conns_offset = (cgsize_t *) malloc ((size_t)(reg->nobjs+1) * sizeof(cgsize_t));
846 if (NULL == conns_offset)
847 cgnsImportFatal ("write_face_region:malloc failed for connectivity offset");
848 conns_offset[0] = 0;
849 }
850
851 for (nc = 0, n = 0; n < reg->nobjs; n++) {
852 elemid = reg->objid[n] >> 3;
853 facenum = (int)(reg->objid[n] & 7) - 1;
854 elem = GetElement (elemid, &pos);
855 nn = get_face_nodes ((pos << 3) | facenum, nodeid);
856 if (elemtype == CGNS_ENUMV(MIXED)) {
857 conns[nc++] = nn == 3 ? CGNS_ENUMV(TRI_3) : CGNS_ENUMV(QUAD_4);
858 conns_offset[n+1] = conns_offset[n] + conns[nc-1] + 1;
859 }
860 for (i = 0; i < nn; i++) {
861 if (NULL == (node = GetNode (nodeid[i], &pos)))
862 cgnsImportFatal ("write_face_region:missing element node");
863 conns[nc++] = pos + i;
864 }
865 }
866
867 if (elemtype == CGNS_ENUMV(MIXED))
868 {
869 if (cg_poly_section_write (cgnsFile, cgnsBase, cgnsZone, reg->name,
870 elemtype, offset, offset + reg->nobjs - 1, 0, conns, conns_offset, &isect))
871 cgnsImportFatal ((char *)cg_get_error());
872 free(conns_offset);
873 }
874 else {
875 if (cg_section_write (cgnsFile, cgnsBase, cgnsZone, reg->name,
876 elemtype, offset, offset + reg->nobjs - 1, 0, conns, &isect))
877 cgnsImportFatal ((char *)cg_get_error());
878 }
879
880 free (conns);
881 return reg->nobjs;
882 }
883
884 /*---------- write_elem_region --------------------------------------
885 * write elements as region
886 *-------------------------------------------------------------------*/
887
write_elem_region(cgnsREGN * reg,cgsize_t offset)888 static cgsize_t write_elem_region (cgnsREGN *reg, cgsize_t offset)
889 {
890 return 0;
891 }
892
893 /*======================================================================
894 * API routines
895 *======================================================================*/
896
897 /*---------- cgnsImportError -------------------------------------------
898 * setup error handler call back
899 *----------------------------------------------------------------------*/
900
cgnsImportError(void (* callback)(char * msg))901 void cgnsImportError (void (*callback)(char *msg))
902 {
903 errcallback = callback;
904 }
905
906 /*---------- cgnsImportFatal -------------------------------------------
907 * write error message and exit
908 *----------------------------------------------------------------------*/
909
cgnsImportFatal(char * errmsg)910 void cgnsImportFatal (char *errmsg)
911 {
912 if (NULL != errcallback)
913 (*errcallback) (errmsg);
914 else if (NULL != errmsg && *errmsg)
915 fprintf (stderr, "%s\n", errmsg);
916 exit (-1);
917 }
918
919 /*---------- cgnsImportSetTol ------------------------------------------
920 * setup tolerance for duplicate node checking
921 *----------------------------------------------------------------------*/
922
cgnsImportSetTol(double tol)923 double cgnsImportSetTol (double tol)
924 {
925 tolerance = tol >= 0.0 ? tol : TOLERANCE;
926 tol = def_tol;
927 def_tol = tolerance;
928 return (tol);
929 }
930
931 /*---------- cgnsImportGetTol ------------------------------------------
932 * return tolerance for duplicate node checking
933 *----------------------------------------------------------------------*/
934
cgnsImportGetTol(int rel)935 double cgnsImportGetTol (int rel)
936 {
937 double tol = def_tol;
938
939 if (rel && num_nodes) {
940 double avgvol = (xmax-xmin) * (ymax-ymin) * (zmax-zmin) /
941 (DOUBLE)num_nodes;
942 tol *= pow (avgvol, 0.33333);
943 }
944 return (tol);
945 }
946
947 /*---------- cgnsImportSetCheck ----------------------------------------
948 * set duplicate node checking on/off
949 *----------------------------------------------------------------------*/
950
cgnsImportSetCheck(int set)951 void cgnsImportSetCheck (int set)
952 {
953 no_check = set;
954 }
955
956 /*---------- cgnsImportRange --------------------------------------------
957 * gets bounding box of node coordinates
958 *-----------------------------------------------------------------------*/
959
cgnsImportRange(double * x1,double * y1,double * z1,double * x2,double * y2,double * z2)960 cgsize_t cgnsImportRange (double *x1, double *y1, double *z1,
961 double *x2, double *y2, double *z2)
962 {
963 *x1 = xmin; *y1 = ymin; *z1 = zmin;
964 *x2 = xmax; *y2 = ymax; *z2 = zmax;
965 return num_nodes;
966 }
967
968 /*---------- cgnsImportCheck --------------------------------------------
969 * check for and remove duplicate nodes
970 *-----------------------------------------------------------------------*/
971
cgnsImportCheck(int rel)972 cgsize_t cgnsImportCheck (int rel)
973 {
974 cgsize_t n, pos, dup_cnt = 0;
975 cgnsNODE *node;
976
977 if (num_nodes < 2)
978 return (0);
979
980 /* set tolerance */
981
982 tolerance = cgnsImportGetTol (rel);
983
984 /* scan list of nodes, and remove duplicates */
985
986 for (n = 0; n < num_nodes; n++) {
987 if (!nodemap[n].mapped) {
988 nodemap[n].node->dist = sqrt ((double)(
989 (nodemap[n].node->x - xmin) * (nodemap[n].node->x - xmin) +
990 (nodemap[n].node->y - ymin) * (nodemap[n].node->y - ymin) +
991 (nodemap[n].node->z - zmin) * (nodemap[n].node->z - zmin)));
992 node = FindNode (nodemap[n].node, &pos);
993 if (NULL == node)
994 AddNode (nodemap[n].node, pos);
995 else if (node != nodemap[n].node) {
996 if (REFS_BIT == (nodemap[n].node->flags & REFS_BIT)) {
997 for (pos = 0; pos < num_nodes; pos++) {
998 if (nodemap[pos].mapped &&
999 nodemap[pos].node == nodemap[n].node)
1000 nodemap[pos].node = node;
1001 }
1002 }
1003 node->flags |=
1004 ((nodemap[n].node->flags & USED_BIT) | REFS_BIT);
1005 free (nodemap[n].node);
1006 nodemap[n].node = node;
1007 dup_cnt++;
1008 }
1009 nodemap[n].mapped = 1;
1010 }
1011 }
1012
1013 /* free duplicate node checking list */
1014
1015 free (node_list);
1016 num_entries = max_entries = 0;
1017 return (dup_cnt);
1018 }
1019
1020 /*---------- cgnsImportMap ----------------------------------------------
1021 * map a node explicitly to another
1022 *-----------------------------------------------------------------------*/
1023
cgnsImportMap(cgsize_t nodeid,cgsize_t mapid)1024 cgsize_t cgnsImportMap (cgsize_t nodeid, cgsize_t mapid)
1025 {
1026 cgsize_t p1, p2, ret;
1027 cgnsNODE *n1 = GetNode (nodeid, &p1);
1028 cgnsNODE *n2 = GetNode (mapid, &p2);
1029
1030 if (NULL == n1 || NULL == n2)
1031 return (0);
1032 if (n1 == n2)
1033 return (n1->id);
1034 ret = CompareNodes (n1, n2) ? -1 : n1->id;
1035 if (REFS_BIT == (n2->flags & REFS_BIT)) {
1036 cgsize_t n;
1037 for (n = 0; n < num_nodes; n++) {
1038 if (nodemap[n].node == n2) {
1039 nodemap[n].node = n1;
1040 nodemap[n].mapped = 1;
1041 }
1042 }
1043 }
1044 else {
1045 nodemap[p2].node = n1;
1046 nodemap[p2].mapped = 1;
1047 }
1048 n1->flags |= ((n2->flags & USED_BIT) | REFS_BIT);
1049 free (n2);
1050 return (ret);
1051 }
1052
1053 /*---------- cgnsImportNode ---------------------------------------------
1054 * import a node
1055 *-----------------------------------------------------------------------*/
1056
cgnsImportNode(cgsize_t nodeid,double x,double y,double z)1057 cgsize_t cgnsImportNode (cgsize_t nodeid, double x, double y, double z)
1058 {
1059 cgsize_t n, pos;
1060
1061 if (nodeid <= 0)
1062 return (0);
1063
1064 /* find min/max bounds */
1065
1066 if (!num_nodes) {
1067 xmin = xmax = x;
1068 ymin = ymax = y;
1069 zmin = zmax = z;
1070 }
1071 else {
1072 if (xmin > x) xmin = x;
1073 if (xmax < x) xmax = x;
1074 if (ymin > y) ymin = y;
1075 if (ymax < y) ymax = y;
1076 if (zmin > z) zmin = z;
1077 if (zmax < z) zmax = z;
1078 }
1079
1080 /* find position to place node id */
1081
1082 if (NULL != GetNode (nodeid, &pos)) {
1083 nodemap[pos].node->x = (DOUBLE)x;
1084 nodemap[pos].node->y = (DOUBLE)y;
1085 nodemap[pos].node->z = (DOUBLE)z;
1086 return (-1);
1087 }
1088
1089 /* malloc/realloc if needed */
1090
1091 if (num_nodes == max_nodes) {
1092 if (!max_nodes)
1093 nodemap = (NODEMAP *) malloc (NODEMAP_INC * sizeof(NODEMAP));
1094 else
1095 nodemap = (NODEMAP *) realloc (nodemap,
1096 (size_t)(max_nodes + NODEMAP_INC) * sizeof(NODEMAP));
1097 if (NULL == nodemap)
1098 cgnsImportFatal (
1099 "cgnsImportNode:malloc failed for node mapping data");
1100 max_nodes += NODEMAP_INC;
1101 }
1102
1103 /* insert new node */
1104
1105 for (n = num_nodes; n > pos; n--) {
1106 nodemap[n].nodeid = nodemap[n-1].nodeid;
1107 nodemap[n].mapped = nodemap[n-1].mapped;
1108 nodemap[n].node = nodemap[n-1].node;
1109 }
1110 nodemap[pos].nodeid = nodeid;
1111 nodemap[pos].mapped = no_check;
1112 nodemap[pos].node = NewNode (nodeid, x, y, z);
1113 if (NULL == nodemap[pos].node)
1114 cgnsImportFatal ("cgnsImportNode:malloc failed for a new node");
1115 num_nodes++;
1116 return (nodeid);
1117 }
1118
1119 /*---------- cgnsImportSetNode ------------------------------------------
1120 * set node coordinates for node ID
1121 *-----------------------------------------------------------------------*/
1122
cgnsImportSetNode(cgsize_t nodeid,double x,double y,double z)1123 cgsize_t cgnsImportSetNode (cgsize_t nodeid, double x, double y, double z)
1124 {
1125 cgsize_t n;
1126 cgnsNODE *node = GetNode (nodeid, &n);
1127
1128 if (NULL != node) {
1129 node->x = (DOUBLE)x;
1130 node->y = (DOUBLE)y;
1131 node->z = (DOUBLE)z;
1132 return (node->id);
1133 }
1134 return (0);
1135 }
1136
1137 /*---------- cgnsImportGetNode ------------------------------------------
1138 * return node coordinates for node ID
1139 *-----------------------------------------------------------------------*/
1140
cgnsImportGetNode(cgsize_t nodeid,double * x,double * y,double * z)1141 cgsize_t cgnsImportGetNode (cgsize_t nodeid, double *x, double *y, double *z)
1142 {
1143 cgsize_t n;
1144 cgnsNODE *node = GetNode (nodeid, &n);
1145
1146 if (NULL != node) {
1147 *x = node->x;
1148 *y = node->y;
1149 *z = node->z;
1150 return (node->id);
1151 }
1152 return (0);
1153 }
1154
1155 /*---------- cgnsImportNodeList -----------------------------------------
1156 * return list of all node ID's
1157 *-----------------------------------------------------------------------*/
1158
cgnsImportNodeList(void)1159 cgsize_t *cgnsImportNodeList (void)
1160 {
1161 cgsize_t n, *nodeids;
1162
1163 nodeids = (cgsize_t *) malloc ((size_t)(num_nodes + 1) * sizeof(cgsize_t));
1164 if (NULL == nodeids)
1165 cgnsImportFatal ("cgnsImportNodeList:malloc failed for node ID list");
1166 nodeids[0] = num_nodes;
1167 for (n = 0; n < num_nodes; n++)
1168 nodeids[n+1] = nodemap[n].nodeid;
1169 return (nodeids);
1170 }
1171
1172 /*---------- cgnsImportAddVariable -------------------------------------
1173 * create a node variable
1174 *----------------------------------------------------------------------*/
1175
cgnsImportAddVariable(char * varname)1176 int cgnsImportAddVariable (char *varname)
1177 {
1178 int n;
1179 cgnsNODE *node;
1180
1181 if (varname == NULL || !*varname) return -1;
1182 for (n = 0; n < num_vars; n++) {
1183 if (0 == strcmp(varname, var_names[n])) return n;
1184 }
1185 if (num_vars)
1186 var_names = (char **) realloc (var_names, (num_vars + 1) * sizeof(char *));
1187 else
1188 var_names = (char **) malloc (sizeof(char *));
1189 if (var_names == NULL)
1190 cgnsImportFatal ("AddVariable:malloc/realloc failed for variable name list");
1191 var_names[num_vars] = (char *) malloc (strlen(varname) + 1);
1192 if (var_names[num_vars] == NULL)
1193 cgnsImportFatal ("AddVariable:malloc failed for variable name");
1194 strcpy(var_names[num_vars++], varname);
1195
1196 for (n = 0; n < num_entries; n++) {
1197 node = node_list[n];
1198 if (num_vars == 1)
1199 node->vars = (DOUBLE *) malloc (sizeof(DOUBLE));
1200 else
1201 node->vars = (DOUBLE *) realloc (node->vars, num_vars * sizeof(DOUBLE));
1202 if (node->vars == NULL)
1203 cgnsImportFatal ("AddVariable:malloc failed for node variables");
1204 node->vars[num_vars-1] = 0.0;
1205 }
1206 return num_vars-1;
1207 }
1208
1209 /*---------- cgnsImportGetVariable -------------------------------------
1210 * get the variable number for a node variable
1211 *----------------------------------------------------------------------*/
1212
cgnsImportGetVariable(char * varname)1213 int cgnsImportGetVariable (char *varname)
1214 {
1215 int n;
1216
1217 if (varname != NULL && *varname) {
1218 for (n = 0; n < num_vars; n++) {
1219 if (0 == strcmp(varname, var_names[n])) return n;
1220 }
1221 }
1222 return -1;
1223 }
1224
1225 /*---------- cgnsImportVariable ----------------------------------------
1226 * set the value of a variable at a node
1227 *----------------------------------------------------------------------*/
1228
cgnsImportVariable(cgsize_t nodeid,int varnum,double val)1229 cgsize_t cgnsImportVariable (cgsize_t nodeid, int varnum, double val)
1230 {
1231 cgsize_t n;
1232 cgnsNODE *node = GetNode (nodeid, &n);
1233
1234 if (NULL == node || varnum < 0 || varnum >= num_vars) return 0;
1235 node->vars[varnum] = (DOUBLE)val;
1236 return node->id;
1237 }
1238
1239 /*---------- cgnsImportElement ------------------------------------------
1240 * import an element
1241 *-----------------------------------------------------------------------*/
1242
cgnsImportElement(cgsize_t elemid,int elemtype,cgsize_t * nodelist)1243 int cgnsImportElement (cgsize_t elemid, int elemtype, cgsize_t *nodelist)
1244 {
1245 int n, ret;
1246 cgsize_t pos;
1247 cgnsNODE *node;
1248 cgnsELEM *elem;
1249
1250 if (elemid <= 0 ||
1251 (elemtype != cgnsELEM_TET && elemtype != cgnsELEM_PYR &&
1252 elemtype != cgnsELEM_WDG && elemtype != cgnsELEM_HEX))
1253 return (0);
1254
1255 /* element not found */
1256
1257 if (NULL == (elem = GetElement (elemid, &pos))) {
1258 ret = elemtype;
1259 elem = NewElement (pos);
1260 }
1261
1262 /* element already exists */
1263
1264 else {
1265 ret = -1;
1266 free (elem->nodeid);
1267 }
1268
1269 /* set element values */
1270
1271 elem->elemid = elemid;
1272 elem->elemtype = elemtype;
1273 elem->nodeid = (cgsize_t *) malloc (elemtype * sizeof(cgsize_t));
1274 if (NULL == elem->nodeid)
1275 cgnsImportFatal (
1276 "cgnsImportElement:malloc failed for a new element");
1277
1278 for (n = 0; n < elemtype; n++) {
1279 if (NULL == (node = GetNode (nodelist[n], &pos))) {
1280 char errmsg[50];
1281 sprintf (errmsg, "cgnsImportElement:element node %ld not found",
1282 (long)nodelist[n]);
1283 cgnsImportFatal (errmsg);
1284 }
1285 elem->nodeid[n] = node->id;
1286 node->flags |= USED_BIT;
1287 }
1288 for (n = 0; n < 6; n++)
1289 elem->facemap[n] = n;
1290
1291 return (ret);
1292 }
1293
1294 /*---------- cgnsImportGetElement ---------------------------------------
1295 * return element for element ID
1296 *-----------------------------------------------------------------------*/
1297
cgnsImportGetElement(cgsize_t elemid,cgsize_t nodeid[])1298 int cgnsImportGetElement (cgsize_t elemid, cgsize_t nodeid[])
1299 {
1300 int n;
1301 cgsize_t pos;
1302 cgnsELEM *elem = GetElement (elemid, &pos);
1303
1304 if (NULL != elem) {
1305 for (n = 0; n < elem->elemtype; n++)
1306 nodeid[n] = elem->nodeid[n];
1307 return (elem->elemtype);
1308 }
1309 return (0);
1310 }
1311
1312 /*---------- cgnsImportElementList --------------------------------------
1313 * return list of all element ID's
1314 *-----------------------------------------------------------------------*/
1315
cgnsImportElementList(void)1316 cgsize_t *cgnsImportElementList (void)
1317 {
1318 cgsize_t n, *elemids;
1319
1320 elemids = (cgsize_t *) malloc ((size_t)(num_elements + 1) * sizeof(int));
1321 if (NULL == elemids)
1322 cgnsImportFatal (
1323 "cgnsImportElementList:malloc failed for element ID list");
1324 elemids[0] = num_elements;
1325 for (n = 0; n < num_elements; n++)
1326 elemids[n+1] = elemlist[n].elemid;
1327 return (elemids);
1328 }
1329
1330 /*---------- cgnsImportGetFace ------------------------------------------
1331 * return element face node ID's
1332 *-----------------------------------------------------------------------*/
1333
cgnsImportGetFace(cgsize_t elemid,int facenum,cgsize_t nodeid[])1334 int cgnsImportGetFace (cgsize_t elemid, int facenum, cgsize_t nodeid[])
1335 {
1336 int nfaces;
1337 cgsize_t pos;
1338 cgnsELEM *elem = GetElement (elemid, &pos);
1339
1340 if (NULL == elem)
1341 return (-1);
1342 switch (elem->elemtype) {
1343 case cgnsELEM_WDG:
1344 nfaces = 5;
1345 break;
1346 case cgnsELEM_HEX:
1347 nfaces = 6;
1348 break;
1349 default:
1350 nfaces = elem->elemtype;
1351 break;
1352 }
1353 if (--facenum < 0 || facenum >= nfaces)
1354 return (0);
1355 return get_face_nodes ((pos << 3) | facenum, nodeid);
1356 }
1357
1358 /*---------- cgnsImportFindFace -----------------------------------------
1359 * return element face number given face node ID's
1360 *-----------------------------------------------------------------------*/
1361
cgnsImportFindFace(cgsize_t elemid,int nnodes,cgsize_t nodeid[])1362 int cgnsImportFindFace (cgsize_t elemid, int nnodes, cgsize_t nodeid[])
1363 {
1364 int i, j, nfaces = 0, noff = 0, mask = 0;
1365 cgsize_t pos;
1366 cgnsELEM *elem = GetElement (elemid, &pos);
1367 static int facemask[4][6] = {
1368 /* tet */
1369 { 7, 11, 14, 13, 0, 0},
1370 /* pyramid */
1371 {15, 19, 22, 28, 25, 0},
1372 /* wedge */
1373 {27, 54, 45, 7, 56, 0},
1374 /* hex */
1375 {15, 51, 102, 204, 153, 240}
1376 };
1377
1378 if (NULL == elem || NULL == nodeid)
1379 return (-1);
1380
1381 switch (elem->elemtype) {
1382 case cgnsELEM_TET:
1383 if (nnodes != 3)
1384 return (-1);
1385 noff = 0;
1386 nfaces = 4;
1387 break;
1388 case cgnsELEM_PYR:
1389 if (nnodes < 3 || nnodes > 4)
1390 return (-1);
1391 noff = 1;
1392 nfaces = 5;
1393 break;
1394 case cgnsELEM_WDG:
1395 if (nnodes < 3 || nnodes > 4)
1396 return (-1);
1397 noff = 2;
1398 nfaces = 5;
1399 break;
1400 case cgnsELEM_HEX:
1401 if (nnodes != 4)
1402 return (-1);
1403 noff = 3;
1404 nfaces = 6;
1405 break;
1406 default:
1407 cgnsImportFatal ("cgnsImportFindFace:invalid element type");
1408 }
1409
1410 for (j = 0; j < nnodes; j++) {
1411 for (i = 0; i < elem->elemtype; i++) {
1412 if (nodeid[j] == elem->nodeid[i])
1413 break;
1414 }
1415 if (i == elem->elemtype)
1416 return (0);
1417 mask |= (1 << i);
1418 }
1419 for (i = 0; i < nfaces; i++) {
1420 if (mask == facemask[noff][i]) {
1421 for (j = 0; j < 6; j++) {
1422 if (i == (int)elem->facemap[j])
1423 return (j + 1);
1424 }
1425 }
1426 }
1427 return (0);
1428 }
1429
1430 /*---------- cgnsImportBegReg --------------------------------------------
1431 * begin a region specification
1432 *-----------------------------------------------------------------------*/
1433
cgnsImportBegReg(char * regname,int regtype)1434 cgsize_t cgnsImportBegReg (char *regname, int regtype)
1435 {
1436 int n;
1437 cgnsREGN *reg;
1438
1439 if (region_id)
1440 cgnsImportEndReg ();
1441
1442 /* initialize region node list */
1443
1444 if (0 == region_max) {
1445 region_list = (cgsize_t *) malloc (REGION_INC * sizeof(cgsize_t));
1446 if (NULL == region_list)
1447 cgnsImportFatal (
1448 "cgnsImportBegReg:malloc failed for region node list");
1449 region_max = REGION_INC;
1450 }
1451
1452 /* initialize region data */
1453
1454 region_id = num_regions + 1;
1455 region_type = regtype;
1456 if (NULL == regname || !*regname)
1457 sprintf (region_name, "%s%ld", REGION_BASE, (long)region_id);
1458 else {
1459 strncpy (region_name, regname, sizeof(region_name));
1460 region_name[sizeof(region_name)-1] = 0;
1461 }
1462 region_nobjs = 0;
1463
1464 if (NULL == (reg = GetRegion (region_name, &n)))
1465 return (0);
1466 if (reg->type != regtype)
1467 cgnsImportFatal ("cgnsImportBegReg:only 1 type allowed for a region");
1468 return (reg->nobjs);
1469 }
1470
1471 /*---------- cgnsImportAddReg -------------------------------------------
1472 * add nodes to the region
1473 *-----------------------------------------------------------------------*/
1474
cgnsImportAddReg(cgsize_t numobjs,cgsize_t * objlist)1475 cgsize_t cgnsImportAddReg (cgsize_t numobjs, cgsize_t *objlist)
1476 {
1477 cgsize_t n, pos;
1478 char errmsg[50];
1479
1480 if (!region_id)
1481 cgnsImportFatal ("cgnsImportAddReg:region not defined");
1482
1483 /* realloc region list array if needed */
1484
1485 if (region_nobjs + numobjs > region_max) {
1486 n = region_nobjs + numobjs - region_max;
1487 if (n < REGION_INC) n = REGION_INC;
1488 region_list = (cgsize_t *) realloc (region_list,
1489 (size_t)(region_max + n) * sizeof(cgsize_t));
1490 if (NULL == region_list)
1491 cgnsImportFatal (
1492 "cgnsImportAddReg:malloc failed for region node list");
1493 region_max += n;
1494 }
1495
1496 /* node region */
1497
1498 if (region_type == cgnsREG_NODES) {
1499 cgnsNODE *node;
1500 for (n = 0; n < numobjs; n++) {
1501 if (NULL == (node = GetNode (objlist[n], &pos))) {
1502 sprintf (errmsg, "cgnsImportAddReg:region node %ld not found",
1503 (long)objlist[n]);
1504 cgnsImportFatal (errmsg);
1505 }
1506 region_list[region_nobjs++] = node->id;
1507 node->flags |= USED_BIT;
1508 }
1509 }
1510
1511 /* face region */
1512
1513 else if (region_type == cgnsREG_FACES) {
1514 int facenum, nfaces;
1515 cgsize_t elemid;
1516 cgnsELEM *elem;
1517 for (n = 0; n < numobjs; n++) {
1518 elemid = objlist[n] >> 3;
1519 facenum = (int)(objlist[n] & 7);
1520 if (NULL == (elem = GetElement (elemid, &pos))) {
1521 sprintf (errmsg, "cgnsImportAddReg:region element %ld not found",
1522 (long)elemid);
1523 cgnsImportFatal (errmsg);
1524 }
1525 if (elem->elemtype == cgnsELEM_WDG)
1526 nfaces = 5;
1527 else if (elem->elemtype == cgnsELEM_HEX)
1528 nfaces = 6;
1529 else
1530 nfaces = elem->elemtype;
1531 if (facenum < 1 || facenum > nfaces)
1532 cgnsImportFatal ("cgnsImportAddReg:region face number out of range");
1533 region_list[region_nobjs++] = objlist[n];
1534 }
1535 }
1536
1537 /* element region */
1538
1539 else if (region_type == cgnsREG_ELEMS) {
1540 for (n = 0; n < numobjs; n++) {
1541 if (NULL == GetElement (objlist[n], &pos)) {
1542 sprintf (errmsg, "cgnsImportAddReg:region element %ld not found",
1543 (long)objlist[n]);
1544 cgnsImportFatal (errmsg);
1545 }
1546 region_list[region_nobjs++] = objlist[n];
1547 }
1548 }
1549
1550 else
1551 cgnsImportFatal ("cgnsImportAddReg:undefined region type");
1552
1553 return (region_nobjs);
1554 }
1555
1556 /*---------- cgnsImportEndReg -------------------------------------------
1557 * end region definition and import region
1558 *-----------------------------------------------------------------------*/
1559
cgnsImportEndReg(void)1560 cgsize_t cgnsImportEndReg (void)
1561 {
1562 int pos;
1563 cgsize_t n;
1564 cgnsREGN *reg;
1565
1566 if (!region_id || !region_nobjs)
1567 return (region_id = 0);
1568 region_id = 0;
1569
1570 /* create a new region */
1571
1572 if (NULL == (reg = GetRegion (region_name, &pos))) {
1573 reg = NewRegion (region_name, pos);
1574 reg->type = region_type;
1575 }
1576 if (0 == reg->nobjs)
1577 reg->objid = (cgsize_t *) malloc ((size_t)region_nobjs * sizeof(cgsize_t));
1578 else
1579 reg->objid = (cgsize_t *) realloc (reg->objid,
1580 (size_t)(reg->nobjs + region_nobjs) * sizeof(cgsize_t));
1581 if (NULL == reg->objid)
1582 cgnsImportFatal (
1583 "cgnsImportRegion:malloc failed for the region object list");
1584
1585 for (n = 0; n < region_nobjs; n++)
1586 reg->objid[n+reg->nobjs] = region_list[n];
1587 reg->nobjs += region_nobjs;
1588 return (reg->nobjs);
1589 }
1590
1591 /*---------- cgnsImportRegion -------------------------------------------
1592 * import a named region
1593 *-----------------------------------------------------------------------*/
1594
cgnsImportRegion(char * regname,int regtype,cgsize_t numobjs,cgsize_t * objlist)1595 cgsize_t cgnsImportRegion (char *regname, int regtype,
1596 cgsize_t numobjs, cgsize_t *objlist)
1597 {
1598 cgnsImportBegReg (regname, regtype);
1599 cgnsImportAddReg (numobjs, objlist);
1600 return (cgnsImportEndReg ());
1601 }
1602
1603 /*---------- cgnsImportRegionList ---------------------------------------
1604 * return a list of all region names
1605 *-----------------------------------------------------------------------*/
1606
cgnsImportRegionList(void)1607 char **cgnsImportRegionList (void)
1608 {
1609 int n, len = 0;
1610 char **namelist, *names;
1611
1612 for (n = 0; n < num_regions; n++)
1613 len += ((int)strlen (reglist[n].name) + 1);
1614 n = num_regions + 1;
1615 namelist = (char **) malloc (len + n * sizeof(char *));
1616 if (NULL == namelist)
1617 cgnsImportFatal (
1618 "cgnsImportRegionList:malloc failed for region name list");
1619 names = (char *) (namelist + n);
1620 for (n = 0; n < num_regions; n++) {
1621 namelist[n] = names;
1622 strcpy (names, reglist[n].name);
1623 names += (strlen (reglist[n].name) + 1);
1624 }
1625 namelist[num_regions] = NULL;
1626 return (namelist);
1627 }
1628
1629 /*---------- cgnsImportGetRegion ----------------------------------------
1630 * get node ID's for a region
1631 *-----------------------------------------------------------------------*/
1632
cgnsImportGetRegion(char * regname)1633 cgsize_t *cgnsImportGetRegion (char *regname)
1634 {
1635 int pos;
1636 cgsize_t n, *objlist;
1637 cgnsREGN *reg;
1638
1639 if (NULL == regname || !*regname ||
1640 NULL == (reg = GetRegion (regname, &pos)))
1641 return (NULL);
1642 objlist = (cgsize_t *) malloc ((size_t)(reg->nobjs + 2) * sizeof(cgsize_t));
1643 if (NULL == objlist)
1644 cgnsImportFatal (
1645 "cgnsImportGetRegion:malloc failed for region object ID list");
1646 objlist[0] = reg->type;
1647 objlist[1] = reg->nobjs;
1648 for (n = 0; n < reg->nobjs; n++)
1649 objlist[n+2] = reg->objid[n];
1650 return (objlist);
1651 }
1652
1653 /*---------- cgnsImportOpen ---------------------------------------------
1654 * open CGNS file
1655 *-----------------------------------------------------------------------*/
1656
cgnsImportOpen(char * filename)1657 int cgnsImportOpen (char *filename)
1658 {
1659 cgnsImportClose ();
1660 if (cg_open (filename, CG_MODE_MODIFY, &cgnsFile) &&
1661 cg_open (filename, CG_MODE_WRITE, &cgnsFile))
1662 cgnsImportFatal ((char *)cg_get_error());
1663 return cgnsFile;
1664 }
1665
1666 /*---------- cgnsImportBase ---------------------------------------------
1667 * set CGNS base
1668 *-----------------------------------------------------------------------*/
1669
cgnsImportBase(char * basename)1670 int cgnsImportBase (char *basename)
1671 {
1672 if (cg_base_write (cgnsFile, basename, 3, 3, &cgnsBase))
1673 cgnsImportFatal ((char *)cg_get_error());
1674 return cgnsBase;
1675 }
1676
1677 /*---------- cgnsImportZone ---------------------------------------------
1678 * set CGNS zone
1679 *-----------------------------------------------------------------------*/
1680
cgnsImportZone(char * zonename)1681 void cgnsImportZone (char *zonename)
1682 {
1683 int n;
1684
1685 for (n = 0; n < num_nodes; n++) {
1686 if (nodemap[n].node != NULL)
1687 free (nodemap[n].node);
1688 }
1689 for (n = 0; n < num_elements; n++) {
1690 if (elemlist[n].nodeid != NULL)
1691 free (elemlist[n].nodeid);
1692 }
1693 if (num_regions) {
1694 for (n = 0; n < num_regions; n++) {
1695 if (reglist[n].objid != NULL)
1696 free (reglist[n].objid);
1697 }
1698 free (reglist);
1699 }
1700 if (num_faces) {
1701 for (n = 0; n < num_faces; n++) {
1702 if (facelist[n] != NULL) free (facelist[n]);
1703 }
1704 free (facelist);
1705 }
1706 num_nodes = num_elements = num_faces = 0;
1707 num_regions = 0;
1708
1709 strncpy (cgnsZoneName, zonename, 32);
1710 cgnsZoneName[32] = 0;
1711 }
1712
1713 /*---------- cgnsImportWrite --------------------------------------------
1714 * write data to the CGNS file
1715 *-----------------------------------------------------------------------*/
1716
cgnsImportWrite(void)1717 int cgnsImportWrite (void)
1718 {
1719 int icoord, isect;
1720 cgsize_t n, nn, nnodes, sizes[3];
1721 cgsize_t nc, nconn, *conns, *conns_offset, pos;
1722 CGNS_ENUMT(ElementType_t) elemtype = CGNS_ENUMV(ElementTypeNull);
1723 #ifdef DOUBLE_PRECISION
1724 CGNS_ENUMT(DataType_t) datatype = CGNS_ENUMV(RealDouble);
1725 double *xyz;
1726 #else
1727 CGNS_ENUMT(DataType_t) datatype = CGNS_ENUMV(RealSingle);
1728 float *xyz;
1729 #endif
1730 cgnsNODE *node;
1731 cgnsELEM *elem;
1732 cgnsREGN *regn;
1733
1734 if (!cgnsFile)
1735 cgnsImportFatal ("cgnsImportWrite:CGNS file not open");
1736 if (region_id)
1737 cgnsImportEndReg ();
1738
1739 /* count the nodes */
1740
1741 nnodes = 0;
1742 for (n = 0; n < num_nodes; n++) {
1743 if (nodemap[n].nodeid == nodemap[n].node->id &&
1744 USED_BIT == (nodemap[n].node->flags & USED_BIT))
1745 nnodes++;
1746 }
1747 if (!nnodes) return 0;
1748
1749 if (!cgnsBase) cgnsImportBase ("Base");
1750 if (!*cgnsZoneName) strcpy (cgnsZoneName, "Zone");
1751
1752 sizes[0] = nnodes;
1753 sizes[1] = num_elements;
1754 sizes[2] = 0;
1755
1756 if (cg_zone_write (cgnsFile, cgnsBase, cgnsZoneName,
1757 sizes, CGNS_ENUMV(Unstructured), &cgnsZone))
1758 cgnsImportFatal ((char *)cg_get_error());
1759
1760 /* write the node list */
1761
1762 #ifdef DOUBLE_PRECISION
1763 xyz = (double *) malloc ((size_t)nnodes * sizeof(double));
1764 #else
1765 xyz = (float *) malloc ((size_t)nnodes * sizeof(float));
1766 #endif
1767 if (NULL == xyz)
1768 cgnsImportFatal ("cgnsImportWrite:malloc failed for nodes");
1769
1770 for (nn = 0, n = 0; n < num_nodes; n++) {
1771 if (nodemap[n].nodeid == nodemap[n].node->id &&
1772 USED_BIT == (nodemap[n].node->flags & USED_BIT))
1773 #ifdef DOUBLE_PRECISION
1774 xyz[nn++] = (double)nodemap[n].node->x;
1775 #else
1776 xyz[nn++] = (float)nodemap[n].node->x;
1777 #endif
1778 }
1779 if (cg_coord_write (cgnsFile, cgnsBase, cgnsZone, datatype,
1780 "CoordinateX", (void *)xyz, &icoord))
1781 cgnsImportFatal ((char *)cg_get_error());
1782
1783 for (nn = 0, n = 0; n < num_nodes; n++) {
1784 if (nodemap[n].nodeid == nodemap[n].node->id &&
1785 USED_BIT == (nodemap[n].node->flags & USED_BIT))
1786 #ifdef DOUBLE_PRECISION
1787 xyz[nn++] = (double)nodemap[n].node->y;
1788 #else
1789 xyz[nn++] = (float)nodemap[n].node->y;
1790 #endif
1791 }
1792 if (cg_coord_write (cgnsFile, cgnsBase, cgnsZone, datatype,
1793 "CoordinateY", (void *)xyz, &icoord))
1794 cgnsImportFatal ((char *)cg_get_error());
1795
1796 for (nn = 0, n = 0; n < num_nodes; n++) {
1797 if (nodemap[n].nodeid == nodemap[n].node->id &&
1798 USED_BIT == (nodemap[n].node->flags & USED_BIT))
1799 #ifdef DOUBLE_PRECISION
1800 xyz[nn++] = (double)nodemap[n].node->z;
1801 #else
1802 xyz[nn++] = (float)nodemap[n].node->z;
1803 #endif
1804 }
1805 if (cg_coord_write (cgnsFile, cgnsBase, cgnsZone, datatype,
1806 "CoordinateZ", (void *)xyz, &icoord))
1807 cgnsImportFatal ((char *)cg_get_error());
1808
1809 /* write variables */
1810
1811 if (num_vars) {
1812 int isol, ifld, nv;
1813 if (cg_sol_write(cgnsFile, cgnsBase, cgnsZone,
1814 "NodeVariables", CGNS_ENUMV(Vertex), &isol))
1815 cgnsImportFatal ((char *)cg_get_error());
1816 for (nv = 0; nv < num_vars; nv++) {
1817 for (nn = 0, n = 0; n < num_nodes; n++) {
1818 if (nodemap[n].nodeid == nodemap[n].node->id &&
1819 USED_BIT == (nodemap[n].node->flags & USED_BIT))
1820 #ifdef DOUBLE_PRECISION
1821 xyz[nn++] = (double)nodemap[n].node->vars[nv];
1822 #else
1823 xyz[nn++] = (float)nodemap[n].node->vars[nv];
1824 #endif
1825 }
1826 if (strlen(var_names[nv]) > 32) var_names[nv][32] = 0;
1827 if (cg_field_write(cgnsFile, cgnsBase, cgnsZone, isol,
1828 datatype, var_names[nv], xyz, &ifld))
1829 cgnsImportFatal ((char *)cg_get_error());
1830 }
1831 }
1832
1833 free (xyz);
1834
1835 /* write the element list */
1836
1837 switch (elemlist->elemtype) {
1838 case cgnsELEM_TET:
1839 elemtype = CGNS_ENUMV(TETRA_4);
1840 break;
1841 case cgnsELEM_PYR:
1842 elemtype = CGNS_ENUMV(PYRA_5);
1843 break;
1844 case cgnsELEM_WDG:
1845 elemtype = CGNS_ENUMV(PENTA_6);
1846 break;
1847 case cgnsELEM_HEX:
1848 elemtype = CGNS_ENUMV(HEXA_8);
1849 break;
1850 }
1851 for (n = 0, elem = elemlist; n < num_elements; n++, elem++) {
1852 if (elem->elemtype != elemlist->elemtype) {
1853 elemtype = CGNS_ENUMV(MIXED);
1854 break;
1855 }
1856 }
1857
1858 if (elemtype == CGNS_ENUMV(MIXED)) {
1859 nconn = 0;
1860 for (n = 0, elem = elemlist; n < num_elements; n++, elem++)
1861 nconn += (1 + elem->elemtype);
1862 }
1863 else
1864 nconn = num_elements * elemlist->elemtype;
1865 conns = (cgsize_t *) malloc ((size_t)nconn * sizeof(cgsize_t));
1866 if (NULL == conns)
1867 cgnsImportFatal ("cgnsImportWrite:malloc failed for element data");
1868
1869 if (elemtype == CGNS_ENUMV(MIXED)) {
1870 conns_offset = (cgsize_t *) malloc ((size_t)(num_elements+1) * sizeof(cgsize_t));
1871 if (NULL == conns_offset)
1872 cgnsImportFatal ("cgnsImportWrite:malloc failed for element offset data");
1873 conns_offset[0] = 0;
1874 }
1875
1876 nc = 0;
1877 for (n = 0, elem = elemlist; n < num_elements; n++, elem++) {
1878 if (elemtype == CGNS_ENUMV(MIXED)) {
1879 switch (elem->elemtype) {
1880 case cgnsELEM_TET :
1881 conns[nc] = CGNS_ENUMV(TETRA_4);
1882 break;
1883 case cgnsELEM_PYR:
1884 conns[nc] = CGNS_ENUMV(PYRA_5);
1885 break;
1886 case cgnsELEM_WDG:
1887 conns[nc] = CGNS_ENUMV(PENTA_6);
1888 break;
1889 case cgnsELEM_HEX:
1890 conns[nc] = CGNS_ENUMV(HEXA_8);
1891 break;
1892 }
1893 conns_offset[n+1] = conns_offset[n] + conns[nc] + 1;
1894 nc++;
1895 }
1896 for (nn = 0; nn < elem->elemtype; nn++) {
1897 if (NULL == (node = GetNode (elem->nodeid[nn], &pos)))
1898 cgnsImportFatal ("cgnsImportWrite:missing element node");
1899 conns[nc++] = pos + 1;
1900 }
1901 }
1902
1903 if (elemtype == CGNS_ENUMV(MIXED)) {
1904 if (cg_poly_section_write (cgnsFile, cgnsBase, cgnsZone, "GridElements",
1905 elemtype, 1, num_elements, 0, conns, conns_offset, &isect))
1906 cgnsImportFatal ((char *)cg_get_error());
1907 free(conns_offset);
1908 }
1909 else {
1910 if (cg_section_write (cgnsFile, cgnsBase, cgnsZone, "GridElements",
1911 elemtype, 1, num_elements, 0, conns, &isect))
1912 cgnsImportFatal ((char *)cg_get_error());
1913 }
1914
1915 free (conns);
1916
1917 /* write the regions */
1918
1919 nn = num_elements + 1;
1920 for (n = 0, regn = reglist; n < num_regions; n++, regn++) {
1921 if (regn->type == cgnsREG_NODES)
1922 nn += write_node_region (regn, nn);
1923 else if (regn->type == cgnsREG_FACES)
1924 nn += write_face_region (regn, nn);
1925 else if (regn->type == cgnsREG_ELEMS)
1926 nn += write_elem_region (regn, nn);
1927 }
1928
1929 return cgnsZone;
1930 }
1931
1932 /*---------- cgnsImportClose --------------------------------------------
1933 * close the CGNS file
1934 *-----------------------------------------------------------------------*/
1935
cgnsImportClose(void)1936 void cgnsImportClose (void)
1937 {
1938 if (cgnsFile) {
1939 cg_close (cgnsFile);
1940 cgnsFile = 0;
1941 }
1942 }
1943
1944