vtkPLOT3DReader.cxx 53.3 KB
Newer Older
Will Schroeder's avatar
Will Schroeder committed
1 2
/*=========================================================================

Ken Martin's avatar
Ken Martin committed
3
  Program:   Visualization Toolkit
Ken Martin's avatar
Ken Martin committed
4
  Module:    vtkPLOT3DReader.cxx
Will Schroeder's avatar
Will Schroeder committed
5

6
  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
7 8
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
Ken Martin's avatar
Ken Martin committed
9

10 11
     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12
     PURPOSE.  See the above copyright notice for more information.
Will Schroeder's avatar
Will Schroeder committed
13 14

=========================================================================*/
Ken Martin's avatar
Ken Martin committed
15
#include "vtkPLOT3DReader.h"
Berk Geveci's avatar
Berk Geveci committed
16

Ken Martin's avatar
Ken Martin committed
17
#include "vtkByteSwap.h"
Berk Geveci's avatar
Berk Geveci committed
18
#include "vtkErrorCode.h"
19
#include "vtkFieldData.h"
20
#include "vtkFloatArray.h"
Berk Geveci's avatar
Berk Geveci committed
21 22
#include "vtkIntArray.h"
#include "vtkObjectFactory.h"
23
#include "vtkPointData.h"
24
#include "vtkStructuredGrid.h"
Berk Geveci's avatar
Berk Geveci committed
25
#include "vtkUnsignedCharArray.h"
26

Berk Geveci's avatar
Berk Geveci committed
27
vtkCxxRevisionMacro(vtkPLOT3DReader, "1.83");
Brad King's avatar
Brad King committed
28
vtkStandardNewMacro(vtkPLOT3DReader);
29

Bill Lorensen's avatar
Bill Lorensen committed
30 31 32 33
#define VTK_RHOINF 1.0
#define VTK_CINF 1.0
#define VTK_PINF ((VTK_RHOINF*VTK_CINF) * (VTK_RHOINF*VTK_CINF) / this->Gamma)
#define VTK_CV (this->R / (this->Gamma-1.0))
Will Schroeder's avatar
Will Schroeder committed
34

Ken Martin's avatar
Ken Martin committed
35
vtkPLOT3DReader::vtkPLOT3DReader()
Will Schroeder's avatar
Will Schroeder committed
36
{
Will Schroeder's avatar
Will Schroeder committed
37 38
  this->XYZFileName = NULL;
  this->QFileName = NULL;
Berk Geveci's avatar
Berk Geveci committed
39 40 41 42 43 44 45
  this->BinaryFile = 1;
  this->HasByteCount = 0;
  this->FileSize = 0;
  this->MultiGrid = 0;
  this->ForceRead = 0;
  this->ByteOrder = FILE_BIG_ENDIAN;
  this->IBlanking = 0;
46
  this->TwoDimensionalGeometry = 0;
Berk Geveci's avatar
Berk Geveci committed
47
  this->DoNotReduceNumberOfOutputs = 1;
Will Schroeder's avatar
Will Schroeder committed
48 49 50 51 52 53 54

  this->R = 1.0;
  this->Gamma = 1.4;
  this->Uvinf = 0.0;
  this->Vvinf = 0.0;
  this->Wvinf = 0.0;

Berk Geveci's avatar
Berk Geveci committed
55
  this->FunctionList = vtkIntArray::New();
56

Berk Geveci's avatar
Berk Geveci committed
57 58 59 60
  this->ScalarFunctionNumber = -1;
  this->SetScalarFunctionNumber(100);
  this->VectorFunctionNumber = -1;
  this->SetVectorFunctionNumber(202);
61 62 63

  this->PointCache = 0;
  this->IBlankCache = 0;
64
} 
Will Schroeder's avatar
Will Schroeder committed
65

Ken Martin's avatar
Ken Martin committed
66
vtkPLOT3DReader::~vtkPLOT3DReader()
Will Schroeder's avatar
Will Schroeder committed
67
{
Berk Geveci's avatar
Berk Geveci committed
68 69
  delete [] this->XYZFileName;
  delete [] this->QFileName;
70
  this->FunctionList->Delete();
71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92
  this->ClearGeometryCache();
}

void vtkPLOT3DReader::ClearGeometryCache()
{
  if ( this->PointCache )
    {
    for ( int g=0; this->PointCache[g]; ++g )
      this->PointCache[g]->UnRegister( this );

    delete [] this->PointCache;
    this->PointCache = 0;
    }

  if ( this->IBlankCache )
    {
    for ( int i=0; this->IBlankCache[i]; ++i )
      this->IBlankCache[i]->UnRegister( this );

    delete [] this->IBlankCache;
    this->IBlankCache = 0;
    }
Will Schroeder's avatar
Will Schroeder committed
93 94
}

Berk Geveci's avatar
Berk Geveci committed
95
int vtkPLOT3DReader::CheckFile(FILE*& fp, const char* fname)
96
{
97 98 99 100 101 102 103 104 105
  if (this->BinaryFile)
    {
    fp = fopen(fname, "rb");
    }
  else
    {
    fp = fopen(fname, "r");
    }
  if ( fp == NULL)
106
    {
Berk Geveci's avatar
Berk Geveci committed
107 108 109
    this->SetErrorCode(vtkErrorCode::FileNotFoundError);
    vtkErrorMacro(<< "File: " << fname << " not found.");
    return VTK_ERROR;
110
    }
Berk Geveci's avatar
Berk Geveci committed
111
  return VTK_OK;
112
}
Charles Law's avatar
Charles Law committed
113

Berk Geveci's avatar
Berk Geveci committed
114
int vtkPLOT3DReader::CheckGeometryFile(FILE*& xyzFp)
115
{
Berk Geveci's avatar
Berk Geveci committed
116
  if ( this->XYZFileName == NULL || this->XYZFileName[0] == '\0'  )
117
    {
Berk Geveci's avatar
Berk Geveci committed
118
    this->SetErrorCode(vtkErrorCode::NoFileNameError);
119
    vtkErrorMacro(<< "Must specify geometry file");
Berk Geveci's avatar
Berk Geveci committed
120
    return VTK_ERROR;
121
    }
Berk Geveci's avatar
Berk Geveci committed
122
  return this->CheckFile(xyzFp, this->XYZFileName);
123 124
}

