1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <string.h>
4 #include <ctype.h>
5 #include <math.h>
6 
7 #include "cgnslib.h"
8 #include "calc.h"
9 #include "vecerr.h"
10 
11 #if CGNS_VERSION < 3000
12 # define Celsius Celcius
13 #endif
14 
15 #ifndef CG_MODE_READ
16 # define CG_MODE_READ   MODE_READ
17 # define CG_MODE_MODIFY MODE_MODIFY
18 #endif
19 
20 #ifdef READ_NODE
21 #include "cgns_header.h"
22 #include "cgns_io.h"
23 #endif
24 
25 #ifndef CGNSTYPES_H
26 # define cgsize_t  int
27 # define cglong_t  long
28 # define cgulong_t unsigned long
29 #endif
30 #ifndef CGNS_ENUMT
31 # define CGNS_ENUMT(e) e
32 # define CGNS_ENUMV(e) e
33 #endif
34 #ifndef CG_MAX_INT32
35 # define CG_MAX_INT32 0x7FFFFFFF
36 #endif
37 
38 int cgnsFile = 0;
39 
40 /*--- base data ---*/
41 
42 int NumBases = 0;
43 int cgnsBase;
44 char BaseName[33];
45 int CellDim, PhyDim;
46 int BaseClass, BaseUnits[5];
47 
48 /*--- zone data ---*/
49 
50 int NumZones = 0;
51 int cgnsZone;
52 char ZoneName[33];
53 int ZoneType;
54 int ZoneDims[6];
55 int ZoneClass, ZoneUnits[5];
56 int GridClass, GridUnits[5];
57 
58 /*--- solution data ---*/
59 
60 int NumSolns = 0;
61 int cgnsSoln;
62 char SolnName[33];
63 int SolnLocation;
64 int SolnDims[3], SolnRind[6];
65 int SolnClass, SolnUnits[5];
66 
67 /*--- field data ---*/
68 
69 int NumReference = 0;
70 Variable *reference;
71 
72 int NumCoordinates = 0;
73 Variable *coordinates;
74 
75 int NumVariables = 0;
76 Variable *variables;
77 
78 /*--- local variables ---*/
79 
80 static int cmdstrlen = 0;
81 static char *cmdstr;
82 static int VectorLen = 100;
83 
84 /* unit specifications */
85 
86 typedef struct {
87     char *name;
88     int type;
89     int value;
90 } UnitSpec;
91 
92 static UnitSpec unitspec[] = {
93     {"cel", 3, CGNS_ENUMV(Celsius}),
94     {"cen", 1, CGNS_ENUMV(Centimeter}),
95     {"cm",  1, CGNS_ENUMV(Centimeter}),
96     {"c",   3, CGNS_ENUMV(Celsius}),
97     {"d",   4, CGNS_ENUMV(Degree}),
98     {"fa",  3, CGNS_ENUMV(Fahrenheit}),
99     {"fo",  1, CGNS_ENUMV(Foot}),
100     {"ft",  1, CGNS_ENUMV(Foot}),
101     {"f",   3, CGNS_ENUMV(Fahrenheit}),
102     {"g",   0, CGNS_ENUMV(Gram}),
103     {"in",  1, CGNS_ENUMV(Inch}),
104     {"ke",  3, CGNS_ENUMV(Kelvin}),
105     {"ki",  0, CGNS_ENUMV(Kilogram}),
106     {"kg",  0, CGNS_ENUMV(Kilogram}),
107     {"k",   3, CGNS_ENUMV(Kelvin}),
108     {"lb",  0, CGNS_ENUMV(PoundMass}),
109     {"me",  1, CGNS_ENUMV(Meter}),
110     {"mi",  1, CGNS_ENUMV(Millimeter}),
111     {"mm",  1, CGNS_ENUMV(Millimeter}),
112     {"m",   1, CGNS_ENUMV(Meter}),
113     {"p",   0, CGNS_ENUMV(PoundMass}),
114     {"rad", 4, CGNS_ENUMV(Radian}),
115     {"ran", 3, CGNS_ENUMV(Rankine}),
116     {"r",   3, CGNS_ENUMV(Rankine}),
117     {"se",  2, CGNS_ENUMV(Second}),
118     {"sl",  0, CGNS_ENUMV(Slug}),
119     {"s",   2, CGNS_ENUMV(Second})
120 };
121 
122 #define NUM_UNITSPEC (sizeof(unitspec)/sizeof(UnitSpec))
123 
124 /*---------- read_node ---------------------------------------------
125  * read a node from the CGNS file
126  *------------------------------------------------------------------*/
127 
128 #ifdef READ_NODE
129 
130 static VECDATA *read_node (char *nodename)
131 {
132     cgns_file *cgfile;
133     int n, bytes, dt, cgio;
134     int ndim;
135     cgsize_t np, dims[CGIO_MAX_DIMENSIONS];
136     char *values;
137     char type[CGIO_MAX_DATATYPE_LENGTH+1];
138     char errmsg[CGIO_MAX_ERROR_LENGTH+1];
139     double rootid, nodeid;
140     VECDATA *vd;
141 
142     static struct dataTypes {
143         char *name;
144         int bytes;
145     } data_types[6] = {
146         {"I4", 4},
147         {"I8", 8},
148         {"U4", 4},
149         {"U8", 8},
150         {"R4", 4},
151         {"R8", 8}
152     };
153 
154     /* get node ID for node */
155 
156     cgfile = cgi_get_file (cgnsFile);
157     cgio = cgfile->cgio;
158     rootid = cgfile->rootid;
159 
160     if (cgio_get_node_id (cgio, rootid, nodename, &nodeid)) {
161         cgio_error_message (errmsg);
162         cgnsCalcFatal (errmsg);
163     }
164 
165     /* get the type of data */
166 
167     if (cgio_get_data_type (cgio, nodeid, type)) {
168         cgio_error_message (errmsg);
169         cgnsCalcFatal (errmsg);
170     }
171     for (n = 0; n < CGIO_MAX_DATATYPE_LENGTH && type[n]; n++) {
172         if (islower (type[n]))
173             type[n] = toupper (type[n]);
174     }
175     for (bytes = 0, dt = 0; dt < 6; dt++) {
176         if (0 == strncmp (type, data_types[dt].name, 2)) {
177             bytes = data_types[dt].bytes;
178             break;
179         }
180     }
181     if (bytes == 0) {
182         sprintf (errmsg, "can't handle data type %s", type);
183         cgnsCalcFatal (errmsg);
184     }
185 
186     /* get data dimensions */
187 
188     if (cgio_get_dimensions (cgio, nodeid, &ndim, dims)) {
189         cgio_error_message (errmsg);
190         cgnsCalcFatal (errmsg);
191     }
192     np = 0;
193     if (ndim > 0) {
194         for (np = 1, n = 0; n < ndim; n++)
195             np *= dims[n];
196     }
197     if (np == 0)
198         cgnsCalcFatal ("no data for node");
199     if (np > CG_MAX_INT32)
200         cgnsCalcFatal ("exceeded 32-bit integer");
201 
202     /* read the data */
203 
204     values = (char *) malloc ((size_t)(np * bytes));
205     if (NULL == values)
206         cgnsCalcFatal ("malloc failed for node data");
207 
208     if (cgio_read_all_data_type (cgio, nodeid, type, values)) {
209         cgio_error_message (errmsg);
210         cgnsCalcFatal (errmsg);
211     }
212 
213     if (np == 1) {
214         vd = vec_create (VEC_VALUE, 0, 1);
215         if (dt == 0) {
216             int *data = (int *)values;
217             vd->f.val = (VECFLOAT)*data;
218         }
219         else if (dt == 1) {
220             cglong_t *data = (cglong_t *)values;
221             vd->f.val = (VECFLOAT)*data;
222         }
223         else if (dt == 2) {
224             unsigned int *data = (unsigned int *)values;
225             vd->f.val = (VECFLOAT)*data;
226         }
227         else if (dt == 3) {
228             cgulong_t *data = (cgulong_t *)values;
229             vd->f.val = (VECFLOAT)*data;
230         }
231         else if (dt == 4) {
232             float *data = (float *)values;
233             vd->f.val = (VECFLOAT)*data;
234         }
235         else {
236             double *data = (double *)values;
237             vd->f.val = (VECFLOAT)*data;
238         }
239     }
240     else {
241         vd = vec_create (VEC_VECTOR, (int)np, 1);
242         if (dt == 0) {
243             int *data = (int *)values;
244             for (n = 0; n < np; n++)
245                 vd->f.vec[n] = (VECFLOAT)data[n];
246         }
247         else if (dt == 1) {
248             cglong_t *data = (cglong_t *)values;
249             for (n = 0; n < np; n++)
250                 vd->f.vec[n] = (VECFLOAT)data[n];
251         }
252         else if (dt == 2) {
253             unsigned int *data = (unsigned int *)values;
254             for (n = 0; n < np; n++)
255                 vd->f.vec[n] = (VECFLOAT)data[n];
256         }
257         else if (dt == 3) {
258             cgulong_t *data = (cgulong_t *)values;
259             for (n = 0; n < np; n++)
260                 vd->f.vec[n] = (VECFLOAT)data[n];
261         }
262         else if (dt == 4) {
263             float *data = (float *)values;
264             for (n = 0; n < np; n++)
265                 vd->f.vec[n] = (VECFLOAT)data[n];
266         }
267         else {
268             double *data = (double *)values;
269             for (n = 0; n < np; n++)
270                 vd->f.vec[n] = (VECFLOAT)data[n];
271         }
272     }
273 
274     free (values);
275     return vd;
276 }
277 
278 #endif
279 
280 /*---------- read_units -----------------------------------------------
281  * read unit specifications
282  *---------------------------------------------------------------------*/
283 
284 static int read_units (int units[5])
285 {
286     int n;
287     CGNS_ENUMT(MassUnits_t) mass;
288     CGNS_ENUMT(LengthUnits_t) length;
289     CGNS_ENUMT(TimeUnits_t) time;
290     CGNS_ENUMT(TemperatureUnits_t) temp;
291     CGNS_ENUMT(AngleUnits_t) angle;
292 
293     if (cg_units_read (&mass, &length, &time, &temp, &angle)) {
294         for (n = 0; n < 5; n++)
295             units[n] = 0;
296         return 0;
297     }
298     units[0] = mass;
299     units[1] = length;
300     units[2] = time;
301     units[3] = temp;
302     units[4] = angle;
303     return 1;
304 }
305 
306 /*---------- read_class -----------------------------------------------
307  * get data class, units and conversion factors
308  *---------------------------------------------------------------------*/
309 
310 static void read_class (Variable *var, int dataclass, int units[5])
311 {
312     int i;
313     CGNS_ENUMT(DataType_t) datatype;
314 
315     if (cg_dataclass_read ((CGNS_ENUMT(DataClass_t) *)&var->dataclass))
316         var->dataclass = dataclass;
317     var->hasunits = read_units (var->units);
318     if (!var->hasunits) {
319         for (i = 0; i < 5; i++)
320             var->units[i] = units[i];
321     }
322     if (cg_conversion_info (&datatype) ||
323        (datatype != CGNS_ENUMV(RealSingle) && datatype != CGNS_ENUMV(RealDouble))) {
324         var->hasconv = 0;
325         var->dataconv[0] = 1.0;
326         var->dataconv[1] = 0.0;
327     }
328     else {
329         var->hasconv = datatype;
330         if (datatype == CGNS_ENUMV(RealSingle)) {
331             float conv[2];
332             if (cg_conversion_read (conv))
333                 cgnsCalcFatal ((char *)cg_get_error());
334             for (i = 0; i < 2; i++)
335                 var->dataconv[i] = conv[i];
336         }
337         else {
338             if (cg_conversion_read (var->dataconv))
339                 cgnsCalcFatal ((char *)cg_get_error());
340         }
341     }
342     if (cg_exponents_info (&datatype) ||
343        (datatype != CGNS_ENUMV(RealSingle) && datatype != CGNS_ENUMV(RealDouble))) {
344         var->hasexp = 0;
345         for (i = 0; i < 5; i++)
346             var->exponent[i] = 0.0;
347     }
348     else {
349         var->hasexp = datatype;
350         if (datatype == CGNS_ENUMV(RealSingle)) {
351             float exp[5];
352             if (cg_exponents_read (exp))
353                 cgnsCalcFatal ((char *)cg_get_error());
354             for (i = 0; i < 5; i++)
355                 var->exponent[i] = exp[i];
356         }
357         else {
358             if (cg_exponents_read (var->exponent))
359                 cgnsCalcFatal ((char *)cg_get_error());
360         }
361     }
362 }
363 
364 /*---------- read_reference --------------------------------------------
365  * read reference conditions
366  *----------------------------------------------------------------------*/
367 
368 static void read_reference (void)
369 {
370     int n, narrays, na, dim;
371     cgsize_t vec[12];
372     CGNS_ENUMT(DataType_t) datatype;
373     char name[33];
374 
375     NumReference = 0;
376     if (cg_goto (cgnsFile, cgnsBase, "ReferenceState_t", 1, "end") ||
377         cg_narrays (&narrays) || narrays < 1)
378         return;
379     for (na = 1; na <= narrays; na++) {
380         if (cg_array_info (na, name, &datatype, &dim, vec))
381             cgnsCalcFatal ((char *)cg_get_error());
382         if (datatype != CGNS_ENUMV(Character) && dim >= 1 && vec[0] >= 1)
383             NumReference++;
384     }
385     if (!NumReference) return;
386     reference = (Variable *) malloc (NumReference * sizeof(Variable));
387     if (NULL == reference)
388         cgnsCalcFatal ("malloc failed for reference variables");
389 
390     for (n = 0, na = 1; na <= narrays; na++) {
391         if (cg_array_info (na, name, &datatype, &dim, vec))
392             cgnsCalcFatal ((char *)cg_get_error());
393         if (datatype != CGNS_ENUMV(Character) && dim >= 1 && vec[0] >= 1) {
394             dim *= vec[0];
395             strcpy (reference[n].name, name);
396             reference[n].type = 0;
397             reference[n].id = na;
398             reference[n].len = dim;
399             reference[n].valid = 1;
400             reference[n].datatype = datatype;
401             if (dim == 1) {
402                 reference[n].vd = vec_create (VEC_VALUE, 0, 0);
403                 if (cg_array_read_as (na, CGNS_ENUMV(RealDouble), &reference[n].vd->f.val))
404                     cgnsCalcFatal ((char *)cg_get_error());
405             }
406             else {
407                 reference[n].vd = vec_create (VEC_VECTOR, dim, 0);
408                 if (cg_array_read_as (na, CGNS_ENUMV(RealDouble), reference[n].vd->f.vec))
409                     cgnsCalcFatal ((char *)cg_get_error());
410             }
411             if (cg_goto (cgnsFile, cgnsBase, "ReferenceState_t", 1,
412                     "DataArray_t", na, "end"))
413                 cgnsCalcFatal ((char *)cg_get_error());
414             read_class (&reference[n], BaseClass, BaseUnits);
415             cg_goto (cgnsFile, cgnsBase, "ReferenceState_t", 1, "end");
416             n++;
417         }
418     }
419 }
420 
421 /*---------- get_variable ----------------------------------------------
422  * return variable data - read if necessary
423  *----------------------------------------------------------------------*/
424 
425 static VECDATA *get_variable (Variable *var)
426 {
427     if (var->vd == NULL) {
428         int n;
429         cgsize_t min[3], max[3];
430         for (n = 0; n < 3; n++) {
431             min[n] = 1;
432             max[n] = SolnDims[n];
433         }
434         var->vd = vec_create (VEC_VECTOR, var->len, 0);
435         if (cg_field_read (cgnsFile, cgnsBase, cgnsZone, cgnsSoln,
436             var->name, CGNS_ENUMV(RealDouble), min, max, var->vd->f.vec))
437             cgnsCalcFatal ((char *)cg_get_error());
438     }
439     return var->vd;
440 }
441 
442 /*---------- get_coordinate --------------------------------------------
443  * return coordinate data - read if necessary
444  *----------------------------------------------------------------------*/
445 
446 static VECDATA *get_coordinate (Variable *var)
447 {
448     if (var->vd == NULL) {
449         int n;
450         cgsize_t min[3], max[3];
451         for (n = 0; n < 3; n++) {
452             min[n] = 1;
453             max[n] = ZoneDims[n];
454         }
455         var->vd = vec_create (VEC_VECTOR, var->len, 0);
456         if (cg_coord_read (cgnsFile, cgnsBase, cgnsZone,
457             var->name, CGNS_ENUMV(RealDouble), min, max, var->vd->f.vec))
458             cgnsCalcFatal ((char *)cg_get_error());
459     }
460     return var->vd;
461 }
462 
463 /*---------- print_error ---------------------------------------------
464  * print error message on parsing error
465  *--------------------------------------------------------------------*/
466 
467 static void print_error (int errnum, char *errmsg, int pos, char *str)
468 {
469     printf (errnum < 0 ? "FATAL:" : "ERROR:");
470     if (NULL != str) {
471         printf ("%s\n      ", str);
472         while (pos-- > 0)
473             putchar ('-');
474         putchar ('^');
475     }
476     printf ("%s\n", errmsg);
477 }
478 
479 /*---------- get_name -----------------------------------------------
480  * get symbol name from string
481  *-------------------------------------------------------------------*/
482 
483 static char *get_name (char **str)
484 {
485     int n;
486     char *p = *str;
487     static char name[SYMNAME_MAXLEN+1];
488 
489     if (*p == '"') {
490         for (++p, n = 0; n < SYMNAME_MAXLEN && *p; n++) {
491             if (*p == '"') break;
492             name[n] = *p++;
493         }
494         if (*p++ != '"') return NULL;
495     }
496     else if (!isalpha (*p) && *p != '_')
497         return (NULL);
498     else {
499         for (n = 0; n < SYMNAME_MAXLEN && *p; n++) {
500             if (!isalnum(*p) && *p != '_') break;
501             name[n] = *p++;
502         }
503     }
504     name[n] = 0;
505     *str = p;
506     return name;
507 }
508 
509 /*---------- print_symbols ------------------------------------------
510  * print symbol names
511  *-------------------------------------------------------------------*/
512 
513 static int print_symbols (VECSYM *sym, void *data)
514 {
515     FILE *fp = (FILE *)data;
516 
517     if (VECSYM_EQUSTR == vecsym_type(sym))
518         fprintf (fp, "%s{%d}\n", vecsym_name(sym), vecsym_nargs(sym));
519     else if (VECSYM_MACRO  == vecsym_type(sym))
520         fprintf (fp, "%s<%d>\n", vecsym_name(sym), vecsym_nargs(sym));
521     else if (VECSYM_FUNC == vecsym_type(sym)) {
522         if (vecsym_nargs(sym) < 0)
523             fprintf (fp, "%s(...)\n", vecsym_name(sym));
524         else
525             fprintf (fp, "%s(%d)\n", vecsym_name(sym), vecsym_nargs(sym));
526     }
527     else if (VECSYM_VECTOR == vecsym_type(sym))
528         fprintf (fp, "%s[%ld]\n", vecsym_name(sym), (long)vecsym_veclen(sym));
529     else
530         fprintf (fp, "%s\n", vecsym_name(sym));
531     return 0;
532 }
533 
534 /*---------- delete_units -------------------------------------------
535  * called when symbol deleted
536  *-------------------------------------------------------------------*/
537 
538 static void delete_units (VECSYM *sym)
539 {
540     if (vecsym_user(sym) != NULL)
541         free (vecsym_user(sym));
542 }
543 
544 /*---------- callback -----------------------------------------------
545  * callback function for vector parser
546  *-------------------------------------------------------------------*/
547 
548 static VECDATA *callback (int check, char **pp, char **err)
549 {
550     int n, type = 0;
551     char *p = *pp, *name;
552 
553     /* check for reading node data */
554 
555 #ifdef READ_NODE
556     if (*p == '{') {
557         char nodename[257];
558         for (n = 0; n < 256; n++) {
559             if (!*++p || *p == '}') break;
560             nodename[n] = *p;
561         }
562         if (n == 256) {
563             *err = "internal node name length exceeded";
564             return NULL;
565         }
566         if (!n || *p != '}') {
567             *err = "incomplete node name specification";
568             return NULL;
569         }
570         *pp = ++p;
571         return read_node (nodename);
572     }
573 #endif
574 
575     /* check for reference or coordinate data */
576 
577     if (*p == '~' || *p == '\'')
578         type = *p++;
579 
580     /* get name */
581 
582     if ((name = get_name (&p)) != NULL) {
583 
584         /* check for variable */
585 
586         if (!type) {
587             for (n = 0; n < NumVariables; n++) {
588                 if (0 == strcmp (variables[n].name, name)) {
589                     *pp = p;
590                     return get_variable (&variables[n]);
591                 }
592             }
593         }
594 
595         /* check for grid coordinates */
596 
597         if (type != '~') {
598             for (n = 0; n < NumCoordinates; n++) {
599                 if (0 == strcmp (coordinates[n].name, name)) {
600                     *pp = p;
601                     return get_coordinate (&coordinates[n]);
602                 }
603             }
604         }
605 
606         /* check for reference quantity */
607 
608         if (type != '\'') {
609             for (n = 0; n < NumReference; n++) {
610                 if (0 == strcmp (reference[n].name, name)) {
611                     *pp = p;
612                     return reference[n].vd;
613                 }
614             }
615         }
616     }
617 
618     return (NULL);
619 }
620 
621 /*---------- parse_units -------------------------------------------
622  * get unit specification
623  *------------------------------------------------------------------*/
624 
625 static Units *parse_units (char **pp)
626 {
627     int n, par, div;
628     char *p = *pp, name[33];
629     float exp;
630     Units units, *u;
631     UnitSpec *us;
632 
633     for (n = 0; n < 5; n++) {
634         units.units[n] = 0;
635         units.exps[n] = 0.0;
636     }
637     par = div = 0;
638     while (1) {
639         while (*p && isspace (*p))
640             p++;
641         if (*p == '*' || *p == '-') {
642             p++;
643             continue;
644         }
645         if (*p == '/') {
646             div++;
647             p++;
648             continue;
649         }
650         if (*p == '(') {
651             par++;
652             p++;
653             continue;
654         }
655         if (*p == ')') {
656             if (!par--) return NULL;
657             if (div) div--;
658             p++;
659             continue;
660         }
661         n = 0;
662         while (*p && isalpha (*p)) {
663             if (n < 32) name[n++] = tolower (*p);
664             p++;
665         }
666         if (!n) break;
667         name[n] = 0;
668         for (n = 0; n < NUM_UNITSPEC; n++) {
669             if (name[0] == unitspec[n].name[0]) break;
670         }
671         us = NULL;
672         while (n < NUM_UNITSPEC) {
673             if (name[0] != unitspec[n].name[0]) break;
674             if (!strncmp (name, unitspec[n].name, strlen(unitspec[n].name))) {
675                 us = &unitspec[n];
676                 break;
677             }
678             n++;
679         }
680         if (us == NULL) return NULL;
681         while (*p && isspace (*p))
682             p++;
683         if (*p == '^') {
684             if (1 != sscanf (++p, "%f%n", &exp, &n)) return NULL;
685             for (p += n; *p && isspace(*p); p++)
686                 ;
687         }
688         else
689             exp = 1.0;
690         units.units[us->type] = us->value;
691         if (div) {
692             units.exps[us->type] -= exp;
693             if (*p != '-' && div > par) div--;
694         }
695         else
696             units.exps[us->type] += exp;
697     }
698     u = (Units *) malloc (sizeof(Units));
699     if (u == NULL)
700         cgnsCalcFatal ("malloc failed for unit specification");
701     for (n = 0; n < 5; n++) {
702         u->units[n] = units.units[n];
703         u->exps[n] = units.exps[n];
704     }
705     *pp = p;
706     return u;
707 }
708 
709 /*---------- free_all ----------------------------------------------
710  * free all data
711  *------------------------------------------------------------------*/
712 
713 static void free_all (void)
714 {
715     int n;
716 
717     if (NumReference) {
718         for (n = 0; n < NumReference; n++)
719             vec_destroy (reference[n].vd);
720         free (reference);
721         NumReference = 0;
722     }
723     if (NumCoordinates) {
724         for (n = 0; n < NumCoordinates; n++)
725             vec_destroy (coordinates[n].vd);
726         free (coordinates);
727         NumCoordinates = 0;
728     }
729     if (NumVariables) {
730         for (n = 0; n < NumVariables; n++)
731             vec_destroy (variables[n].vd);
732         free (variables);
733         NumVariables = 0;
734     }
735 }
736 
737 /*---------- cgnsCalcFatal -----------------------------------------
738  * terminate with error message
739  *------------------------------------------------------------------*/
740 
741 void cgnsCalcFatal (char *errmsg)
742 {
743     cgnsCalcDone ();
744     if (NULL != errmsg && *errmsg) {
745         if (NULL == vec_errhandler)
746             print_error (-1, errmsg, 0, NULL);
747         else
748             (*vec_errhandler) (-1, errmsg, 0, NULL);
749     }
750     exit (-1);
751 }
752 
753 /*---------- cgnsCalcError -----------------------------------------
754  * print error message
755  *------------------------------------------------------------------*/
756 
757 void cgnsCalcError (char *errmsg)
758 {
759     if (NULL != errmsg && *errmsg) {
760         if (NULL == vec_errhandler)
761             print_error (0, errmsg, 0, NULL);
762         else
763             (*vec_errhandler) (0, errmsg, 0, NULL);
764     }
765 }
766 
767 /*---------- cgnsCalcReset -----------------------------------------
768  * reset calculator (symbol table)
769  *------------------------------------------------------------------*/
770 
771 void cgnsCalcReset (void)
772 {
773     sym_free ();
774 #ifdef EXTERN_FUNCS
775     add_funcs ();
776 #endif
777     VectorLen = 100;
778 }
779 
780 /*---------- cgnsCalcInit ------------------------------------------
781  * load CGNS file and initialize
782  *------------------------------------------------------------------*/
783 
784 int cgnsCalcInit (char *cgnsfile, int modify,
785                   void (*errhandler)(int,char *,int,char *))
786 {
787     cgnsCalcDone ();
788 
789     /* set up error handler */
790 
791     if (NULL == errhandler)
792         vec_errhandler = print_error;
793     else
794         vec_errhandler = errhandler;
795 
796     /* callback to delete unit data */
797 
798     sym_delfunc = delete_units;
799 
800     /* open CGNS file */
801 
802     if (modify) {
803         if (cg_open (cgnsfile, CG_MODE_MODIFY, &cgnsFile))
804             cgnsCalcFatal ("couldn't open file in modify mode");
805     }
806     else {
807         if (cg_open (cgnsfile, CG_MODE_READ, &cgnsFile))
808             cgnsCalcFatal ("couldn't open file in read mode");
809     }
810     if (cg_nbases (cgnsFile, &NumBases))
811         cgnsCalcFatal ((char *)cg_get_error());
812     if (NumBases < 1)
813         cgnsCalcFatal ("no bases found in CGNS file");
814 
815     cgnsCalcBase (1);
816 
817     /* return number of bases */
818 
819     return NumBases;
820 }
821 
822 /*---------- cgnsCalcDone ------------------------------------------
823  * close CGNS file
824  *------------------------------------------------------------------*/
825 
826 void cgnsCalcDone (void)
827 {
828     if (cgnsFile) {
829         cg_close (cgnsFile);
830         cgnsFile = 0;
831     }
832     free_all ();
833     cgnsBase = NumBases = 0;
834     cgnsZone = NumZones = 0;
835     cgnsSoln = NumSolns = 0;
836 }
837 
838 /*---------- cgnsCalcBase ------------------------------------------
839  * set base for calculations
840  *------------------------------------------------------------------*/
841 
842 int cgnsCalcBase (int base)
843 {
844     if (base < 1 || base > NumBases)
845         cgnsCalcFatal ("invalid base specified");
846     if (cg_base_read (cgnsFile, base, BaseName, &CellDim, &PhyDim))
847         cgnsCalcFatal ((char *)cg_get_error());
848     free_all ();
849     cgnsBase = base;
850     cgnsZone = NumZones = 0;
851     cgnsSoln = NumSolns = 0;
852 
853     /* read base class and units */
854 
855     if (cg_goto (cgnsFile, cgnsBase, "end"))
856         cgnsCalcFatal ((char *)cg_get_error());
857     if (cg_dataclass_read ((CGNS_ENUMT(DataClass_t) *)&BaseClass))
858         BaseClass = 0;
859     read_units (BaseUnits);
860 
861     /* read reference conditions */
862 
863     read_reference ();
864 
865     /* get number of zones and initialize */
866 
867     if (cg_nzones (cgnsFile, cgnsBase, &NumZones))
868         cgnsCalcFatal ((char *)cg_get_error());
869     if (NumZones) cgnsCalcZone (1);
870 
871     return NumZones;
872 }
873 
874 /*---------- cgnsCalcZone ------------------------------------------
875  * set zone for calculations
876  *------------------------------------------------------------------*/
877 
878 int cgnsCalcZone (int zone)
879 {
880     int n, nc, size = 1;
881     cgsize_t dims[9];
882     CGNS_ENUMT(DataType_t) datatype;
883     char name[33];
884 
885     if (zone < 1 || zone > NumZones)
886         cgnsCalcFatal ("invalid zone specified");
887     if (cg_zone_read (cgnsFile, cgnsBase, zone, ZoneName, dims) ||
888         cg_zone_type (cgnsFile, cgnsBase, zone, (CGNS_ENUMT(ZoneType_t) *)&ZoneType))
889         cgnsCalcFatal ((char *)cg_get_error());
890     cgnsZone = zone;
891     cgnsSoln = NumSolns = 0;
892 
893     for (n = 0; n < 6; n++)
894         ZoneDims[n] = 1;
895     if (ZoneType == CGNS_ENUMV(Structured)) {
896         for (n = 0; n < CellDim; n++) {
897             ZoneDims[n] = dims[n];
898             ZoneDims[n+3] = dims[n+CellDim];
899             size *= dims[n];
900         }
901     }
902     else if (ZoneType == CGNS_ENUMV(Unstructured)) {
903         ZoneDims[0] = dims[0];
904         ZoneDims[3] = dims[1];
905         size = dims[0];
906     }
907     else
908         cgnsCalcFatal ("invalid zone type");
909     VectorLen = size;
910 
911     /* free-up previous data */
912 
913     if (NumCoordinates) {
914         for (n = 0; n < NumCoordinates; n++)
915             vec_destroy (coordinates[n].vd);
916         free (coordinates);
917         NumCoordinates = 0;
918     }
919     if (NumVariables) {
920         for (n = 0; n < NumVariables; n++)
921             vec_destroy (variables[n].vd);
922         free (variables);
923         NumVariables = 0;
924     }
925 
926     /* read zone class and units */
927 
928     if (cg_goto (cgnsFile, cgnsBase, "Zone_t", cgnsZone, "end"))
929         cgnsCalcFatal ((char *)cg_get_error());
930     if (cg_dataclass_read ((CGNS_ENUMT(DataClass_t) *)&ZoneClass))
931         ZoneClass = BaseClass;
932     if (!read_units (ZoneUnits)) {
933         for (n = 0; n < 5; n++)
934             ZoneUnits[n] = BaseUnits[n];
935     }
936     GridClass = ZoneClass;
937     for (n = 0; n < 5; n++)
938         GridUnits[n] = ZoneUnits[n];
939 
940     /* get coordinate info */
941 
942     if (cg_ncoords (cgnsFile, cgnsBase, cgnsZone, &nc))
943         cgnsCalcFatal ((char *)cg_get_error());
944     if (nc > 0) {
945 
946         /* read grid class and units */
947 
948         if (cg_goto (cgnsFile, cgnsBase, "Zone_t", cgnsZone,
949             "GridCoordinates_t", 1, "end"))
950             cgnsCalcFatal ((char *)cg_get_error());
951         if (cg_dataclass_read ((CGNS_ENUMT(DataClass_t) *)&GridClass))
952             GridClass = ZoneClass;
953         if (!read_units (GridUnits)) {
954             for (n = 0; n < 5; n++)
955                 GridUnits[n] = ZoneUnits[n];
956         }
957 
958         /* read coordinates */
959 
960         NumCoordinates = nc;
961         coordinates = (Variable *) malloc (NumCoordinates * sizeof(Variable));
962         if (NULL == coordinates)
963             cgnsCalcFatal ("malloc failed for coordinate info");
964 
965         for (n = 0; n < NumCoordinates; n++) {
966             if (cg_coord_info (cgnsFile, cgnsBase, cgnsZone, n+1,
967                 &datatype, name))
968                 cgnsCalcFatal ((char *)cg_get_error());
969             strcpy (coordinates[n].name, name);
970             coordinates[n].type = 1;
971             coordinates[n].id = n + 1;
972             coordinates[n].len = size;
973             coordinates[n].valid = 1;
974             coordinates[n].datatype = datatype;
975             if (cg_goto (cgnsFile, cgnsBase, "Zone_t", cgnsZone,
976                 "GridCoordinates_t", 1, "DataArray_t", n+1, "end"))
977                 cgnsCalcFatal ((char *)cg_get_error());
978             read_class (&coordinates[n], GridClass, GridUnits);
979             coordinates[n].vd = NULL;
980         }
981     }
982 
983     /* get number of solutions and initialize */
984 
985     if (cg_nsols (cgnsFile, cgnsBase, cgnsZone, &NumSolns))
986         cgnsCalcFatal ((char *)cg_get_error());
987     if (NumSolns) cgnsCalcSoln (1);
988 
989     return NumSolns;
990 }
991 
992 /*---------- cgnsCalcSoln ------------------------------------------
993  * set solution for calculations
994  *------------------------------------------------------------------*/
995 
996 int cgnsCalcSoln (int soln)
997 {
998     int i, n, size, nflds;
999     CGNS_ENUMT(DataType_t) datatype;
1000     char name[33];
1001 
1002     if (soln < 1 || soln > NumSolns)
1003         cgnsCalcFatal ("invalid solution specified");
1004     if (cg_sol_info (cgnsFile, cgnsBase, cgnsZone, soln, SolnName,
1005         (CGNS_ENUMT(GridLocation_t) *)&SolnLocation))
1006         cgnsCalcFatal ((char *)cg_get_error());
1007     cgnsSoln = soln;
1008 
1009     for (n = 0; n < 3; n++)
1010         SolnDims[n] = 1;
1011     if (ZoneType == CGNS_ENUMV(Structured)) {
1012         size = 1;
1013         if (SolnLocation == CGNS_ENUMV(CellCenter)) {
1014             if (cg_goto (cgnsFile, cgnsBase, "Zone_t", cgnsZone,
1015                 "FlowSolution_t", cgnsSoln, "end"))
1016                 cgnsCalcFatal ((char *)cg_get_error());
1017             if (cg_rind_read (SolnRind)) {
1018                 for (n = 0; n < 6; n++)
1019                     SolnRind[n] = 0;
1020             }
1021             for (i = 0, n = 0; n < CellDim; n++, i += 2) {
1022                 SolnDims[n] = ZoneDims[n] - 1 + SolnRind[i] + SolnRind[i+1];
1023                 size *= SolnDims[n];
1024             }
1025         }
1026         else {
1027             for (n = 0; n < 6; n++)
1028                 SolnRind[n] = 0;
1029             for (n = 0; n < CellDim; n++) {
1030                 SolnDims[n] = ZoneDims[n];
1031                 size *= SolnDims[n];
1032             }
1033         }
1034     }
1035     else {
1036         for (n = 0; n < 6; n++)
1037             SolnRind[n] = 0;
1038         size = SolnLocation == CGNS_ENUMV(CellCenter) ? ZoneDims[3] : ZoneDims[0];
1039         SolnDims[0] = size;
1040     }
1041     VectorLen = size;
1042 
1043     /* free-up previous data */
1044 
1045     if (NumVariables) {
1046         for (n = 0; n < NumVariables; n++)
1047             vec_destroy (variables[n].vd);
1048         free (variables);
1049         NumVariables = 0;
1050     }
1051 
1052     /* read solution class and units */
1053 
1054     if (cg_goto (cgnsFile, cgnsBase, "Zone_t", cgnsZone,
1055         "FlowSolution_t", cgnsSoln, "end"))
1056         cgnsCalcFatal ((char *)cg_get_error());
1057     if (cg_dataclass_read ((CGNS_ENUMT(DataClass_t) *)&SolnClass))
1058         SolnClass = ZoneClass;
1059     if (!read_units (SolnUnits)) {
1060         for (n = 0; n < 5; n++)
1061             SolnUnits[n] = ZoneUnits[n];
1062     }
1063 
1064     /* get field info */
1065 
1066     if (cg_nfields (cgnsFile, cgnsBase, cgnsZone, cgnsSoln, &nflds))
1067         cgnsCalcFatal ((char *)cg_get_error());
1068     if (nflds > 0) {
1069         NumVariables = nflds;
1070         variables = (Variable *) malloc (NumVariables * sizeof(Variable));
1071         if (NULL == variables)
1072             cgnsCalcFatal ("malloc failed for field info");
1073         for (n = 0; n < NumVariables; n++) {
1074             if (cg_field_info (cgnsFile, cgnsBase, cgnsZone, cgnsSoln,
1075                     n+1, &datatype, name))
1076                 cgnsCalcFatal ((char *)cg_get_error());
1077             strcpy (variables[n].name, name);
1078             variables[n].type = 2;
1079             variables[n].id = n + 1;
1080             variables[n].len = size;
1081             variables[n].valid = 1;
1082             variables[n].datatype = datatype;
1083             if (cg_goto (cgnsFile, cgnsBase, "Zone_t", cgnsZone,
1084                 "FlowSolution_t", cgnsSoln, "DataArray_t", n+1, "end"))
1085                 cgnsCalcFatal ((char *)cg_get_error());
1086             read_class (&variables[n], SolnClass, SolnUnits);
1087             variables[n].vd = NULL;
1088         }
1089     }
1090 
1091     return NumVariables;
1092 }
1093 
1094 /*---------- cgnsCalcCheck -----------------------------------------
1095  * check expression
1096  *------------------------------------------------------------------*/
1097 
1098 int cgnsCalcCheck (char *expression)
1099 {
1100     int length = (int)strlen (expression);
1101     char *p, *cmd, *name;
1102     Units *units;
1103 
1104     if (length > cmdstrlen) {
1105         if (cmdstrlen == 0) {
1106             cmdstrlen = length > 256 ? length : 256;
1107             cmdstr = (char *) malloc (cmdstrlen + 1);
1108         }
1109         else {
1110             cmdstrlen = length + 32;
1111             cmdstr = (char *) realloc (cmdstr, cmdstrlen + 1);
1112         }
1113         if (cmdstr == NULL)
1114             cgnsCalcFatal ("malloc failed for cmdstr");
1115     }
1116     cmd = strcpy (cmdstr, expression);
1117 
1118     for (p = cmd + strlen(cmd) - 1; p >= cmd && isspace(*p); p--)
1119         ;
1120     *++p = 0;
1121     while (*cmd && isspace(*cmd))
1122         cmd++;
1123 
1124     if (!*cmd) return (0);
1125 
1126     p = cmd;
1127     if ((name = get_name (&p)) != NULL) {
1128         int nargs = -1;
1129         char *equ;
1130 
1131         for (equ = p; *equ && isspace(*equ); equ++)
1132             ;
1133 
1134         /* check for units */
1135 
1136         if ('[' == *equ) {
1137             equ++;
1138             units = parse_units (&equ);
1139             if (units == NULL || *equ != ']') {
1140                 cgnsCalcError ("bad units specification");
1141                 if (units != NULL) free (units);
1142                 return 0;
1143             }
1144             free (units);
1145             while (*++equ && isspace (*equ))
1146                 ;
1147             if (!*equ) {
1148                 if (find_symbol (name, 0) != NULL) return 1;
1149             }
1150         }
1151 
1152         /* check for equation */
1153 
1154         if ('(' == *equ) {
1155             char *arg = equ;
1156             while (*++arg && (isspace (*arg) || isdigit(*arg)))
1157                 ;
1158             if (')' == *arg) {
1159                 nargs = atoi (equ + 1);
1160                 for (equ = arg+1; *equ && isspace (*equ); equ++)
1161                     ;
1162             }
1163         }
1164 
1165         if ('=' == *equ && '=' != *++equ) {
1166             for (cmd = equ; *cmd && isspace(*cmd); cmd++)
1167                 ;
1168             if (nargs > 9) {
1169                 cgnsCalcError ("invalid number of equation arguments");
1170                 return 0;
1171             }
1172         }
1173     }
1174 
1175     return vec_check (cmd, VectorLen, callback);
1176 }
1177 
1178 /*---------- cgnsCalcCommand ---------------------------------------
1179  * parse expression and return results
1180  *------------------------------------------------------------------*/
1181 
1182 VECSYM *cgnsCalcCommand (char *expression)
1183 {
1184     int n, length = (int)strlen (expression);
1185     char *p, *cmd, *name, sym[SYMNAME_MAXLEN+1];
1186     VECDATA *vd;
1187     Units *units = NULL;
1188 
1189     if (length > cmdstrlen) {
1190         if (cmdstrlen == 0) {
1191             cmdstrlen = length > 256 ? length : 256;
1192             cmdstr = (char *) malloc (cmdstrlen + 1);
1193         }
1194         else {
1195             cmdstrlen = length + 32;
1196             cmdstr = (char *) realloc (cmdstr, cmdstrlen + 1);
1197         }
1198         if (cmdstr == NULL)
1199             cgnsCalcFatal ("malloc failed for cmdstr");
1200     }
1201     cmd = strcpy (cmdstr, expression);
1202 
1203     /* skip leading and trailing space */
1204 
1205     for (p = cmd + strlen(cmd) - 1; p >= cmd && isspace(*p); p--)
1206         ;
1207     *++p = 0;
1208     while (*cmd && isspace(*cmd))
1209         cmd++;
1210 
1211     /* empty string */
1212 
1213     if (!*cmd) return (NULL);
1214 
1215     /* check for defining a new symbol */
1216 
1217     p = cmd;
1218     strcpy (sym, "_temp_");
1219 
1220     if ((name = get_name (&p)) != NULL) {
1221         int nargs = -1;
1222         char *equ;
1223 
1224         for (equ = p; *equ && isspace(*equ); equ++)
1225             ;
1226 
1227         /* check for units */
1228 
1229         if ('[' == *equ) {
1230             equ++;
1231             units = parse_units (&equ);
1232             if (units == NULL || *equ != ']') {
1233                 cgnsCalcError ("bad units specification");
1234                 if (units != NULL) free (units);
1235                 return (NULL);
1236             }
1237             while (*++equ && isspace (*equ))
1238                 ;
1239             if (!*equ) {
1240                 VECSYM *symbol = find_symbol (name, 0);
1241                 if (symbol != NULL) {
1242                     if (vecsym_user(symbol) != NULL)
1243                         free (vecsym_user(symbol));
1244                     vecsym_user(symbol) = units;
1245                     return (symbol);
1246                 }
1247             }
1248         }
1249 
1250         /* check for equation */
1251 
1252         if ('(' == *equ) {
1253             char *arg = equ;
1254             while (*++arg && (isspace (*arg) || isdigit(*arg)))
1255                 ;
1256             if (')' == *arg) {
1257                 nargs = atoi (equ + 1);
1258                 for (equ = arg+1; *equ && isspace (*equ); equ++)
1259                     ;
1260             }
1261         }
1262 
1263         if ('=' == *equ && '=' != *++equ) {
1264             strcpy (sym, name);
1265             for (cmd = equ; *cmd && isspace(*cmd); cmd++)
1266                 ;
1267 
1268             /* add equation as string */
1269 
1270             if (nargs >= 0) {
1271                 if (nargs > 9) {
1272                     cgnsCalcError ("invalid number of equation arguments");
1273                     return (NULL);
1274                 }
1275                 n = sym_addequ (sym, nargs, cmd, units);
1276                 if (n) {
1277                     cgnsCalcError (sym_errmsg (n));
1278                     return (NULL);
1279                 }
1280                 return (find_symbol (sym, 1));
1281             }
1282         }
1283     }
1284 
1285     vd = vec_parse (cmd, VectorLen, callback);
1286 
1287     if (NULL == vd) {
1288         if (NULL != units) free (units);
1289         return (NULL);
1290     }
1291 
1292     /* add to symbol table */
1293 
1294     if (VEC_VALUE == vd->type)
1295         n = sym_addval (sym, vd->f.val, units);
1296     else
1297         n = sym_addvec (sym, vd->len, vd->f.vec, units);
1298     vec_destroy (vd);
1299     if (!n)
1300         return (find_symbol (sym, 0));
1301     cgnsCalcError (sym_errmsg (n));
1302     return (NULL);
1303 }
1304 
1305 /*---------- cgnsCalcVarGet ----------------------------------------
1306  * return a variable
1307  *------------------------------------------------------------------*/
1308 
1309 Variable *cgnsCalcVarGet (char *varname)
1310 {
1311     int n, type = 0;
1312     char *name = varname;
1313 
1314     if (name == NULL || !*name) return NULL;
1315     if (*name == '~' || *name == '\'')
1316         type = *name++;
1317 
1318     /* check for variable */
1319 
1320     if (!type) {
1321         for (n = 0; n < NumVariables; n++) {
1322             if (0 == strcmp (variables[n].name, name))
1323                 return &variables[n];
1324         }
1325     }
1326 
1327     /* check for grid coordinates */
1328 
1329     if (type != '~') {
1330         for (n = 0; n < NumCoordinates; n++) {
1331             if (0 == strcmp (coordinates[n].name, name))
1332                 return &coordinates[n];
1333         }
1334     }
1335 
1336     /* check for reference quantity */
1337 
1338     if (type != '\'') {
1339         for (n = 0; n < NumReference; n++) {
1340             if (0 == strcmp (reference[n].name, name))
1341                 return &reference[n];
1342         }
1343     }
1344     return NULL;
1345 }
1346 
1347 /*---------- cgnsCalcVarList ---------------------------------------
1348  * print variables
1349  *------------------------------------------------------------------*/
1350 
1351 void cgnsCalcVarList (FILE *fp)
1352 {
1353     int n;
1354 
1355     if (NULL == fp)
1356         fp = stdout;
1357 
1358     if (NumReference) {
1359         fprintf (fp, "=== Reference ===\n");
1360         for (n = 0; n < NumReference; n++) {
1361             if (reference[n].len > 1)
1362                 fprintf (fp, "~%s[%d]\n", reference[n].name, reference[n].len);
1363             else
1364                 fprintf (fp, "~%s\n", reference[n].name);
1365         }
1366     }
1367 
1368     if (NumCoordinates) {
1369         fprintf (fp, "=== Coordinates ===\n");
1370         for (n = 0; n < NumCoordinates; n++)
1371             fprintf (fp, "'%s[%d]\n", coordinates[n].name, coordinates[n].len);
1372     }
1373 
1374     if (NumVariables) {
1375         fprintf (fp, "=== Solution ===\n");
1376         for (n = 0; n < NumVariables; n++)
1377             fprintf (fp, "%s[%d]\n", variables[n].name, variables[n].len);
1378     }
1379 }
1380 
1381 /*---------- cgnsCalcSymList ---------------------------------------
1382  * print list of symbols
1383  *------------------------------------------------------------------*/
1384 
1385 void cgnsCalcSymList (FILE *fp)
1386 {
1387     if (NULL == fp)
1388         fp = stdout;
1389     fprintf (fp, "=== Symbols ===\n");
1390     sym_list (0, print_symbols, fp);
1391 }
1392 
1393