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 (&reglist[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 (&reglist[hi]);
403     }
404 
405     while (1) {
406         *pos = (lo + hi) >> 1;
407         if (0 == (cmp = strcmp (name, reglist[*pos].name)))
408             return (&reglist[*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 (&reglist[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