Berk Geveci's avatar
Berk Geveci committed
125
int vtkPLOT3DReader::CheckSolutionFile(FILE*& qFp)
Will Schroeder's avatar
Will Schroeder committed
126
{
Berk Geveci's avatar
Berk Geveci committed
127
  if ( this->QFileName == NULL || this->QFileName[0] == '\0' )
128
    {
Berk Geveci's avatar
Berk Geveci committed
129
    this->SetErrorCode(vtkErrorCode::NoFileNameError);
Ken Martin's avatar
Ken Martin committed
130
    vtkErrorMacro(<< "Must specify geometry file");
Berk Geveci's avatar
Berk Geveci committed
131
    return VTK_ERROR;
Ken Martin's avatar
Ken Martin committed
132
    }
Berk Geveci's avatar
Berk Geveci committed
133 134 135 136 137 138 139
  return this->CheckFile(qFp, this->QFileName);
}

// Skip Fortran style byte count.
void vtkPLOT3DReader::SkipByteCount(FILE* fp)
{
  if (this->BinaryFile && this->HasByteCount)
140
    {
Berk Geveci's avatar
Berk Geveci committed
141 142
    int tmp;
    fread(&tmp, sizeof(int), 1, fp);
143
    }
Berk Geveci's avatar
Berk Geveci committed
144
}
145

Berk Geveci's avatar
Berk Geveci committed
146 147 148 149
// Read a block of ints (ascii or binary) and return number read.
int vtkPLOT3DReader::ReadIntBlock(FILE* fp, int n, int* block)
{
  if (this->BinaryFile)
150
    {
Berk Geveci's avatar
Berk Geveci committed
151 152
    int retVal=static_cast<int>(fread(block, sizeof(int), n, fp));
    if (this->ByteOrder == FILE_LITTLE_ENDIAN)
LYMB Demo's avatar
LYMB Demo committed
153
      {
Berk Geveci's avatar
Berk Geveci committed
154
      vtkByteSwap::Swap4LERange(block, n);
LYMB Demo's avatar
LYMB Demo committed
155
      }
156
    else
Ken Martin's avatar
Ken Martin committed
157
      {
Berk Geveci's avatar
Berk Geveci committed
158
      vtkByteSwap::Swap4BERange(block, n);
159
      }
Berk Geveci's avatar
Berk Geveci committed
160 161 162 163 164 165
    return retVal;
    }
  else
    {
    int count = 0;
    for(int i=0; i<n; i++)
166
      {
Berk Geveci's avatar
Berk Geveci committed
167 168
      int num = fscanf(fp, "%d", &(block[i]));
      if ( num > 0 )
169
        {
Berk Geveci's avatar
Berk Geveci committed
170 171 172 173 174
        count++;
        }
      else
        {
        return 0;
175 176
        }
      }
Berk Geveci's avatar
Berk Geveci committed
177
    return count;
178
    }
Berk Geveci's avatar
Berk Geveci committed
179
}
180

Berk Geveci's avatar
Berk Geveci committed
181 182 183
int vtkPLOT3DReader::ReadFloatBlock(FILE* fp, int n, float* block)
{
  if (this->BinaryFile)
184
    {
Berk Geveci's avatar
Berk Geveci committed
185 186
    int retVal=static_cast<int>(fread(block, sizeof(float), n, fp));
    if (this->ByteOrder == FILE_LITTLE_ENDIAN)
187
      {
Berk Geveci's avatar
Berk Geveci committed
188
      vtkByteSwap::Swap4LERange(block, n);
LYMB Demo's avatar
LYMB Demo committed
189
      }
190
    else
Ken Martin's avatar
Ken Martin committed
191
      {
Berk Geveci's avatar
Berk Geveci committed
192
      vtkByteSwap::Swap4BERange(block, n);
193
      }
Berk Geveci's avatar
Berk Geveci committed
194
    return retVal;
195
    }
Berk Geveci's avatar
Berk Geveci committed
196
  else
197
    {
Berk Geveci's avatar
Berk Geveci committed
198 199
    int count = 0;
    for(int i=0; i<n; i++)
200
      {
Berk Geveci's avatar
Berk Geveci committed
201 202 203 204 205 206 207 208 209
      int num = fscanf(fp, "%f", &(block[i]));
      if ( num > 0 )
        {
        count++;
        }
      else
        {
        return 0;
        }
210
      }
Berk Geveci's avatar
Berk Geveci committed
211
    return count;
212
    }
Berk Geveci's avatar
Berk Geveci committed
213
}
214

Berk Geveci's avatar
Berk Geveci committed
215 216 217 218 219 220 221 222
// Read a block of floats (ascii or binary) and return number read.
void vtkPLOT3DReader::CalculateFileSize(FILE* fp)
{
  long curPos = ftell(fp);
  fseek(fp, 0, SEEK_END);
  this->FileSize = ftell(fp);
  fseek(fp, curPos, SEEK_SET);
}
Will Schroeder's avatar
Will Schroeder committed
223

224

Berk Geveci's avatar
Berk Geveci committed
225 226 227
// Estimate the size of a grid (binary file only)
long vtkPLOT3DReader::EstimateSize(int ni, int nj, int nk)
{
228 229 230 231 232 233 234 235 236 237 238
  long size; // the header portion, 3 ints
  if (!this->TwoDimensionalGeometry)
    {
    size = 3*4;
    size += ni*nj*nk*3*4; // x, y, z
    }
  else
    {
    size = 2*4;
    size += ni*nj*nk*2*4; // x, y, z
    }
Berk Geveci's avatar
Berk Geveci committed
239
  if (this->HasByteCount)
240
    {
Berk Geveci's avatar
Berk Geveci committed
241
    size += 2*4; // the byte counts
242
    }
Berk Geveci's avatar
Berk Geveci committed
243
  if (this->IBlanking)
244
    {
Berk Geveci's avatar
Berk Geveci committed
245
    size += ni*nj*nk*4;
246 247
    }

Berk Geveci's avatar
Berk Geveci committed
248
  return size;
Will Schroeder's avatar
Will Schroeder committed
249 250
}

251
int vtkPLOT3DReader::CanReadBinaryFile(const char* fname)
Will Schroeder's avatar
Will Schroeder committed
252
{
Berk Geveci's avatar
Berk Geveci committed
253 254 255
  FILE* xyzFp;

  if (!fname || fname[0] == '\0')
256
    {
Berk Geveci's avatar
Berk Geveci committed
257
    return 0;
258
    }
259

Berk Geveci's avatar
Berk Geveci committed
260
  if ( this->CheckFile(xyzFp, fname) != VTK_OK)
261
    {
Berk Geveci's avatar
Berk Geveci committed
262
    return 0;
263
    }
264 265 266

  this->CalculateFileSize(xyzFp);

Berk Geveci's avatar
Berk Geveci committed
267 268 269
  int numOutputs = this->GetNumberOfOutputsInternal(xyzFp, 1);
  fclose(xyzFp);
  if (numOutputs != 0)
270 271 272
    {
    return 1;
    }
Berk Geveci's avatar
Berk Geveci committed
273 274
  return 0;
}
275

