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