1 /*========================================================================= 2 3 Program: Visualization Toolkit 4 Module: $RCSfile: vtkPhastaReader.cxx,v $ 5 6 Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen 7 All rights reserved. 8 See Copyright.txt or http://www.kitware.com/Copyright.htm for details. 9 10 This software is distributed WITHOUT ANY WARRANTY; without even 11 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 12 PURPOSE. See the above copyright notice for more information. 13 14 =========================================================================*/ 15 #include "vtkPhastaReader.h" 16 17 #include "vtkByteSwap.h" 18 #include "vtkCellType.h" //added for constants such as VTK_TETRA etc... 19 #include "vtkDataArray.h" 20 #include "vtkIntArray.h" 21 #include "vtkDoubleArray.h" 22 #include "vtkFloatArray.h" 23 #include "vtkInformation.h" 24 #include "vtkInformationVector.h" 25 #include "vtkObjectFactory.h" 26 #include "vtkPointData.h" 27 #include "vtkCellData.h" 28 #include "vtkPointSet.h" 29 #include "vtkSmartPointer.h" 30 #include "vtkUnstructuredGrid.h" 31 //change 32 #include "vtkStreamingDemandDrivenPipeline.h" 33 #include "vtkPVXMLElement.h" 34 #include "vtkPVXMLParser.h" 35 36 #include "vtkCellData.h" 37 #include "vtkFieldData.h" 38 #include "vtkMultiBlockDataSet.h" 39 #include "vtkMultiPieceDataSet.h" 40 41 //change end 42 vtkCxxRevisionMacro(vtkPhastaReader, "$Revision: 1.11 $"); 43 vtkStandardNewMacro(vtkPhastaReader); 44 45 #include <vtkstd/map> 46 #include <vtkstd/vector> 47 #include <vtkstd/string> 48 #include <vtksys/ios/sstream> 49 50 //CHANGE//////////////////////////////////////////////////////////////// 51 #include "rdtsc.h" 52 #define clockRate 2670000000.0 53 unsigned long long start, end; 54 //double opentime_total = 0.0; 55 int LAST_FILE_ID; 56 57 void startTimer(unsigned long long* start) { 58 *start = rdtsc(); 59 } 60 61 void endTimer(unsigned long long* end) { 62 *end = rdtsc(); 63 } 64 65 void computeTime(unsigned long long* start, unsigned long long* end) { 66 double time = (double)((*end-*start)/clockRate); 67 opentime_total += time; 68 } 69 70 71 #define VERSION_INFO_HEADER_SIZE 8192 72 #define DB_HEADER_SIZE 1024 73 #define TWO_MEGABYTE 2097152 74 #define ENDIAN_TEST_NUMBER 12180 // Troy's Zip Code!! 75 #define MAX_PHASTA_FILES 64 76 #define MAX_PHASTA_FILE_NAME_LENGTH 1024 77 #define MAX_FIELDS_NUMBER 48 78 #define MAX_FIELDS_NAME_LENGTH 128 79 #define DefaultMHSize (4*1024*1024) 80 int MasterHeaderSize = DefaultMHSize; 81 int diff_endian = 0; 82 long long counter = 0; 83 84 enum PhastaIO_Errors 85 { 86 MAX_PHASTA_FILES_EXCEEDED = -1, 87 UNABLE_TO_OPEN_FILE = -2, 88 NOT_A_MPI_FILE = -3, 89 GPID_EXCEEDED = -4, 90 DATA_TYPE_ILLEGAL = -5, 91 }; 92 93 //CHANGE END//////////////////////////////////////////////////////////// 94 struct vtkPhastaReaderInternal 95 { 96 struct FieldInfo 97 { 98 int StartIndexInPhastaArray; 99 int NumberOfComponents; 100 int DataDependency; // 0-nodal, 1-elemental 101 vtkstd::string DataType; // "int" or "double" 102 vtkstd::string PhastaFieldTag; 103 104 FieldInfo() : StartIndexInPhastaArray(-1), NumberOfComponents(-1), DataDependency(-1), DataType(""), PhastaFieldTag("") 105 { 106 } 107 }; 108 109 typedef vtkstd::map<vtkstd::string, FieldInfo> FieldInfoMapType; 110 FieldInfoMapType FieldInfoMap; 111 }; 112 113 114 // Begin of copy from phastaIO 115 116 //CHANGE//////////////////////////////////////////////////////////////// 117 118 /***********************************************************************/ 119 /***************** NEW PHASTA IO CODE STARTS HERE **********************/ 120 /***********************************************************************/ 121 122 int partID_counter; 123 124 typedef struct 125 { 126 bool Wrong_Endian; /* default to false */ 127 128 char filename[MAX_PHASTA_FILE_NAME_LENGTH]; /* defafults to 1024 */ 129 int nppp; 130 int nPPF; 131 int nFiles; 132 int nFields; 133 unsigned long long my_offset; 134 135 char * master_header; 136 double * double_chunk; 137 int * int_chunk; 138 double * read_double_chunk; 139 int * read_int_chunk; 140 unsigned long long **my_offset_table; 141 unsigned long long **my_read_table; 142 int field_count; 143 int part_count; 144 int read_field_count; 145 int read_part_count; 146 int GPid; 147 int start_id; 148 unsigned long long next_start_address; 149 150 int myrank; 151 int numprocs; 152 int local_myrank; 153 int local_numprocs; 154 } phastaio_file_t; 155 156 //default: Paraview disabled 157 158 typedef struct 159 { 160 int fileID; 161 int nppf, nfields; 162 int GPid; 163 int read_field_count; 164 char * masterHeader; 165 unsigned long long **offset_table; 166 unsigned long long my_offset; 167 168 }serial_file; 169 170 serial_file *SerialFile; 171 phastaio_file_t *PhastaIOActiveFiles[MAX_PHASTA_FILES]; 172 int PhastaIONextActiveIndex = 0; /* indicates next index to allocate */ 173 174 //CHANGE END//////////////////////////////////////////////////////////// 175 176 #define swap_char(A,B) { ucTmp = A; A = B ; B = ucTmp; } 177 178 vtkstd::map< int , char* > LastHeaderKey; 179 vtkstd::vector< FILE* > fileArray; 180 vtkstd::vector< int > byte_order; 181 vtkstd::vector< int > header_type; 182 int DataSize=0; 183 int LastHeaderNotFound = 0; 184 int Wrong_Endian = 0 ; 185 int Strict_Error = 0 ; 186 int binary_format = 0; 187 188 // the caller has the responsibility to delete the returned string 189 char* vtkPhastaReader::StringStripper( const char istring[] ) 190 { 191 int length = strlen( istring ); 192 char* dest = new char [ length + 1 ]; 193 strcpy( dest, istring ); 194 dest[ length ] = '\0'; 195 196 if ( char* p = strpbrk( dest, " ") ) 197 { 198 *p = '\0'; 199 } 200 201 return dest; 202 } 203 204 int vtkPhastaReader::cscompare( const char teststring[], 205 const char targetstring[] ) 206 { 207 208 char* s1 = const_cast<char*>(teststring); 209 char* s2 = const_cast<char*>(targetstring); 210 211 while( *s1 == ' ') { s1++; } 212 while( *s2 == ' ') { s2++; } 213 while( ( *s1 ) 214 && ( *s2 ) 215 && ( *s2 != '?') 216 && ( tolower( *s1 )==tolower( *s2 ) ) ) 217 { 218 s1++; 219 s2++; 220 while( *s1 == ' ') { s1++; } 221 while( *s2 == ' ') { s2++; } 222 } 223 if ( !( *s1 ) || ( *s1 == '?') ) 224 { 225 return 1; 226 } 227 else 228 { 229 return 0; 230 } 231 } 232 233 void vtkPhastaReader::isBinary( const char iotype[] ) 234 { 235 char* fname = StringStripper( iotype ); 236 if ( cscompare( fname, "binary" ) ) 237 { 238 binary_format = 1; 239 } 240 else 241 { 242 binary_format = 0; 243 } 244 delete [] fname; 245 246 } 247 248 size_t vtkPhastaReader::typeSize( const char typestring[] ) 249 { 250 char* ts1 = StringStripper( typestring ); 251 252 if ( cscompare( "integer", ts1 ) ) 253 { 254 delete [] ts1; 255 return sizeof(int); 256 } 257 else if ( cscompare( "double", ts1 ) ) 258 { 259 delete [] ts1; 260 return sizeof( double ); 261 } 262 else if ( cscompare( "float", ts1 ) ) 263 { 264 delete [] ts1; 265 return sizeof( float ); 266 } 267 else 268 { 269 delete [] ts1; 270 fprintf(stderr,"unknown type : %s\n",ts1); 271 return 0; 272 } 273 } 274 275 int vtkPhastaReader::readHeader( FILE* fileObject, 276 const char phrase[], 277 int* params, 278 int expect ) 279 { 280 char* text_header; 281 char* token; 282 char Line[1024]; 283 char junk; 284 int FOUND = 0 ; 285 int real_length; 286 int skip_size, integer_value; 287 int rewind_count=0; 288 289 if( !fgets( Line, 1024, fileObject ) && feof( fileObject ) ) 290 { 291 rewind( fileObject ); 292 clearerr( fileObject ); 293 rewind_count++; 294 fgets( Line, 1024, fileObject ); 295 } 296 297 while( !FOUND && ( rewind_count < 2 ) ) 298 { 299 if ( ( Line[0] != '\n' ) && ( real_length = strcspn( Line, "#" )) ) 300 { 301 text_header = new char [ real_length + 1 ]; 302 strncpy( text_header, Line, real_length ); 303 text_header[ real_length ] =static_cast<char>(NULL); 304 token = strtok ( text_header, ":" ); 305 if( cscompare( phrase , token ) ) 306 { 307 FOUND = 1 ; 308 token = strtok( NULL, " ,;<>" ); 309 skip_size = atoi( token ); 310 int i; 311 for( i=0; i < expect && ( token = strtok( NULL," ,;<>") ); i++) 312 { 313 params[i] = atoi( token ); 314 } 315 if ( i < expect ) 316 { 317 fprintf(stderr,"Expected # of ints not found for: %s\n",phrase ); 318 } 319 } 320 else if ( cscompare(token,"byteorder magic number") ) 321 { 322 if ( binary_format ) 323 { 324 fread((void*)&integer_value,sizeof(int),1,fileObject); 325 fread( &junk, sizeof(char), 1 , fileObject ); 326 if ( 362436 != integer_value ) 327 { 328 Wrong_Endian = 1; 329 } 330 } 331 else 332 { 333 fscanf(fileObject, "%d\n", &integer_value ); 334 } 335 } 336 else 337 { 338 /* some other header, so just skip over */ 339 token = strtok( NULL, " ,;<>" ); 340 skip_size = atoi( token ); 341 if ( binary_format) 342 { 343 fseek( fileObject, skip_size, SEEK_CUR ); 344 } 345 else 346 { 347 for( int gama=0; gama < skip_size; gama++ ) 348 { 349 fgets( Line, 1024, fileObject ); 350 } 351 } 352 } 353 delete [] text_header; 354 } 355 356 if ( !FOUND ) 357 { 358 if( !fgets( Line, 1024, fileObject ) && feof( fileObject ) ) 359 { 360 rewind( fileObject ); 361 clearerr( fileObject ); 362 rewind_count++; 363 fgets( Line, 1024, fileObject ); 364 } 365 } 366 } 367 368 if ( !FOUND ) 369 { 370 fprintf(stderr, "Error: Cound not find: %s\n", phrase); 371 return 1; 372 } 373 return 0; 374 } 375 376 void vtkPhastaReader::SwapArrayByteOrder_( void* array, 377 int nbytes, 378 int nItems ) 379 { 380 /* This swaps the byte order for the array of nItems each 381 of size nbytes , This will be called only locally */ 382 int i,j; 383 unsigned char ucTmp; 384 unsigned char* ucDst = (unsigned char*)array; 385 386 for(i=0; i < nItems; i++) 387 { 388 for(j=0; j < (nbytes/2); j++) 389 { 390 swap_char( ucDst[j] , ucDst[(nbytes - 1) - j] ); 391 } 392 ucDst += nbytes; 393 } 394 } 395 396 //CHANGE/////////////////////////////////////////////////////// 397 398 void vtkPhastaReader::queryphmpiio_(const char filename[],int *nfields, int *nppf) 399 { 400 401 FILE * fileHandle; 402 char* fname = StringStripper( filename ); 403 404 fileHandle = fopen (fname,"rb"); 405 if (fileHandle == NULL ) { 406 printf("\n File %s doesn't exist! Please check!\n",fname); 407 exit(1); 408 } 409 else 410 { 411 SerialFile =(serial_file *)calloc( 1, sizeof( serial_file) ); 412 413 SerialFile->masterHeader = (char *)malloc(MasterHeaderSize); 414 fread(SerialFile->masterHeader,1,MasterHeaderSize,fileHandle); 415 416 char read_out_tag[MAX_FIELDS_NAME_LENGTH]; 417 char * token; 418 int magic_number; 419 memcpy( read_out_tag, 420 SerialFile->masterHeader, 421 MAX_FIELDS_NAME_LENGTH-1 ); 422 423 if ( cscompare ("MPI_IO_Tag",read_out_tag) ) 424 { 425 // Test endianess ... 426 memcpy ( &magic_number, 427 SerialFile->masterHeader+sizeof("MPI_IO_Tag :"), 428 sizeof(int) ); 429 430 if ( magic_number != ENDIAN_TEST_NUMBER ) 431 { 432 diff_endian = 1; 433 } 434 435 char version[MAX_FIELDS_NAME_LENGTH/4]; 436 int mhsize; 437 438 memcpy(version, 439 SerialFile->masterHeader + MAX_FIELDS_NAME_LENGTH/2, 440 MAX_FIELDS_NAME_LENGTH/4 - 1); //TODO: why -1? 441 442 if( cscompare ("version",version) ) 443 { 444 // if there is "version" tag in the file, then it is newer format 445 // read master header size from here, otherwise use default 446 // TODO: if version is "1", we know mhsize is at 3/4 place... 447 448 token = strtok(version, ":"); 449 token = strtok(NULL, " ,;<>" ); 450 int iversion = atoi(token); 451 452 if( iversion == 1) { 453 memcpy( &mhsize, 454 SerialFile->masterHeader + MAX_FIELDS_NAME_LENGTH/4*3 + sizeof("mhsize : ")-1, 455 sizeof(int)); 456 if ( diff_endian) 457 SwapArrayByteOrder_(&mhsize, sizeof(int), 1); 458 free(SerialFile->masterHeader); 459 SerialFile->masterHeader = (char *)malloc(mhsize); 460 fseek(fileHandle, 0, SEEK_SET); 461 fread(SerialFile->masterHeader,1,mhsize,fileHandle); 462 } 463 //TODO: check if this is a valid int?? 464 MasterHeaderSize = mhsize; 465 } 466 else { // else it's version 0's format w/o version tag, implicating MHSize=4M 467 MasterHeaderSize = DefaultMHSize; 468 //printf("-----> version = 0; mhsize = %d\n", MasterHeaderSize); 469 } 470 471 // END OF CHANGE FOR VERSION 472 // 473 memcpy( read_out_tag, 474 SerialFile->masterHeader+MAX_FIELDS_NAME_LENGTH+1, 475 MAX_FIELDS_NAME_LENGTH ); 476 477 // Read in # fields ... 478 token = strtok ( read_out_tag, ":" ); 479 token = strtok( NULL," ,;<>" ); 480 *nfields = atoi( token ); 481 SerialFile->nfields=*nfields; 482 483 memcpy( read_out_tag, 484 SerialFile->masterHeader+ 485 *nfields * MAX_FIELDS_NAME_LENGTH + 486 MAX_FIELDS_NAME_LENGTH * 2, 487 MAX_FIELDS_NAME_LENGTH); 488 489 token = strtok ( read_out_tag, ":" ); 490 token = strtok( NULL," ,;<>" ); 491 *nppf = atoi( token ); 492 SerialFile->nppf=*nppf; 493 } 494 else 495 { 496 printf("The file you opened is not new format, please check!\n"); 497 } 498 499 fclose(fileHandle); 500 } 501 delete [] fname; 502 503 } 504 505 void vtkPhastaReader::finalizephmpiio_( int *fileDescriptor ) 506 { 507 //printf("total open time is %lf\n", opentime_total); 508 // free master header, offset table [][], and serial file struc 509 free( SerialFile->masterHeader); 510 int j; 511 for ( j = 0; j < SerialFile->nfields; j++ ) 512 { 513 free( SerialFile->offset_table[j] ); 514 } 515 free( SerialFile->offset_table); 516 free( SerialFile ); 517 } 518 519 char* StrReverse(char* str) 520 { 521 char *temp, *ptr; 522 int len, i; 523 524 temp=str; 525 for(len=0; *temp !='\0';temp++, len++); 526 527 ptr=(char*)malloc(sizeof(char)*(len+1)); 528 529 for(i=len-1; i>=0; i--) 530 ptr[len-i-1]=str[i]; 531 532 ptr[len]='\0'; 533 return ptr; 534 } 535 536 537 //CHANGE END////////////////////////////////////////////////// 538 539 void vtkPhastaReader::openfile( const char filename[], 540 const char mode[], 541 int* fileDescriptor ) 542 //CHANGE//////////////////////////////////////////////////// 543 { 544 //printf("in open(): counter = %ld\n", counter++); 545 546 FILE* file=NULL ; 547 *fileDescriptor = 0; 548 char* fname = StringStripper( filename ); 549 char* imode = StringStripper( mode ); 550 551 int string_length = strlen( fname ); 552 char* buffer = (char*) malloc ( string_length+1 ); 553 strcpy ( buffer, fname ); 554 buffer[ string_length ] = '\0'; 555 556 char* tempbuf = StrReverse(buffer); 557 free(buffer); 558 buffer = tempbuf; 559 560 //printf("buffer is %s\n",buffer); 561 562 char* st2 = strtok ( buffer, "." ); 563 //st2 = strtok (NULL, "."); 564 565 //printf("st2 is %s\n",st2); 566 567 string_length = strlen(st2); 568 char* buffer2 = (char*)malloc(string_length+1); 569 strcpy(buffer2,st2); 570 buffer2[string_length]='\0'; 571 572 char* tempbuf2 = StrReverse(buffer2); 573 free(buffer2); 574 buffer2 = tempbuf2; 575 //printf("buffer2 is %s\n",buffer2); 576 577 SerialFile->fileID = atoi(buffer2); 578 if ( char* p = strpbrk(buffer, "@") ) 579 *p = '\0'; 580 581 startTimer(&start); 582 if ( cscompare( "read", imode ) ) file = fopen(fname, "rb" ); 583 else if( cscompare( "write", imode ) ) file = fopen(fname, "wb" ); 584 else if( cscompare( "append", imode ) ) file = fopen(fname, "ab" ); 585 endTimer(&end); 586 computeTime(&start, &end); 587 588 589 if ( !file ){ 590 fprintf(stderr,"unable to open file : %s\n",fname ) ; 591 } else { 592 fileArray.push_back( file ); 593 byte_order.push_back( false ); 594 header_type.push_back( sizeof(int) ); 595 *fileDescriptor = fileArray.size(); 596 } 597 598 //////////////////////////////////////////////// 599 //unsigned long long **header_table; 600 SerialFile->offset_table = ( unsigned long long ** )calloc(SerialFile->nfields, 601 sizeof(unsigned long long *)); 602 603 int j; 604 for ( j = 0; j < SerialFile->nfields; j++ ) 605 { 606 SerialFile->offset_table[j]=( unsigned long long * ) calloc( SerialFile->nppf , 607 sizeof( unsigned long long)); 608 } 609 610 // Read in the offset table ... 611 for ( j = 0; j < SerialFile->nfields; j++ ) 612 { 613 614 memcpy( SerialFile->offset_table[j], 615 SerialFile->masterHeader + 616 VERSION_INFO_HEADER_SIZE + 617 j * SerialFile->nppf * sizeof(unsigned long long), 618 SerialFile->nppf * sizeof(unsigned long long) ); 619 620 if(diff_endian) { 621 SwapArrayByteOrder_( SerialFile->offset_table[j], 622 sizeof(unsigned long long int), 623 SerialFile->nppf); 624 } 625 // Swap byte order if endianess is different ... 626 /*if ( PhastaIOActiveFiles[i]->Wrong_Endian ) 627 { 628 SwapArrayByteOrder_( PhastaIOActiveFiles[i]->my_read_table[j], 629 sizeof(long long int), 630 PhastaIOActiveFiles[i]->nppp ); 631 } 632 */ 633 } 634 635 //////////////////////////////////////////////// 636 delete [] fname; 637 delete [] imode; 638 //free(fname); 639 //free(imode); 640 free(buffer); 641 free(buffer2); 642 } 643 644 //CHANGE END//////////////////////////////////////////////// 645 void vtkPhastaReader::closefile( int* fileDescriptor, 646 const char mode[] ) 647 //CHANGE/////////////////////////////////////////////// 648 { 649 char* imode = StringStripper( mode ); 650 651 if( cscompare( "write", imode ) 652 || cscompare( "append", imode ) ) { 653 fflush( fileArray[ *fileDescriptor - 1 ] ); 654 } 655 656 fclose( fileArray[ *fileDescriptor - 1 ] ); 657 delete [] imode; 658 } 659 660 661 //CHANGE END/////////////////////////////////////////// 662 663 664 void vtkPhastaReader::readheader( int* fileDescriptor, 665 const char keyphrase[], 666 void* valueArray, 667 int* nItems, 668 const char datatype[], 669 const char iotype[] ) 670 //CHANGE//////////////////////////////////////////////////// 671 { 672 int filePtr = *fileDescriptor - 1; 673 FILE* fileObject; 674 int* valueListInt; 675 676 if ( *fileDescriptor < 1 || *fileDescriptor > (int)fileArray.size() ) { 677 fprintf(stderr,"No file associated with Descriptor %d\n",*fileDescriptor); 678 fprintf(stderr,"openfile_ function has to be called before \n") ; 679 fprintf(stderr,"acessing the file\n ") ; 680 fprintf(stderr,"fatal error: cannot continue, returning out of call\n"); 681 return; 682 } 683 684 LastHeaderKey[ filePtr ] = const_cast< char* >( keyphrase ); 685 LastHeaderNotFound = false; 686 687 fileObject = fileArray[ filePtr ] ; 688 Wrong_Endian = byte_order[ filePtr ]; 689 690 isBinary( iotype ); 691 typeSize( datatype ); //redundant call, just avoid a compiler warning. 692 693 // right now we are making the assumption that we will only write integers 694 // on the header line. 695 696 valueListInt = static_cast< int* >( valueArray ); 697 698 ///////////////////////////////////////////////////////// 699 int j; 700 bool FOUND = false ; 701 unsigned int skip_size; 702 char * token; 703 char readouttag[MAX_FIELDS_NUMBER][MAX_FIELDS_NAME_LENGTH]; 704 705 int string_length = strlen( keyphrase ); 706 char* buffer = (char*) malloc ( string_length+1 ); 707 strcpy ( buffer, keyphrase ); 708 buffer[ string_length ] = '\0'; 709 710 char* st2 = strtok ( buffer, "@" ); 711 st2 = strtok (NULL, "@"); 712 SerialFile->GPid = atoi(st2); 713 if ( char* p = strpbrk(buffer, "@") ) 714 *p = '\0'; 715 716 //printf("field is %s and nfields is %d\n",keyphrase,SerialFile->nfields); 717 718 for ( j = 0; j<SerialFile->nfields; j++ ) 719 { 720 memcpy( readouttag[j], 721 SerialFile->masterHeader + j*MAX_FIELDS_NAME_LENGTH+MAX_FIELDS_NAME_LENGTH*2+1, 722 MAX_FIELDS_NAME_LENGTH-1 ); 723 } 724 725 for ( j = 0; j<SerialFile->nfields; j++ ) 726 { 727 token = strtok ( readouttag[j], ":" ); 728 729 if ( cscompare( buffer, token ) ) 730 { 731 SerialFile->read_field_count = j; 732 FOUND = true; 733 break; 734 } 735 } 736 if (!FOUND) 737 { 738 printf("Not found %s \n",keyphrase); 739 return; 740 } 741 742 int read_part_count = SerialFile->GPid - ( SerialFile->fileID - 1 ) * SerialFile->nppf - 1; 743 SerialFile->my_offset = SerialFile->offset_table[SerialFile->read_field_count][read_part_count]; 744 745 746 //printf("GP id is %d and fileID is %d and nppf is %d; ",SerialFile->GPid,SerialFile->fileID,SerialFile->nppf); 747 //printf("read field count is %d and read part count is %d; ",SerialFile->read_field_count,read_part_count); 748 749 char read_out_header[MAX_FIELDS_NAME_LENGTH]; 750 fseek(fileObject, SerialFile->my_offset+1, SEEK_SET); 751 fread( read_out_header, 1, MAX_FIELDS_NAME_LENGTH-1, fileObject ); 752 753 754 token = strtok ( read_out_header, ":" ); 755 756 if( cscompare( keyphrase , token ) ) 757 { 758 FOUND = true ; 759 token = strtok( NULL, " ,;<>" ); 760 skip_size = atoi( token ); 761 for( j=0; j < *nItems && ( token = strtok( NULL," ,;<>") ); j++ ) 762 valueListInt[j] = atoi( token ); 763 //printf("$$Keyphrase is %s Value list [0] is %d \n",keyphrase,valueListInt[0] ); 764 if ( j < *nItems ) 765 { 766 fprintf( stderr, "Expected # of ints not found for: %s\n", keyphrase ); 767 } 768 } 769 770 ///////////////////////////////////////////////////////// 771 772 byte_order[ filePtr ] = Wrong_Endian ; 773 774 //if ( ierr ) LastHeaderNotFound = true; 775 776 free(buffer); 777 778 return; 779 } 780 781 //CHANGE END//////////////////////////////////////////////// 782 783 void vtkPhastaReader::readdatablock( int* fileDescriptor, 784 const char keyphrase[], 785 void* valueArray, 786 int* nItems, 787 const char datatype[], 788 const char iotype[] ) 789 //CHANGE////////////////////////////////////////////////////// 790 { 791 792 int filePtr = *fileDescriptor - 1; 793 FILE* fileObject; 794 char junk; 795 796 if ( *fileDescriptor < 1 || *fileDescriptor > (int)fileArray.size() ) { 797 fprintf(stderr,"No file associated with Descriptor %d\n",*fileDescriptor); 798 fprintf(stderr,"openfile_ function has to be called before \n") ; 799 fprintf(stderr,"acessing the file\n ") ; 800 fprintf(stderr,"fatal error: cannot continue, returning out of call\n"); 801 return; 802 } 803 804 // error check.. 805 // since we require that a consistant header always preceed the data block 806 // let us check to see that it is actually the case. 807 808 if ( ! cscompare( LastHeaderKey[ filePtr ], keyphrase ) ) { 809 fprintf(stderr, "Header not consistant with data block\n"); 810 fprintf(stderr, "Header: %s\n", LastHeaderKey[ filePtr ] ); 811 fprintf(stderr, "DataBlock: %s\n ", keyphrase ); 812 fprintf(stderr, "Please recheck read sequence \n"); 813 if( Strict_Error ) { 814 fprintf(stderr, "fatal error: cannot continue, returning out of call\n"); 815 return; 816 } 817 } 818 819 if ( LastHeaderNotFound ) return; 820 821 fileObject = fileArray[ filePtr ]; 822 Wrong_Endian = byte_order[ filePtr ]; 823 //printf("in readdatablock(): wrong_endian = %d\n", Wrong_Endian); 824 825 size_t type_size = typeSize( datatype ); 826 int nUnits = *nItems; 827 isBinary( iotype ); 828 829 if ( binary_format ) { 830 fseek(fileObject, SerialFile->my_offset+DB_HEADER_SIZE, SEEK_SET); 831 832 fread( valueArray, type_size, nUnits, fileObject ); 833 //fread( &junk, sizeof(char), 1 , fileObject ); 834 //if ( Wrong_Endian ) SwapArrayByteOrder_( valueArray, type_size, nUnits ); 835 if ( diff_endian ) 836 SwapArrayByteOrder_( valueArray, type_size, nUnits ); // fj 837 } else { 838 839 char* ts1 = StringStripper( datatype ); 840 if ( cscompare( "integer", ts1 ) ) { 841 for( int n=0; n < nUnits ; n++ ) 842 fscanf(fileObject, "%d\n",(int*)((int*)valueArray+n) ); 843 } else if ( cscompare( "double", ts1 ) ) { 844 for( int n=0; n < nUnits ; n++ ) 845 fscanf(fileObject, "%lf\n",(double*)((double*)valueArray+n) ); 846 } 847 delete [] ts1; 848 } 849 return; 850 } 851 852 // End of copy from phastaIO 853 854 855 vtkPhastaReader::vtkPhastaReader() 856 { 857 //this->DebugOn(); // TODO: comment out this line to turn off debug 858 this->GeometryFileName = NULL; 859 this->FieldFileName = NULL; 860 this->SetNumberOfInputPorts(0); 861 this->Internal = new vtkPhastaReaderInternal; 862 863 ////////// 864 this->Parser = 0; 865 this->FileName = 0; 866 ////////// 867 868 } 869 870 vtkPhastaReader::~vtkPhastaReader() 871 { 872 if (this->GeometryFileName) 873 { 874 delete [] this->GeometryFileName; 875 } 876 if (this->FieldFileName) 877 { 878 delete [] this->FieldFileName; 879 } 880 delete this->Internal; 881 882 //////////////////// 883 if (this->Parser) 884 this->Parser->Delete(); 885 //////////////////// 886 887 } 888 889 void vtkPhastaReader::ClearFieldInfo() 890 { 891 this->Internal->FieldInfoMap.clear(); 892 } 893 894 void vtkPhastaReader::SetFieldInfo(const char* paraviewFieldTag, 895 const char* phastaFieldTag, 896 int index, 897 int numOfComps, 898 int dataDependency, 899 const char* dataType) 900 { 901 //printf("In P setfino\n"); 902 903 //CHANGE///////////////////////// 904 partID_counter=0; 905 906 //CHANGE END///////////////////// 907 908 vtkPhastaReaderInternal::FieldInfo &info = 909 this->Internal->FieldInfoMap[paraviewFieldTag]; 910 911 info.PhastaFieldTag = phastaFieldTag; 912 info.StartIndexInPhastaArray = index; 913 info.NumberOfComponents = numOfComps; 914 info.DataDependency = dataDependency; 915 info.DataType = dataType; 916 } 917 918 int vtkPhastaReader::RequestData(vtkInformation*, 919 vtkInformationVector**, 920 vtkInformationVector* outputVector) 921 { 922 vtkDebugMacro("In P RequestData"); 923 924 int firstVertexNo = 0; 925 int fvn = 0; 926 int noOfNodes, noOfCells, noOfDatas; 927 928 929 // get the data object 930 //TODO Just Testing 931 vtkSmartPointer<vtkInformation> outInfo = 932 outputVector->GetInformationObject(0); 933 934 //change/////// This part not working 935 // get the current piece being requested 936 int piece = 937 outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER()); 938 939 //printf("piece is %d\n",piece); 940 941 //partID_counter++; 942 partID_counter=PART_ID; 943 //change end//////////////// 944 945 vtkUnstructuredGrid *output = vtkUnstructuredGrid::SafeDownCast( 946 outInfo->Get(vtkDataObject::DATA_OBJECT())); 947 948 //printf("This geom file is %s\n",this->GeometryFileName); 949 950 951 int numPieces=NUM_PIECES, numFiles=NUM_FILES, timeStep=TIME_STEP; 952 //int numPieces=8, numFiles=4, timeStep=50400; 953 954 int numPiecesPerFile = numPieces/numFiles; 955 int fileID; 956 fileID = int((partID_counter-1)/numPiecesPerFile)+1; 957 958 vtkDebugMacro(<< "FILE_PATH: " << FILE_PATH); 959 // FILE_PATH is set to be the path of .pht file ? 960 //sprintf(this->FieldFileName,"%s%s.%d.%d",FILE_PATH,"/restart-dat",timeStep,fileID); // this is Ning's version 961 962 // the file id of fieldfilename need to be changed to file id now 963 char* str = this->FieldFileName; 964 for (int i = strlen(this->FieldFileName); i >= 0; i--) { 965 if(str[i] != '.') 966 str[i] = 0; 967 else 968 break; 969 } 970 sprintf(str, "%s%d", str, FILE_ID); 971 vtkDebugMacro(<<"tweaked FieldFileName="<<this->FieldFileName); 972 /////////////////////////////////////////////////////////// 973 974 vtkPoints *points; 975 976 output->Allocate(10000, 2100); 977 978 points = vtkPoints::New(); 979 980 vtkDebugMacro(<<"Reading Phasta file..."); 981 982 if(!this->GeometryFileName || !this->FieldFileName ) 983 { 984 vtkErrorMacro(<<"All input parameters not set."); 985 return 0; 986 } 987 vtkDebugMacro(<< "Updating ensa with ...."); 988 vtkDebugMacro(<< "Geom File : " << this->GeometryFileName); 989 vtkDebugMacro(<< "Field File : " << this->FieldFileName); 990 991 fvn = firstVertexNo; 992 this->ReadGeomFile(this->GeometryFileName, firstVertexNo, points, noOfNodes, noOfCells); 993 /* set the points over here, this is because vtkUnStructuredGrid 994 only insert points once, next insertion overwrites the previous one */ 995 // acbauer is not sure why the above comment is about... 996 output->SetPoints(points); 997 points->Delete(); 998 999 if (!this->Internal->FieldInfoMap.size()) 1000 { 1001 vtkDataSetAttributes* field = output->GetPointData(); 1002 this->ReadFieldFile(this->FieldFileName, fvn, field, noOfNodes); 1003 } 1004 else 1005 { 1006 this->ReadFieldFile(this->FieldFileName, fvn, output, noOfDatas); 1007 } 1008 1009 // if there exists point arrays called coordsX, coordsY and coordsZ, 1010 // create another array of point data and set the output to use this 1011 vtkPointData* pointData = output->GetPointData(); 1012 vtkDoubleArray* coordsX = vtkDoubleArray::SafeDownCast( 1013 pointData->GetArray("coordsX")); 1014 vtkDoubleArray* coordsY = vtkDoubleArray::SafeDownCast( 1015 pointData->GetArray("coordsY")); 1016 vtkDoubleArray* coordsZ = vtkDoubleArray::SafeDownCast( 1017 pointData->GetArray("coordsZ")); 1018 if(coordsX && coordsY && coordsZ) 1019 { 1020 vtkIdType numPoints = output->GetPoints()->GetNumberOfPoints(); 1021 if(numPoints != coordsX->GetNumberOfTuples() || 1022 numPoints != coordsY->GetNumberOfTuples() || 1023 numPoints != coordsZ->GetNumberOfTuples() ) 1024 { 1025 vtkWarningMacro("Wrong number of points for moving mesh. Using original points."); 1026 return 0; 1027 } 1028 vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New(); 1029 points->DeepCopy(output->GetPoints()); 1030 for(vtkIdType i=0;i<numPoints;i++) 1031 { 1032 points->SetPoint(i, coordsX->GetValue(i), coordsY->GetValue(i), 1033 coordsZ->GetValue(i)); 1034 } 1035 output->SetPoints(points); 1036 } 1037 1038 //printf("end of RequestData():"); 1039 //printf("counter = %ld\n", counter++); 1040 1041 vtkDebugMacro("end of P RequestData() # of pieces is "<<numPieces << ", partID_counter = "<< partID_counter); 1042 vtkDebugMacro("-------> total open time is "<< opentime_total); 1043 1044 return 1; 1045 } 1046 1047 /* firstVertexNo is useful when reading multiple geom files and coalescing 1048 them into one, ReadGeomfile can then be called repeatedly from Execute with 1049 firstVertexNo forming consecutive series of vertex numbers */ 1050 1051 void vtkPhastaReader::ReadGeomFile(char* geomFileName, 1052 int &firstVertexNo, 1053 vtkPoints *points, 1054 int &num_nodes, 1055 int &num_cells) 1056 { 1057 vtkDebugMacro("in P ReadGeomFile(): partID="<<partID_counter); 1058 1059 /* variables for vtk */ 1060 vtkUnstructuredGrid *output = this->GetOutput(); 1061 double *coordinates; 1062 vtkIdType *nodes; 1063 int cell_type; 1064 1065 // int num_tpblocks; 1066 /* variables for the geom data file */ 1067 /* nodal information */ 1068 // int byte_order; 1069 //int data[11], data1[7]; 1070 int dim; 1071 int num_int_blocks; 1072 double *pos; 1073 //int *nlworkdata; 1074 /* element information */ 1075 int num_elems,num_vertices,num_per_line; 1076 int *connectivity = NULL; 1077 1078 1079 /* misc variables*/ 1080 int i, j,k,item; 1081 int geomfile; 1082 1083 //CHANGE/////////////////////////////////////// 1084 int getNfiles, getNPPF; 1085 queryphmpiio_(geomFileName,&getNfiles,&getNPPF); 1086 char fieldName[255]; 1087 1088 1089 //CHANGE END/////////////////////////////////// 1090 1091 openfile(geomFileName,"read",&geomfile); 1092 //geomfile = fopen(GeometryFileName,"rb"); 1093 1094 if(!geomfile) 1095 { 1096 vtkErrorMacro(<<"Cannot open file " << geomFileName); 1097 //return; 1098 } 1099 1100 int expect; 1101 int array[10]; 1102 expect = 1; 1103 1104 /* read number of nodes */ 1105 1106 ///CHANGE///////////////////////////////////////////////////// 1107 bzero((void*)fieldName,255); 1108 sprintf(fieldName,"%s@%d","number of nodes",partID_counter); 1109 ///CHANGE END////////////////////////////////////////////////// 1110 readheader(&geomfile,fieldName,array,&expect,"integer","binary"); 1111 //readheader(&geomfile,"number of nodes",array,&expect,"integer","binary"); 1112 vtkDebugMacro("after readheader(), fieldName=" << fieldName << ", geomfile (file desc) = " << geomfile); 1113 num_nodes = array[0]; 1114 1115 /* read number of elements */ 1116 1117 ///CHANGE///////////////////////////////////////////////////// 1118 bzero((void*)fieldName,255); 1119 sprintf(fieldName,"%s@%d","number of interior elements",partID_counter); 1120 ///CHANGE END////////////////////////////////////////////////// 1121 readheader(&geomfile,fieldName,array,&expect,"integer","binary"); 1122 1123 /*readheader(&geomfile, 1124 "number of interior elements", 1125 array, 1126 &expect, 1127 "integer", 1128 "binary"); 1129 */ 1130 num_elems = array[0]; 1131 num_cells = array[0]; 1132 1133 /* read number of interior */ 1134 1135 ///CHANGE///////////////////////////////////////////////////// 1136 bzero((void*)fieldName,255); 1137 sprintf(fieldName,"%s@%d","number of interior tpblocks",partID_counter); 1138 ///CHANGE END////////////////////////////////////////////////// 1139 readheader(&geomfile,fieldName,array,&expect,"integer","binary"); 1140 1141 1142 /*readheader(&geomfile, 1143 "number of interior tpblocks", 1144 array, 1145 &expect, 1146 "integer", 1147 "binary"); 1148 */ 1149 num_int_blocks = array[0]; 1150 1151 vtkDebugMacro ( << "Nodes: " << num_nodes 1152 << "Elements: " << num_elems 1153 << "tpblocks: " << num_int_blocks ); 1154 1155 /* read coordinates */ 1156 expect = 2; 1157 1158 ///CHANGE///////////////////////////////////////////////////// 1159 bzero((void*)fieldName,255); 1160 sprintf(fieldName,"%s@%d","co-ordinates",partID_counter); 1161 ///CHANGE END////////////////////////////////////////////////// 1162 readheader(&geomfile,fieldName,array,&expect,"double","binary"); 1163 1164 1165 //readheader(&geomfile,"co-ordinates",array,&expect,"double","binary"); 1166 // TEST ******************* 1167 num_nodes=array[0]; 1168 // TEST ******************* 1169 if(num_nodes !=array[0]) 1170 { 1171 vtkErrorMacro(<<"Ambigous information in geom.data file, number of nodes does not match the co-ordinates size. Nodes: " << num_nodes << " Coordinates: " << array[0]); 1172 return; 1173 } 1174 dim = array[1]; 1175 1176 1177 /* read the coordinates */ 1178 1179 coordinates = new double [dim]; 1180 if(coordinates == NULL) 1181 { 1182 vtkErrorMacro(<<"Unable to allocate memory for nodal info"); 1183 return; 1184 } 1185 1186 pos = new double [num_nodes*dim]; 1187 if(pos == NULL) 1188 { 1189 vtkErrorMacro(<<"Unable to allocate memory for nodal info"); 1190 return; 1191 } 1192 1193 item = num_nodes*dim; 1194 1195 //CHANGE 1196 //readdatablock(&geomfile,"co-ordinates",pos,&item,"double","binary"); 1197 //CHANGE END 1198 readdatablock(&geomfile,fieldName,pos,&item,"double","binary"); 1199 1200 for(i=0;i<num_nodes;i++) 1201 { 1202 for(j=0;j<dim;j++) 1203 { 1204 coordinates[j] = pos[j*num_nodes + i]; 1205 } 1206 switch(dim) 1207 { 1208 case 1: 1209 points->InsertPoint(i+firstVertexNo,coordinates[0],0,0); 1210 break; 1211 case 2: 1212 points->InsertPoint(i+firstVertexNo,coordinates[0],coordinates[1],0); 1213 break; 1214 case 3: 1215 points->InsertNextPoint(coordinates); 1216 break; 1217 default: 1218 vtkErrorMacro(<<"Unrecognized dimension in "<< geomFileName) 1219 return; 1220 } 1221 } 1222 1223 /* read the connectivity information */ 1224 expect = 7; 1225 1226 for(k=0;k<num_int_blocks;k++) 1227 { 1228 1229 ///CHANGE///////////////////////////////////////////////////// 1230 bzero((void*)fieldName,255); 1231 sprintf(fieldName,"%s%d@%d","connectivity interior",k+1,partID_counter); 1232 ///CHANGE END////////////////////////////////////////////////// 1233 1234 readheader(&geomfile, 1235 fieldName, 1236 array, 1237 &expect, 1238 "integer", 1239 "binary"); 1240 1241 /*readheader(&geomfile, 1242 "connectivity interior", 1243 array, 1244 &expect, 1245 "integer", 1246 "binary"); 1247 */ 1248 1249 /* read information about the block*/ 1250 num_elems = array[0]; 1251 num_vertices = array[1]; 1252 num_per_line = array[3]; 1253 connectivity = new int [num_elems*num_per_line]; 1254 1255 if(connectivity == NULL) 1256 { 1257 vtkErrorMacro(<<"Unable to allocate memory for connectivity info"); 1258 return; 1259 } 1260 1261 item = num_elems*num_per_line; 1262 /*readdatablock(&geomfile, 1263 "connectivity interior", 1264 connectivity, 1265 &item, 1266 "integer", 1267 "binary"); 1268 */ 1269 1270 readdatablock(&geomfile, 1271 fieldName, 1272 connectivity, 1273 &item, 1274 "integer", 1275 "binary"); 1276 1277 /* insert cells */ 1278 for(i=0;i<num_elems;i++) 1279 { 1280 nodes = new vtkIdType[num_vertices]; 1281 1282 //connectivity starts from 1 so node[j] will never be -ve 1283 for(j=0;j<num_vertices;j++) 1284 { 1285 nodes[j] = connectivity[i+num_elems*j] + firstVertexNo - 1; 1286 } 1287 1288 /* 1 is subtracted from the connectivity info to reflect that in vtk 1289 vertex numbering start from 0 as opposed to 1 in geomfile */ 1290 1291 // find out element type 1292 switch(num_vertices) 1293 { 1294 case 4: 1295 cell_type = VTK_TETRA; 1296 break; 1297 case 5: 1298 cell_type = VTK_PYRAMID; 1299 break; 1300 case 6: 1301 cell_type = VTK_WEDGE; 1302 break; 1303 case 8: 1304 cell_type = VTK_HEXAHEDRON; 1305 1306 break; 1307 default: 1308 vtkErrorMacro(<<"Unrecognized CELL_TYPE in "<< geomFileName) 1309 return; 1310 } 1311 1312 /* insert the element */ 1313 output->InsertNextCell(cell_type,num_vertices,nodes); 1314 delete [] nodes; 1315 } 1316 } 1317 // update the firstVertexNo so that next slice/partition can be read 1318 firstVertexNo = firstVertexNo + num_nodes; 1319 1320 // clean up 1321 closefile(&geomfile,"read"); 1322 finalizephmpiio_(&geomfile); 1323 delete [] coordinates; 1324 delete [] pos; 1325 delete [] connectivity; 1326 1327 } // end of ReadGeomFile 1328 1329 void vtkPhastaReader::ReadFieldFile(char* fieldFileName, 1330 int, 1331 vtkDataSetAttributes *field, 1332 int &noOfNodes) 1333 { 1334 1335 vtkDebugMacro("In P ReadFieldFile (vtkDataSetAttr), readheader etc calls need to be updated.."); 1336 1337 int i, j; 1338 int item; 1339 double *data; 1340 int fieldfile; 1341 1342 int getNfiles, getNPPF; 1343 queryphmpiio_(fieldFileName,&getNfiles,&getNPPF); 1344 1345 openfile(fieldFileName,"read",&fieldfile); 1346 //fieldfile = fopen(FieldFileName,"rb"); 1347 1348 if(!fieldfile) 1349 { 1350 vtkErrorMacro(<<"Cannot open file " << FieldFileName) 1351 return; 1352 } 1353 int array[10], expect; 1354 1355 /* read the solution */ 1356 vtkDoubleArray* pressure = vtkDoubleArray::New(); 1357 pressure->SetName("pressure"); 1358 vtkDoubleArray* velocity = vtkDoubleArray::New(); 1359 velocity->SetName("velocity"); 1360 velocity->SetNumberOfComponents(3); 1361 vtkDoubleArray* temperature = vtkDoubleArray::New(); 1362 temperature->SetName("temperature"); 1363 1364 expect = 3; 1365 readheader(&fieldfile,"solution",array,&expect,"double","binary"); 1366 noOfNodes = array[0]; 1367 this->NumberOfVariables = array[1]; 1368 1369 vtkDoubleArray* sArrays[4]; 1370 for (i=0; i<4; i++) 1371 { 1372 sArrays[i] = 0; 1373 } 1374 item = noOfNodes*this->NumberOfVariables; 1375 data = new double[item]; 1376 if(data == NULL) 1377 { 1378 vtkErrorMacro(<<"Unable to allocate memory for field info"); 1379 return; 1380 } 1381 1382 readdatablock(&fieldfile,"solution",data,&item,"double","binary"); 1383 1384 for (i=5; i<this->NumberOfVariables; i++) 1385 { 1386 int idx=i-5; 1387 sArrays[idx] = vtkDoubleArray::New(); 1388 vtksys_ios::ostringstream aName; 1389 aName << "s" << idx+1 << ends; 1390 sArrays[idx]->SetName(aName.str().c_str()); 1391 sArrays[idx]->SetNumberOfTuples(noOfNodes); 1392 } 1393 1394 pressure->SetNumberOfTuples(noOfNodes); 1395 velocity->SetNumberOfTuples(noOfNodes); 1396 temperature->SetNumberOfTuples(noOfNodes); 1397 for(i=0;i<noOfNodes;i++) 1398 { 1399 pressure->SetTuple1(i, data[i]); 1400 velocity->SetTuple3(i, 1401 data[noOfNodes + i], 1402 data[2*noOfNodes + i], 1403 data[3*noOfNodes + i]); 1404 temperature->SetTuple1(i, data[4*noOfNodes + i]); 1405 for(j=5; j<this->NumberOfVariables; j++) 1406 { 1407 sArrays[j-5]->SetTuple1(i, data[j*noOfNodes + i]); 1408 } 1409 } 1410 1411 field->AddArray(pressure); 1412 field->SetActiveScalars("pressure"); 1413 pressure->Delete(); 1414 field->AddArray(velocity); 1415 field->SetActiveVectors("velocity"); 1416 velocity->Delete(); 1417 field->AddArray(temperature); 1418 temperature->Delete(); 1419 1420 for (i=5; i<this->NumberOfVariables; i++) 1421 { 1422 int idx=i-5; 1423 field->AddArray(sArrays[idx]); 1424 sArrays[idx]->Delete(); 1425 } 1426 1427 // clean up 1428 closefile(&fieldfile,"read"); 1429 finalizephmpiio_(&fieldfile); 1430 delete [] data; 1431 1432 } //closes ReadFieldFile (vtkDataSetAttr) 1433 1434 1435 void vtkPhastaReader::ReadFieldFile(char* fieldFileName, 1436 int, 1437 vtkUnstructuredGrid *output, 1438 int &noOfDatas) 1439 { 1440 1441 vtkDebugMacro("In P ReadFieldFile (vtkUnstructuredGrid): fieldFileName="<<fieldFileName << ", partID_counter=" << partID_counter); 1442 1443 int i, j, numOfVars; 1444 int item; 1445 int fieldfile; 1446 1447 //printf("In P readfieldfile 2 and %d\n",partID_counter); 1448 //printf("read filed file name is %s\n",fieldFileName); 1449 1450 //CHANGE/////////////////////////////////////// 1451 int getNfiles, getNPPF; 1452 queryphmpiio_(fieldFileName,&getNfiles,&getNPPF); 1453 char fieldName[255]; 1454 1455 //CHANGE END/////////////////////////////////// 1456 1457 openfile(fieldFileName,"read",&fieldfile); 1458 //fieldfile = fopen(FieldFileName,"rb"); 1459 1460 //printf("open handle is %d\n",fieldfile); 1461 1462 if(!fieldfile) 1463 { 1464 vtkErrorMacro(<<"Cannot open file " << FieldFileName) 1465 //return; 1466 } 1467 int array[10], expect; 1468 1469 int activeScalars = 0, activeTensors = 0; 1470 1471 vtkPhastaReaderInternal::FieldInfoMapType::iterator it = this->Internal->FieldInfoMap.begin(); 1472 vtkPhastaReaderInternal::FieldInfoMapType::iterator itend = this->Internal->FieldInfoMap.end(); 1473 for(; it!=itend; it++) 1474 { 1475 const char* paraviewFieldTag = it->first.c_str(); 1476 const char* phastaFieldTag = it->second.PhastaFieldTag.c_str(); 1477 int index = it->second.StartIndexInPhastaArray; 1478 int numOfComps = it->second.NumberOfComponents; 1479 int dataDependency = it->second.DataDependency; 1480 const char* dataType = it->second.DataType.c_str(); 1481 1482 1483 vtkDataSetAttributes* field; 1484 if(dataDependency) 1485 field = output->GetCellData(); 1486 else 1487 field = output->GetPointData(); 1488 1489 // void *data; 1490 int dtype; // (0=double, 1=float) 1491 vtkDataArray *dataArray; 1492 /* read the field data */ 1493 if(strcmp(dataType,"double")==0) 1494 { 1495 dataArray = vtkDoubleArray::New(); 1496 dtype = 0; 1497 } 1498 else if(strcmp(dataType,"float")==0) 1499 { 1500 dataArray = vtkFloatArray::New(); 1501 dtype=1; 1502 } 1503 else 1504 { 1505 vtkErrorMacro("Data type [" << dataType <<"] NOT supported"); 1506 continue; 1507 } 1508 1509 dataArray->SetName(paraviewFieldTag); 1510 dataArray->SetNumberOfComponents(numOfComps); 1511 1512 expect = 3; 1513 1514 //printf("version 1: hello 3 in P\n"); 1515 1516 ///CHANGE///////////////////////////////////////////////////// 1517 bzero((void*)fieldName,255); 1518 sprintf(fieldName,"%s@%d",phastaFieldTag,partID_counter); 1519 ///CHANGE END////////////////////////////////////////////////// 1520 readheader(&fieldfile,fieldName,array,&expect,dataType,"binary"); 1521 1522 //printf("fieldfile is %d expect is %d and datatype is %s\n",fieldfile,expect,dataType); 1523 1524 //printf("original fieldname is %s and fieldname is %s\n",phastaFieldTag,fieldName); 1525 //printf("array 0 is %d and array 1 is %d\n",array[0],array[1]); 1526 1527 //readheader(&fieldfile,phastaFieldTag,array,&expect,dataType,"binary"); 1528 noOfDatas = array[0]; 1529 this->NumberOfVariables = array[1]; 1530 numOfVars = array[1]; 1531 dataArray->SetNumberOfTuples(noOfDatas); 1532 1533 //printf("hello 4 in P\n"); 1534 1535 if(index<0 || index>numOfVars-1) 1536 { 1537 vtkErrorMacro("index ["<<index<<"] is out of range [num. of vars.:"<<numOfVars<<"] for field [paraview field tag:"<<paraviewFieldTag<<", phasta field tag:"<<phastaFieldTag<<"]"); 1538 1539 dataArray->Delete(); 1540 continue; 1541 } 1542 1543 //printf("hello 5 in P\n"); 1544 1545 if(numOfComps<0 || index+numOfComps>numOfVars) 1546 { 1547 vtkErrorMacro("index ["<<index<<"] with num. of comps. ["<<numOfComps<<"] is out of range [num. of vars.:"<<numOfVars<<"] for field [paraview field tag:"<<paraviewFieldTag<<", phasta field tag:"<<phastaFieldTag<<"]"); 1548 1549 dataArray->Delete(); 1550 continue; 1551 } 1552 1553 //printf("hello 6 in P\n"); 1554 1555 item = numOfVars*noOfDatas; 1556 1557 //printf("item is %d and dtype is %d and numofvars is %d num of data is %d\n",item,dtype,numOfVars,noOfDatas); 1558 1559 if (dtype==0) 1560 { //data is type double 1561 1562 double *data; 1563 data = new double[item]; 1564 1565 if(data == NULL) 1566 { 1567 vtkErrorMacro(<<"Unable to allocate memory for field info"); 1568 1569 dataArray->Delete(); 1570 continue; 1571 } 1572 1573 //printf("hello 8 in P\n"); 1574 1575 //////////CHANGE///////////////////// 1576 readdatablock(&fieldfile,fieldName,data,&item,dataType,"binary"); 1577 ///////////CHANGE END/////////////// 1578 1579 //readdatablock(&fieldfile,phastaFieldTag,data,&item,dataType,"binary"); 1580 1581 //printf("hello 9 in P\n"); 1582 1583 switch(numOfComps) 1584 { 1585 case 1 : 1586 { 1587 int offset = index*noOfDatas; 1588 1589 if(!activeScalars) 1590 field->SetActiveScalars(paraviewFieldTag); 1591 else 1592 activeScalars = 1; 1593 for(i=0;i<noOfDatas;i++) 1594 { 1595 dataArray->SetTuple1(i, data[offset+i]); 1596 } 1597 } 1598 break; 1599 case 3 : 1600 { 1601 int offset[3]; 1602 for(j=0;j<3;j++) 1603 offset[j] = (index+j)*noOfDatas; 1604 1605 if(!activeScalars) 1606 field->SetActiveVectors(paraviewFieldTag); 1607 else 1608 activeScalars = 1; 1609 for(i=0;i<noOfDatas;i++) 1610 { 1611 dataArray->SetTuple3(i, 1612 data[offset[0]+i], 1613 data[offset[1]+i], 1614 data[offset[2]+i]); 1615 } 1616 } 1617 break; 1618 case 9 : 1619 { 1620 int offset[9]; 1621 for(j=0;j<9;j++) 1622 offset[j] = (index+j)*noOfDatas; 1623 1624 if(!activeTensors) 1625 field->SetActiveTensors(paraviewFieldTag); 1626 else 1627 activeTensors = 1; 1628 for(i=0;i<noOfDatas;i++) 1629 { 1630 dataArray->SetTuple9(i, 1631 data[offset[0]+i], 1632 data[offset[1]+i], 1633 data[offset[2]+i], 1634 data[offset[3]+i], 1635 data[offset[4]+i], 1636 data[offset[5]+i], 1637 data[offset[6]+i], 1638 data[offset[7]+i], 1639 data[offset[8]+i]); 1640 } 1641 } 1642 break; 1643 default: 1644 vtkErrorMacro("number of components [" << numOfComps <<"] NOT supported"); 1645 1646 dataArray->Delete(); 1647 delete [] data; 1648 continue; 1649 } 1650 1651 // clean up 1652 delete [] data; 1653 1654 } 1655 else if(dtype==1) 1656 { // data is type float 1657 1658 float *data; 1659 data = new float[item]; 1660 1661 if(data == NULL) 1662 { 1663 vtkErrorMacro(<<"Unable to allocate memory for field info"); 1664 1665 dataArray->Delete(); 1666 continue; 1667 } 1668 1669 //////////CHANGE///////////////////// 1670 readdatablock(&fieldfile,fieldName,data,&item,dataType,"binary"); 1671 //////////CHANGE END///////////////// 1672 1673 //readdatablock(&fieldfile,phastaFieldTag,data,&item,dataType,"binary"); 1674 1675 switch(numOfComps) 1676 { 1677 case 1 : 1678 { 1679 int offset = index*noOfDatas; 1680 1681 if(!activeScalars) 1682 field->SetActiveScalars(paraviewFieldTag); 1683 else 1684 activeScalars = 1; 1685 for(i=0;i<noOfDatas;i++) 1686 { 1687 //double tmpval = (double) data[offset+i]; 1688 //dataArray->SetTuple1(i, tmpval); 1689 dataArray->SetTuple1(i, data[offset+i]); 1690 } 1691 } 1692 break; 1693 case 3 : 1694 { 1695 int offset[3]; 1696 for(j=0;j<3;j++) 1697 offset[j] = (index+j)*noOfDatas; 1698 1699 if(!activeScalars) 1700 field->SetActiveVectors(paraviewFieldTag); 1701 else 1702 activeScalars = 1; 1703 for(i=0;i<noOfDatas;i++) 1704 { 1705 //double tmpval[3]; 1706 //for(j=0;j<3;j++) 1707 // tmpval[j] = (double) data[offset[j]+i]; 1708 //dataArray->SetTuple3(i, 1709 // tmpval[0], tmpval[1], tmpval[3]); 1710 dataArray->SetTuple3(i, 1711 data[offset[0]+i], 1712 data[offset[1]+i], 1713 data[offset[2]+i]); 1714 } 1715 } 1716 break; 1717 case 9 : 1718 { 1719 int offset[9]; 1720 for(j=0;j<9;j++) 1721 offset[j] = (index+j)*noOfDatas; 1722 1723 if(!activeTensors) 1724 field->SetActiveTensors(paraviewFieldTag); 1725 else 1726 activeTensors = 1; 1727 for(i=0;i<noOfDatas;i++) 1728 { 1729 //double tmpval[9]; 1730 //for(j=0;j<9;j++) 1731 // tmpval[j] = (double) data[offset[j]+i]; 1732 //dataArray->SetTuple9(i, 1733 // tmpval[0], 1734 // tmpval[1], 1735 // tmpval[2], 1736 // tmpval[3], 1737 // tmpval[4], 1738 // tmpval[5], 1739 // tmpval[6], 1740 // tmpval[7], 1741 // tmpval[8]); 1742 dataArray->SetTuple9(i, 1743 data[offset[0]+i], 1744 data[offset[1]+i], 1745 data[offset[2]+i], 1746 data[offset[3]+i], 1747 data[offset[4]+i], 1748 data[offset[5]+i], 1749 data[offset[6]+i], 1750 data[offset[7]+i], 1751 data[offset[8]+i]); 1752 } 1753 } 1754 break; 1755 default: 1756 vtkErrorMacro("number of components [" << numOfComps <<"] NOT supported"); 1757 1758 dataArray->Delete(); 1759 delete [] data; 1760 continue; 1761 } 1762 1763 // clean up 1764 delete [] data; 1765 1766 } 1767 else 1768 { 1769 vtkErrorMacro("Data type [" << dataType <<"] NOT supported"); 1770 continue; 1771 } 1772 1773 1774 field->AddArray(dataArray); 1775 1776 // clean up 1777 dataArray->Delete(); 1778 //delete [] data; 1779 } 1780 1781 // close up 1782 closefile(&fieldfile,"read"); 1783 finalizephmpiio_(&fieldfile); 1784 1785 }//closes ReadFieldFile (vtkUnstructuredGrid) 1786 1787 void vtkPhastaReader::PrintSelf(ostream& os, vtkIndent indent) 1788 { 1789 this->Superclass::PrintSelf(os,indent); 1790 1791 os << indent << "GeometryFileName: " 1792 << (this->GeometryFileName?this->GeometryFileName:"(none)") 1793 << endl; 1794 os << indent << "FieldFileName: " 1795 << (this->FieldFileName?this->FieldFileName:"(none)") 1796 << endl; 1797 } 1798