Berk Geveci's avatar
Berk Geveci committed
276 277 278
int vtkPLOT3DReader::GetNumberOfOutputs()
{
  FILE* xyzFp;
279

Berk Geveci's avatar
Berk Geveci committed
280
  if ( this->CheckGeometryFile(xyzFp) != VTK_OK)
281
    {
Berk Geveci's avatar
Berk Geveci committed
282
    return 0;
283
    }
284
  this->CalculateFileSize(xyzFp);
Berk Geveci's avatar
Berk Geveci committed
285 286 287
  int numOutputs = this->GetNumberOfOutputsInternal(xyzFp, 1);
  fclose(xyzFp);
  if (numOutputs != 0)
288
    {
Berk Geveci's avatar
Berk Geveci committed
289
    return numOutputs;
290
    }
Berk Geveci's avatar
Berk Geveci committed
291
  return 1;
Will Schroeder's avatar
Will Schroeder committed
292 293
}

294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331
int vtkPLOT3DReader::GenerateDefaultConfiguration()
{
  FILE* xyzFp;
  
  if ( this->CheckGeometryFile(xyzFp) != VTK_OK)
    {
    return 0;
    }
  char buf[1024];
  fread(buf, 1, 1024, xyzFp);
  int retVal = this->VerifySettings(buf, 1024);
  fclose(xyzFp);
  return retVal;
}

void vtkPLOT3DReader::ReadIntBlockV(char** buf, int n, int* block)
{
  memcpy(block, *buf, sizeof(int)*n);

  if (this->ByteOrder == FILE_LITTLE_ENDIAN)
    {
    vtkByteSwap::Swap4LERange(block, n);
    }
  else
    {
    vtkByteSwap::Swap4BERange(block, n);
    }
  *buf += sizeof(int);
}

void vtkPLOT3DReader::SkipByteCountV(char** buf)
{
  if (this->HasByteCount)
    {
    *buf += sizeof(int);
    }
}

Berk Geveci's avatar
Berk Geveci committed
332
int vtkPLOT3DReader::VerifySettings(char* buf, int vtkNotUsed(bufSize))
333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398
{
  int numGrid=0;

  if ( this->MultiGrid )
    {
    this->SkipByteCountV(&buf);
    this->ReadIntBlockV(&buf, 1, &numGrid);
    this->SkipByteCountV(&buf);
    }
  else
    {
    numGrid=1;
    }
  cout << "Num. grids: " << numGrid << endl;

  int retVal=1;

  long fileSize = 0;
  // Size of number of grids information.
  if ( this->MultiGrid )
    {
    fileSize += 4; // numGrids
    if (this->HasByteCount)
      {
      fileSize += 4*4; // byte counts for the header
      }
    }

  // Add the size of each grid.
  this->SkipByteCountV(&buf);
  for(int i=0; i<numGrid; i++)
    {
    int ni, nj, nk;
    this->ReadIntBlockV(&buf, 1, &ni);
    cout << "Grid " << i << " ni " << ni << endl;
    this->ReadIntBlockV(&buf, 1, &nj);
    cout << "Grid " << i << " nj " << nj << endl;
    if (!this->TwoDimensionalGeometry)
      {
      this->ReadIntBlockV(&buf, 1, &nk);
      cout << "Grid " << i << " nk " << nk << endl;
      }
    else
      {
      nk = 1;
      }
    fileSize += this->EstimateSize(ni, nj, nk);
    // If this number is larger than the file size, there
    // is something wrong.
    if ( fileSize > this->FileSize )
      {
      retVal = 0;
      break;
      }
    }
  this->SkipByteCountV(&buf);
  // If this number is different than the file size, there
  // is something wrong.
  if ( fileSize != this->FileSize )
    {
    retVal = 0;
    }

  return retVal;
}

Berk Geveci's avatar
Berk Geveci committed
399 400
// Read the header and return the number of grids.
int vtkPLOT3DReader::GetNumberOfOutputsInternal(FILE* xyzFp, int verify)
401
{
Berk Geveci's avatar
Berk Geveci committed
402
  int numGrid=0;
403
  int numOutputs;
Berk Geveci's avatar
Berk Geveci committed
404 405

  if ( this->MultiGrid )
406
    {
Berk Geveci's avatar
Berk Geveci committed
407 408 409
    this->SkipByteCount(xyzFp);
    this->ReadIntBlock(xyzFp, 1, &numGrid);
    this->SkipByteCount(xyzFp);
410 411 412
    }
  else
    {
Berk Geveci's avatar
Berk Geveci committed
413
    numGrid=1;
414
    }
Berk Geveci's avatar
Berk Geveci committed
415 416 417 418 419 420 421

  if (!verify)
    {
    // We were told not the verify the number of grid. Just return it.
    numOutputs = numGrid;
    }
  else
422
    {
Berk Geveci's avatar
Berk Geveci committed
423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448
    // We were told to make sure that the file can really contain
    // the number of grid in the header (we can only check this
    // if file is binary)
    int error=0;
    if ( this->BinaryFile )
      {
      // Store the beginning of first grid.
      long pos = ftell(xyzFp);

      long fileSize = 0;
      // Size of number of grids information.
      if ( this->MultiGrid )
        {
        fileSize += 4; // numGrids
        if (this->HasByteCount)
          {
          fileSize += 4*4; // byte counts for the header
          }
        }
      // Add the size of each grid.
      this->SkipByteCount(xyzFp);
      for(int i=0; i<numGrid; i++)
        {
        int ni, nj, nk;
        this->ReadIntBlock(xyzFp, 1, &ni);
        this->ReadIntBlock(xyzFp, 1, &nj);
449 450 451 452 453 454 455 456
        if (!this->TwoDimensionalGeometry)
          {
          this->ReadIntBlock(xyzFp, 1, &nk);
          }
        else
          {
          nk = 1;
          }
Berk Geveci's avatar
Berk Geveci committed
457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477
        fileSize += this->EstimateSize(ni, nj, nk);
        // If this number is larger than the file size, there
        // is something wrong.
        if ( fileSize > this->FileSize )
          {
          error = 1;
          break;
          }
        }
      this->SkipByteCount(xyzFp);
      // If this number is different than the file size, there
      // is something wrong.
      if ( fileSize != this->FileSize && !this->ForceRead)
        {
        this->SetErrorCode(vtkErrorCode::FileFormatError);
        error = 1;
        }

      fseek(xyzFp, pos, SEEK_SET);
      }
    else
478
      {
Berk Geveci's avatar
Berk Geveci committed
479 480 481 482
      if (numGrid == 0)
        {
        this->SetErrorCode(vtkErrorCode::FileFormatError);
        }
483 484
      }
    
Berk Geveci's avatar
Berk Geveci committed
485 486
    // Now set the number of outputs.
    if (!error && numGrid != 0)
487
      {
Berk Geveci's avatar
Berk Geveci committed
488 489 490 491 492 493 494 495 496 497 498 499 500 501 502
      if ( !this->DoNotReduceNumberOfOutputs || 
           numGrid > this->NumberOfOutputs )
        {
        this->SetNumberOfOutputs(numGrid);
        }
      for (int i=1; i<numGrid; i++)
        {
        if (!this->Outputs[i])
          {
          vtkStructuredGrid* sg = vtkStructuredGrid::New();
          this->SetNthOutput(i, sg);
          sg->Delete();
          }
        }
      numOutputs = numGrid;
503
      }
Berk Geveci's avatar
Berk Geveci committed
504
    else
505
      {
Berk Geveci's avatar
Berk Geveci committed
506
      numOutputs = 0;
507 508
      }
    }
Berk Geveci's avatar
Berk Geveci committed
509 510

  return numOutputs;
511 512
}

Berk Geveci's avatar
Berk Geveci committed
513
int vtkPLOT3DReader::ReadGeometryHeader(FILE* fp)
Will Schroeder's avatar
Will Schroeder committed
514
{
Berk Geveci's avatar
Berk Geveci committed
515 516
  int numGrid = this->GetNumberOfOutputsInternal(fp, 1);
  int i;
517
  vtkDebugMacro("Geometry number of grids: " << numGrid);
Berk Geveci's avatar
Berk Geveci committed
518
  if ( numGrid == 0 )
519
    {
Berk Geveci's avatar
Berk Geveci committed
520 521
    // Bad file, set all extents to invalid.
    for (i=0; i<this->NumberOfOutputs; i++)
Bill Lorensen's avatar
Bill Lorensen committed
522
      {
Berk Geveci's avatar
Berk Geveci committed
523
      this->GetOutput(i)->SetWholeExtent(0, -1, 0, -1, 0, -1);
Bill Lorensen's avatar
Bill Lorensen committed
524
      }
Berk Geveci's avatar
Berk Geveci committed
525
    return VTK_ERROR;
526
    }
Berk Geveci's avatar
Berk Geveci committed
527 528 529 530

  // Read and set extents of all outputs.
  this->SkipByteCount(fp);
  for(i=0; i<numGrid; i++)
531
    {
Berk Geveci's avatar
Berk Geveci committed
532 533 534
    int ni, nj, nk;
    this->ReadIntBlock(fp, 1, &ni);
    this->ReadIntBlock(fp, 1, &nj);
535 536 537 538 539 540 541 542 543 544
    if (!this->TwoDimensionalGeometry)
      {
      this->ReadIntBlock(fp, 1, &nk);
      }
    else
      {
      nk = 1;
      }
    vtkDebugMacro("Geometry, block " << i << " dimensions: "
                  << ni << " " << nj << " " << nk);
Berk Geveci's avatar
Berk Geveci committed
545
    this->GetOutput(i)->SetWholeExtent(0, ni-1, 0, nj-1, 0, nk-1);
546
    }
Berk Geveci's avatar
Berk Geveci committed
547
  this->SkipByteCount(fp);
548 549 550 551 552 553 554 555 556 557 558

  if ( !this->PointCache )
    {
    this->PointCache = new vtkFloatArray*[ this->NumberOfOutputs + 1 ];
    this->IBlankCache = new vtkUnsignedCharArray* [ this->NumberOfOutputs + 1 ];
    for ( int g=0; g < this->NumberOfOutputs+1; ++g )
      {
      this->PointCache[g] = 0;
      this->IBlankCache[g] = 0;
      }
    }
Berk Geveci's avatar
Berk Geveci committed
559 560
  return VTK_OK;
}
561

Berk Geveci's avatar
Berk Geveci committed
562 563 564
int vtkPLOT3DReader::ReadQHeader(FILE* fp)
{
  int numGrid = this->GetNumberOfOutputsInternal(fp, 0);
565
  vtkDebugMacro("Q number of grids: " << numGrid);
Berk Geveci's avatar
Berk Geveci committed
566
  if ( numGrid == 0 )
567
    {
Berk Geveci's avatar
Berk Geveci committed
568
    return VTK_ERROR;
569
    }
Berk Geveci's avatar
Berk Geveci committed
570 571 572

  this->SkipByteCount(fp);
  for(int i=0; i<numGrid; i++)
Will Schroeder's avatar
Will Schroeder committed
573
    {
Berk Geveci's avatar
Berk Geveci committed
574 575 576 577
    int ni, nj, nk;
    this->ReadIntBlock(fp, 1, &ni);
    this->ReadIntBlock(fp, 1, &nj);
    this->ReadIntBlock(fp, 1, &nk);
578 579
    vtkDebugMacro("Q, block " << i << " dimensions: "
                  << ni << " " << nj << " " << nk);
580

Berk Geveci's avatar
Berk Geveci committed
581 582 583
    int extent[6];
    this->GetOutput(i)->GetWholeExtent(extent);
    if ( extent[1] != ni-1 || extent[3] != nj-1 || extent[5] != nk-1)
Ken Martin's avatar
Ken Martin committed
584
      {
Berk Geveci's avatar
Berk Geveci committed
585 586 587 588
      this->SetErrorCode(vtkErrorCode::FileFormatError);
      vtkErrorMacro("Geometry and data dimensions do not match. "
                    "Data file may be corrupt.");
      return VTK_ERROR;
589
      }
Will Schroeder's avatar
Will Schroeder committed
590
    }
Berk Geveci's avatar
Berk Geveci committed
591 592 593
  this->SkipByteCount(fp);
  return VTK_OK;
}
Will Schroeder's avatar
Will Schroeder committed
594

595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620
void vtkPLOT3DReader::SetXYZFileName( const char* name )
{
  if ( this->XYZFileName && ! strcmp( this->XYZFileName, name ) )
    {
    return;
    }

  if ( this->XYZFileName )
    {
    delete [] this->XYZFileName;
    }

  if ( name )
    {
    this->XYZFileName = new char [ strlen( name ) + 1 ];
    strcpy( this->XYZFileName, name );
    }
  else
    {
    this->XYZFileName = 0;
    }

  this->ClearGeometryCache();
  this->Modified();
}

Berk Geveci's avatar
Berk Geveci committed
621 622 623
void vtkPLOT3DReader::SetScalarFunctionNumber(int num)
{
  if ( this->ScalarFunctionNumber == num)
Will Schroeder's avatar
Will Schroeder committed
624
    {
Berk Geveci's avatar
Berk Geveci committed
625
    return;
626
    }
Berk Geveci's avatar
Berk Geveci committed
627
  if (num >= 0)
628
    {
Berk Geveci's avatar
Berk Geveci committed
629 630 631
    // If this function is not in the list, add it.
    int found=0;
    for (int i=0; i < this->FunctionList->GetNumberOfTuples(); i++ )
Bill Lorensen's avatar
Bill Lorensen committed
632
      {
Berk Geveci's avatar
Berk Geveci committed
633 634
      if ( this->FunctionList->GetValue(i) == num )
        {
635
        found=1;
Berk Geveci's avatar
Berk Geveci committed
636 637 638 639 640
        }
      }
    if (!found)
      {
      this->AddFunction(num);
Bill Lorensen's avatar
Bill Lorensen committed
641
      }
642
    }
Berk Geveci's avatar
Berk Geveci committed
643 644
  this->ScalarFunctionNumber = num;
}
645

Berk Geveci's avatar
Berk Geveci committed
646 647 648
void vtkPLOT3DReader::SetVectorFunctionNumber(int num)
{
  if ( this->VectorFunctionNumber == num)
649
    {
Berk Geveci's avatar
Berk Geveci committed
650
    return;
Will Schroeder's avatar
Will Schroeder committed
651
    }
Berk Geveci's avatar
Berk Geveci committed
652
  if (num >= 0)
653
    {
Berk Geveci's avatar
Berk Geveci committed
654 655 656
    // If this function is not in the list, add it.
    int found=0;
    for (int i=0; i < this->FunctionList->GetNumberOfTuples(); i++ )
657
      {
Berk Geveci's avatar
Berk Geveci committed
658 659
      if ( this->FunctionList->GetValue(i) == num )
        {
660
        found=1;
Berk Geveci's avatar
Berk Geveci committed
661 662 663 664 665
        }
      }
    if (!found)
      {
      this->AddFunction(num);
666 667
      }
    }
Berk Geveci's avatar
Berk Geveci committed
668 669
  this->VectorFunctionNumber = num;
}
670

Berk Geveci's avatar
Berk Geveci committed
671 672 673
void vtkPLOT3DReader::RemoveFunction(int fnum)
{
  for (int i=0; i < this->FunctionList->GetNumberOfTuples(); i++ )
674
    {
Berk Geveci's avatar
Berk Geveci committed
675
    if ( this->FunctionList->GetValue(i) == fnum )
Bill Lorensen's avatar
Bill Lorensen committed
676
      {
Berk Geveci's avatar
Berk Geveci committed
677 678
      this->FunctionList->SetValue(i,-1);
      this->Modified();
Bill Lorensen's avatar
Bill Lorensen committed
679
      }
680
    }
Berk Geveci's avatar
Berk Geveci committed
681
}
682

Berk Geveci's avatar
Berk Geveci committed
683 684 685
void vtkPLOT3DReader::ExecuteInformation()
{
  FILE* xyzFp;
686

Berk Geveci's avatar
Berk Geveci committed
687 688 689 690
  if ( this->CheckGeometryFile(xyzFp) != VTK_OK)
    {
    return;
    }
691

Berk Geveci's avatar
Berk Geveci committed
692 693
  this->CalculateFileSize(xyzFp);
  this->ReadGeometryHeader(xyzFp);
694

Berk Geveci's avatar
Berk Geveci committed
695
  fclose(xyzFp);
Will Schroeder's avatar
Will Schroeder committed
696 697
}

Berk Geveci's avatar
Berk Geveci committed
698
void vtkPLOT3DReader::Execute()
699
{
Berk Geveci's avatar
Berk Geveci committed
700 701 702 703
  this->SetErrorCode(vtkErrorCode::NoError);

  int i;
  int ndim, nx, ny, nz;
704
  int numberOfDims;
Berk Geveci's avatar
Berk Geveci committed
705 706
  vtkIdType index;

707 708
  // Don't read the geometry if we already have it!
  if ( (!this->PointCache) || (!this->PointCache[0]) )
709
    {
710 711 712 713 714
    FILE* xyzFp;
    if ( this->CheckGeometryFile(xyzFp) != VTK_OK)
      {
      return;
      }
715

716 717 718 719 720 721
    if ( this->ReadGeometryHeader(xyzFp) != VTK_OK )
      {
      vtkErrorMacro("Error reading geometry file.");
      fclose(xyzFp);
      return;
      }
722

723
    if (!this->TwoDimensionalGeometry)
Bill Lorensen's avatar
Bill Lorensen committed
724
      {
725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755
      numberOfDims = 3;
      }
    else
      {
      numberOfDims = 2;
      }
  
    for(i=0; i<this->NumberOfOutputs; i++)
      {

      // Read the geometry of this grid.
      this->SkipByteCount(xyzFp);

      vtkStructuredGrid* nthOutput = this->GetOutput(i);
      int dims[6];
      nthOutput->GetWholeExtent(dims);
      nthOutput->SetExtent(dims);
      nthOutput->GetDimensions(dims);
      this->PointCache[i] = vtkFloatArray::New();
      this->PointCache[i]->SetNumberOfComponents(3);
      this->PointCache[i]->SetNumberOfTuples( dims[0]*dims[1]*dims[2] );

      vtkPoints* points = vtkPoints::New();
      points->SetData(this->PointCache[i]);
      nthOutput->SetPoints(points);
      points->Delete();
      this->PointCache[i]->Register( this );
      this->PointCache[i]->Delete();

      float coord;
      for(ndim=0; ndim < numberOfDims; ndim++)
Berk Geveci's avatar
Berk Geveci committed
756
        {
757
        for(nz=0; nz < dims[2]; nz++)
Berk Geveci's avatar
Berk Geveci committed
758
          {
759
          for(ny=0; ny < dims[1]; ny++)
Berk Geveci's avatar
Berk Geveci committed
760
            {
761
            for(nx=0; nx < dims[0]; nx++)
Berk Geveci's avatar
Berk Geveci committed
762
              {
763
              if ( this->ReadFloatBlock(xyzFp, 1, &coord) == 0 )
Berk Geveci's avatar
Berk Geveci committed
764
                {
765 766 767 768 769
                vtkErrorMacro("Encountered premature end-of-file while reading "
                              "the geometry file (or the file is corrupt).");
                this->SetErrorCode(vtkErrorCode::PrematureEndOfFileError);
                // We need to generate output (otherwise, this filter will
                // keep executing). So we produce all 0's
Ken Martin's avatar
Ken Martin committed
770
                double nullpt[3] = {0.0, 0.0, 0.0};
771 772 773 774 775 776 777
                vtkIdType ipts, npts=this->PointCache[i]->GetNumberOfTuples();
                for(ipts=0; ipts < npts; ipts++)
                  {
                  this->PointCache[i]->SetTuple(ipts, nullpt);
                  }
                fclose(xyzFp);
                return;
Berk Geveci's avatar
Berk Geveci committed
778
                }
779 780
              index = nz*dims[0]*dims[1]+ny*dims[0]+nx;
              this->PointCache[i]->SetComponent(index, ndim, coord);
Berk Geveci's avatar
Berk Geveci committed
781 782 783 784
              }
            }
          }
        }
785

786
      if (this->TwoDimensionalGeometry)
787
        {
788 789 790 791 792
        vtkIdType ipts, npts=this->PointCache[i]->GetNumberOfTuples();
        for(ipts=0; ipts < npts; ipts++)
          {
          this->PointCache[i]->SetComponent(ipts, 2, 0);
          }
793 794
        }

795
      if (this->IBlanking)
Berk Geveci's avatar
Berk Geveci committed
796
        {
797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816
        this->IBlankCache[i] = vtkUnsignedCharArray::New();
        this->IBlankCache[i]->SetNumberOfComponents(1);
        this->IBlankCache[i]->SetNumberOfTuples( dims[0]*dims[1]*dims[2] );
        this->IBlankCache[i]->SetName("Visibility");
        int* ib = new int[dims[0]*dims[1]*dims[2]];
        if ( this->ReadIntBlock(xyzFp, dims[0]*dims[1]*dims[2], ib) == 0)
          {
          vtkErrorMacro("Encountered premature end-of-file while reading "
                        "the q file (or the file is corrupt).");
          this->SetErrorCode(vtkErrorCode::PrematureEndOfFileError);
          fclose(xyzFp);
          return;
          }
        vtkIdType ipts, npts=this->IBlankCache[i]->GetNumberOfTuples();
        unsigned char* ib2 = this->IBlankCache[i]->GetPointer(0);
        for (ipts=0; ipts<npts; ipts++)
          {
          ib2[ipts] = ib[ipts];
          }
        delete[] ib;
817
        nthOutput->SetPointVisibilityArray(this->IBlankCache[i]);
818 819
        this->IBlankCache[i]->Register( this );
        this->IBlankCache[i]->Delete();
Berk Geveci's avatar
Berk Geveci committed
820
        }
821
      this->SkipByteCount(xyzFp);
Berk Geveci's avatar
Berk Geveci committed
822 823
      }

824
    fclose(xyzFp);
825
    }
826 827 828
  else
    {
    numberOfDims = this->TwoDimensionalGeometry ? 2 : 3;
829

830 831 832 833 834 835 836 837 838 839 840 841 842 843
    for(i=0; i<this->NumberOfOutputs; i++)
      {
      vtkStructuredGrid* nthOutput = this->GetOutput(i);
      int dims[6];
      nthOutput->GetWholeExtent(dims);
      nthOutput->SetExtent(dims);

      vtkPoints* points = vtkPoints::New();
      points->SetData(this->PointCache[i]);
      nthOutput->SetPoints(points);
      points->Delete();

      if (this->IBlanking)
        {
844
        nthOutput->SetPointVisibilityArray(this->IBlankCache[i]);
845 846
        }
      }
Berk Geveci's avatar
Berk Geveci committed
847
    }
Berk Geveci's avatar
Berk Geveci committed
848 849 850

  // Now read the solution.
  if (this->QFileName && this->QFileName[0] != '\0')
851
    {
Berk Geveci's avatar
Berk Geveci committed
852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904
    FILE* qFp;
    if ( this->CheckSolutionFile(qFp) != VTK_OK)
      {
      return;
      }
    
    if ( this->ReadQHeader(qFp) != VTK_OK )
      {
      fclose(qFp);
      return;
      }

    for(i=0; i<this->NumberOfOutputs; i++)
      {
      vtkStructuredGrid* nthOutput = this->GetOutput(i);

      float fsmach, alpha, re, time;

      this->SkipByteCount(qFp);
      this->ReadFloatBlock(qFp, 1, &fsmach);
      this->ReadFloatBlock(qFp, 1, &alpha);
      this->ReadFloatBlock(qFp, 1, &re);
      this->ReadFloatBlock(qFp, 1, &time);
      this->SkipByteCount(qFp);

      // Save the properties first
      vtkFloatArray* properties = vtkFloatArray::New();
      properties->SetName("Properties");
      properties->SetNumberOfTuples(4);
      properties->SetTuple1(0, fsmach); 
      properties->SetTuple1(1, alpha); 
      properties->SetTuple1(2, re); 
      properties->SetTuple1(3, time); 
      nthOutput->GetFieldData()->AddArray(properties);
      properties->Delete();
      
      int dims[6];
      nthOutput->GetWholeExtent(dims);
      nthOutput->SetExtent(dims);
      nthOutput->GetDimensions(dims);

      this->SkipByteCount(qFp);

      vtkFloatArray* density = vtkFloatArray::New();
      density->SetNumberOfComponents(1);
      density->SetNumberOfTuples( dims[0]*dims[1]*dims[2] );
      density->SetName("Density");
      float* dens = density->GetPointer(0);
      if ( this->ReadFloatBlock(qFp, dims[0]*dims[1]*dims[2], dens) == 0)
        {
        vtkErrorMacro("Encountered premature end-of-file while reading "
                      "the q file (or the file is corrupt).");
        this->SetErrorCode(vtkErrorCode::PrematureEndOfFileError);
Berk Geveci's avatar
Berk Geveci committed
905
        fclose(qFp);
Berk Geveci's avatar
Berk Geveci committed
906 907 908 909 910 911 912 913 914 915 916
        return;
        }
      nthOutput->GetPointData()->AddArray(density);
      density->Delete();

      vtkFloatArray* momentum = vtkFloatArray::New();
      momentum->SetNumberOfComponents(3);
      momentum->SetNumberOfTuples( dims[0]*dims[1]*dims[2] );
      momentum->SetName("Momentum");

      float comp;
917
      for(ndim=0; ndim < numberOfDims; ndim++)
Berk Geveci's avatar
Berk Geveci committed
918 919 920 921 922 923 924 925 926 927 928
        {
        for(nz=0; nz < dims[2]; nz++)
          {
          for(ny=0; ny < dims[1]; ny++)
            {
            for(nx=0; nx < dims[0]; nx++)
              {
              if ( this->ReadFloatBlock(qFp, 1, &comp) == 0 )
                {
                vtkErrorMacro("Encountered premature end-of-file while "
                              "reading the q file (or the file is corrupt).");
Berk Geveci's avatar
Berk Geveci committed
929
                fclose(qFp);
Berk Geveci's avatar
Berk Geveci committed
930 931 932 933 934 935 936 937
                return;
                }
              index = nz*dims[0]*dims[1]+ny*dims[0]+nx;
              momentum->SetComponent(index, ndim, comp);
              }
            }
          }
        }
938 939 940 941 942 943 944 945
      if (this->TwoDimensionalGeometry)
        {
        vtkIdType ipts, npts=momentum->GetNumberOfTuples();
        for(ipts=0; ipts < npts; ipts++)
          {
          momentum->SetComponent(ipts, 2, 0);
          }
        }
Berk Geveci's avatar
Berk Geveci committed
946 947 948 949 950 951 952 953 954 955 956 957 958

      nthOutput->GetPointData()->AddArray(momentum);
      momentum->Delete();

      vtkFloatArray* se = vtkFloatArray::New();
      se->SetNumberOfComponents(1);
      se->SetNumberOfTuples( dims[0]*dims[1]*dims[2] );
      se->SetName("StagnationEnergy");
      float* sen = se->GetPointer(0);
      if (this->ReadFloatBlock(qFp, dims[0]*dims[1]*dims[2], sen) == 0)
        {
        vtkErrorMacro("Encountered premature end-of-file while reading "
                      "the q file (or the file is corrupt).");
Berk Geveci's avatar
Berk Geveci committed
959
        fclose(qFp);
Berk Geveci's avatar
Berk Geveci committed
960 961 962 963 964 965 966 967 968
        return;
        }
      nthOutput->GetPointData()->AddArray(se);
      se->Delete();

      
      if ( this->FunctionList->GetNumberOfTuples() > 0 )
        {
        int fnum;
969
        for (int tup=0; tup < this->FunctionList->GetNumberOfTuples(); tup++)
Berk Geveci's avatar
Berk Geveci committed
970
          {
971
          if ( (fnum=this->FunctionList->GetValue(tup)) >= 0 )
Berk Geveci's avatar
Berk Geveci committed
972 973 974 975 976 977 978 979 980 981
            {
            this->MapFunction(fnum, nthOutput);
            }
          }
        }
      this->AssignAttribute(this->ScalarFunctionNumber, nthOutput,
                            vtkDataSetAttributes::SCALARS);
      this->AssignAttribute(this->VectorFunctionNumber, nthOutput,
                            vtkDataSetAttributes::VECTORS);
      }
Berk Geveci's avatar
Berk Geveci committed
982
    fclose(qFp);
983 984 985 986
    }

}

987
// Various PLOT3D functions.....................
Berk Geveci's avatar
Berk Geveci committed
988
void vtkPLOT3DReader::MapFunction(int fNumber, vtkStructuredGrid* output)
Will Schroeder's avatar
Will Schroeder committed
989
{
990 991 992 993 994 995
  switch (fNumber)
    {
    case 100: //Density
      break;

    case 110: //Pressure
Berk Geveci's avatar
Berk Geveci committed
996
      this->ComputePressure(output);
997 998 999
      break;

    case 120: //Temperature
Berk Geveci's avatar
Berk Geveci committed
1000
      this->ComputeTemperature(output);
1001 1002 1003
      break;

    case 130: //Enthalpy
Berk Geveci's avatar
Berk Geveci committed
1004
      this->ComputeEnthalpy(output);
1005 1006 1007 1008
      break;

    case 140: //Internal Energy
      break;
Will Schroeder's avatar
Will Schroeder committed
1009

1010
    case 144: //Kinetic Energy
Berk Geveci's avatar
Berk Geveci committed
1011
      this->ComputeKineticEnergy(output);
1012
      break;
Will Schroeder's avatar
Will Schroeder committed
1013

1014
    case 153: //Velocity Magnitude
Berk Geveci's avatar
Berk Geveci committed
1015
      this->ComputeVelocityMagnitude(output);
1016
      break;
Will Schroeder's avatar
Will Schroeder committed
1017

1018 1019
    case 163: //Stagnation energy
      break;
Will Schroeder's avatar
Will Schroeder committed
1020

1021
    case 170: //Entropy
Berk Geveci's avatar
Berk Geveci committed
1022
      this->ComputeEntropy(output);
1023 1024 1025
      break;

    case 184: //Swirl
Berk Geveci's avatar
Berk Geveci committed
1026
      this->ComputeSwirl(output);
1027 1028 1029
      break;

    case 200: //Velocity
Berk Geveci's avatar
Berk Geveci committed
1030
      this->ComputeVelocity(output);
1031 1032 1033
      break;

    case 201: //Vorticity
Berk Geveci's avatar
Berk Geveci committed
1034
      this->ComputeVorticity(output);
1035 1036 1037 1038 1039 1040
      break;

    case 202: //Momentum
      break;

    case 210: //PressureGradient
Berk Geveci's avatar
Berk Geveci committed
1041
      this->ComputePressureGradient(output);
1042 1043 1044
      break;

    default:
Ken Martin's avatar
Ken Martin committed
1045
      vtkErrorMacro(<<"No function number " << fNumber);
1046 1047 1048
    }
}

Berk Geveci's avatar
Berk Geveci committed
1049 1050
void vtkPLOT3DReader::AssignAttribute(int fNumber, vtkStructuredGrid* output,
                                  int attributeType)
1051
{
Berk Geveci's avatar
Berk Geveci committed
1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111 1112 1113 1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131
  switch (fNumber)
    {
    case -1:  //empty mapping
      output->GetPointData()->SetActiveAttribute(0, 
                                                 attributeType);
      break;

    case 100: //Density
      output->GetPointData()->SetActiveAttribute("Density", 
                                                 attributeType);
      break;

    case 110: //Pressure
      output->GetPointData()->SetActiveAttribute("Pressure", 
                                                 attributeType);
      break;

    case 120: //Temperature
      output->GetPointData()->SetActiveAttribute("Temperature", 
                                                 attributeType);
      break;

    case 130: //Enthalpy
      output->GetPointData()->SetActiveAttribute("Enthalpy", 
                                                 attributeType);
      break;

    case 140: //Internal Energy
      output->GetPointData()->SetActiveAttribute("StagnationEnergy", 
                                                 attributeType);
      break;

    case 144: //Kinetic Energy
      output->GetPointData()->SetActiveAttribute("KineticEnergy", 
                                                 attributeType);
      break;

    case 153: //Velocity Magnitude
      output->GetPointData()->SetActiveAttribute("VelocityMagnitude", 
                                                 attributeType);
      break;

    case 163: //Stagnation energy
      output->GetPointData()->SetActiveAttribute("StagnationEnergy", 
                                                 attributeType);
      break;

    case 170: //Entropy
      output->GetPointData()->SetActiveAttribute("Entropy", 
                                                 attributeType);
      break;

    case 184: //Swirl
      output->GetPointData()->SetActiveAttribute("Swirl", 
                                                 attributeType);
      break;

    case 200: //Velocity
      output->GetPointData()->SetActiveAttribute("Velocity", 
                                                 attributeType);
      break;

    case 201: //Vorticity
      output->GetPointData()->SetActiveAttribute("Vorticity", 
                                                 attributeType);
      break;

    case 202: //Momentum
      output->GetPointData()->SetActiveAttribute("Momentum", 
                                                 attributeType);
      break;

    case 210: //PressureGradient
      output->GetPointData()->SetActiveAttribute("PressureGradient", 
                                                 attributeType);
      break;

    default:
      vtkErrorMacro(<<"No function number " << fNumber);
    }
Will Schroeder's avatar
Will Schroeder committed
1132 1133
}

Berk Geveci's avatar
Berk Geveci committed
1134
void vtkPLOT3DReader::ComputeTemperature(vtkStructuredGrid* output)
Will Schroeder's avatar
Will Schroeder committed
1135
{
Ken Martin's avatar
Ken Martin committed
1136
  double *m, e, rr, u, v, w, v2, p, d, rrgas;
Berk Geveci's avatar
Berk Geveci committed
1137
  vtkIdType i;
1138
  vtkFloatArray *temperature;
1139 1140 1141

  //  Check that the required data is available
  //
Berk Geveci's avatar
Berk Geveci committed
1142 1143 1144 1145 1146 1147 1148
  vtkPointData* outputPD = output->GetPointData();
  vtkDataArray* density = outputPD->GetArray("Density");
  vtkDataArray* momentum = outputPD->GetArray("Momentum");
  vtkDataArray* energy = outputPD->GetArray("StagnationEnergy");
  
  if ( density == NULL || momentum == NULL || 
       energy == NULL )
Will Schroeder's avatar
Will Schroeder committed
1149
    {
Ken Martin's avatar
Ken Martin committed
1150
    vtkErrorMacro(<<"Cannot compute temperature");
Will Schroeder's avatar
Will Schroeder committed
1151 1152 1153
    return;
    }

Berk Geveci's avatar
Berk Geveci committed
1154
  vtkIdType numPts = density->GetNumberOfTuples();
1155
  temperature = vtkFloatArray::New();
Berk Geveci's avatar
Berk Geveci committed
1156
  temperature->SetNumberOfTuples(numPts);
1157 1158 1159

  //  Compute the temperature
  //
Will Schroeder's avatar
Will Schroeder committed
1160
  rrgas = 1.0 / this->R;
Berk Geveci's avatar
Berk Geveci committed
1161
  for (i=0; i < numPts; i++) 
Will Schroeder's avatar
Will Schroeder committed
1162
    {
Berk Geveci's avatar
Berk Geveci committed
1163
    d = density->GetComponent(i,0);
Will Schroeder's avatar
Will Schroeder committed
1164
    d = (d != 0.0 ? d : 1.0);
Berk Geveci's avatar
Berk Geveci committed
1165 1166
    m = momentum->GetTuple(i);
    e = energy->GetComponent(i,0);
Will Schroeder's avatar
Will Schroeder committed
1167 1168 1169 1170 1171 1172
    rr = 1.0 / d;
    u = m[0] * rr;        
    v = m[1] * rr;        
    w = m[2] * rr;        
    v2 = u*u + v*v + w*w;
    p = (this->Gamma-1.) * (e - 0.5 * d * v2);
1173
    temperature->SetValue(i, p*rr*rrgas);
Will Schroeder's avatar
Will Schroeder committed
1174
  }
1175

1176 1177
  temperature->SetName("Temperature");
  outputPD->AddArray(temperature);
1178
  
Ken Martin's avatar
Ken Martin committed
1179
  temperature->Delete();
Ken Martin's avatar
Ken Martin committed
1180
  vtkDebugMacro(<<"Created temperature scalar");
Will Schroeder's avatar
Will Schroeder committed
1181 1182
}

Berk Geveci's avatar
Berk Geveci committed
1183
void vtkPLOT3DReader::ComputePressure(vtkStructuredGrid* output)
Will Schroeder's avatar
Will Schroeder committed
1184
{
Ken Martin's avatar
Ken Martin committed
1185
  double *m, e, u, v, w, v2, p, d, rr;
Berk Geveci's avatar
Berk Geveci committed
1186
  vtkIdType i;
1187
  vtkFloatArray *pressure;
1188 1189 1190

  //  Check that the required data is available
  //
Berk Geveci's avatar
Berk Geveci committed
1191 1192 1193 1194 1195 1196
  vtkPointData* outputPD = output->GetPointData();
  vtkDataArray* density = outputPD->GetArray("Density");
  vtkDataArray* momentum = outputPD->GetArray("Momentum");
  vtkDataArray* energy = outputPD->GetArray("StagnationEnergy");
  if ( density == NULL || momentum == NULL || 
       energy == NULL )
1197
    {
Ken Martin's avatar
Ken Martin committed
1198
    vtkErrorMacro(<<"Cannot compute pressure");
1199 1200 1201
    return;
    }

Berk Geveci's avatar
Berk Geveci committed
1202
  vtkIdType numPts = density->GetNumberOfTuples();
1203
  pressure = vtkFloatArray::New();
Berk Geveci's avatar
Berk Geveci committed
1204
  pressure->SetNumberOfTuples(numPts);
1205 1206 1207

  //  Compute the pressure
  //
Berk Geveci's avatar
Berk Geveci committed
1208
  for (i=0; i < numPts; i++) 
1209
    {
Berk Geveci's avatar
Berk Geveci committed
1210
    d = density->GetComponent(i,0);
1211
    d = (d != 0.0 ? d : 1.0);
Berk Geveci's avatar
Berk Geveci committed
1212 1213
    m = momentum->GetTuple(i);
    e = energy->GetComponent(i,0);