vtkAMREnzoReader.cxx 13.2 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
/*=========================================================================

 Program:   Visualization Toolkit
 Module:    vtkAMREnzoReader.cxx

 Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
 All rights reserved.
 See Copyright.txt or http://www.kitware.com/Copyright.htm for details.

 This software is distributed WITHOUT ANY WARRANTY; without even
 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
 PURPOSE.  See the above copyright notice for more information.

 =========================================================================*/
15
16
17

#include "vtkAMREnzoReader.h"
#include "vtkObjectFactory.h"
18
#include "vtkOverlappingAMR.h"
19
20
21
22
23
24
25
26
#include "vtkUniformGrid.h"
#include "vtkDataArraySelection.h"
#include "vtkDataArray.h"
#include "vtkPolyData.h"
#include "vtkAMRUtilities.h"
#include "vtkIndent.h"
#include "vtksys/SystemTools.hxx"

27
28
29
30
31
32
33
34
35
36
37
38
#include "vtkDataSet.h"
#include "vtkCellData.h"
#include "vtkIntArray.h"
#include "vtkLongArray.h"
#include "vtkShortArray.h"
#include "vtkFloatArray.h"
#include "vtkDoubleArray.h"
#include "vtkLongLongArray.h"
#include "vtkUnsignedIntArray.h"
#include "vtkUnsignedCharArray.h"
#include "vtkUnsignedShortArray.h"

39
40
41
#define H5_USE_16_API
#include <hdf5.h>

George Zagaris's avatar
George Zagaris committed
42
#include <sstream>
43
44
#include <vector>
#include <string>
45
46
#include <cassert>

47
#include "vtkAMREnzoReaderInternal.h"
48
49
50
51

vtkStandardNewMacro(vtkAMREnzoReader);

//-----------------------------------------------------------------------------
52
53
vtkAMREnzoReader::vtkAMREnzoReader()
{
54
  this->Internal = new vtkEnzoReaderInternal();
55
56
57
  this->Initialize();
}

58
//-----------------------------------------------------------------------------
59
60
vtkAMREnzoReader::~vtkAMREnzoReader()
{
61
62
  delete this->Internal;
  this->Internal=NULL;
George Zagaris's avatar
George Zagaris committed
63
64
65
66
67

  this->BlockMap.clear();

  if( this->FileName )
    {
68
69
    delete [] this->FileName;
    this->FileName = NULL;
George Zagaris's avatar
George Zagaris committed
70
    }
71
72
73
74
75
76
}

//-----------------------------------------------------------------------------
void vtkAMREnzoReader::PrintSelf( std::ostream &os, vtkIndent indent )
{
  this->Superclass::PrintSelf( os, indent );
77
78
}

George Zagaris's avatar
George Zagaris committed
79
80
81
82
83
84
85
86
87
88
89
90
//-----------------------------------------------------------------------------
int vtkAMREnzoReader::GetIndexFromArrayName( std::string arrayName )
{
  char stringIdx = arrayName.at( arrayName.size()-2 );
  return( atoi( &stringIdx ) );
}

//-----------------------------------------------------------------------------
double vtkAMREnzoReader::GetConversionFactor( const std::string name )
{
  if( this->label2idx.find( name ) != this->label2idx.end() )
    {
91
92
93
94
95
96
97
98
99
    int idx = this->label2idx[ name ];
    if( this->conversionFactors.find( idx ) != this->conversionFactors.end() )
      {
      return( this->conversionFactors[ idx ] );
      }
    else
      {
      return( 1.0 );
      }
George Zagaris's avatar
George Zagaris committed
100
    }
101
  return( 1.0 );
George Zagaris's avatar
George Zagaris committed
102
103
104
105
106
107
108
109
110
111
112
113
114
}

//-----------------------------------------------------------------------------
void vtkAMREnzoReader::ParseLabel(
    const std::string labelString, int &idx, std::string &label)
{

  std::vector< std::string > strings;

  std::istringstream iss( labelString );
  std::string word;
  while ( iss >> word )
    {
115
116
117
118
    if( !vtksys::SystemTools::StringStartsWith( word.c_str(), "=") )
      {
      strings.push_back( word );
      }
George Zagaris's avatar
George Zagaris committed
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
    }

  idx   = this->GetIndexFromArrayName( strings[0] );
  label = strings[ strings.size()-1 ];
}

//-----------------------------------------------------------------------------
void vtkAMREnzoReader::ParseCFactor(
    const std::string labelString, int &idx, double &factor )
{
  std::vector< std::string > strings;

  std::istringstream iss( labelString );
  std::string word;
  while( iss >> word )
    {
135
136
137
138
    if( ! vtksys::SystemTools::StringStartsWith( word.c_str(), "=") )
      {
      strings.push_back( word );
      }
George Zagaris's avatar
George Zagaris committed
139
140
141
142
143
144
145
146
147
148
149
150
151
152
    }

  idx    = this->GetIndexFromArrayName( strings[0] );
  factor = atof( strings[ strings.size()-1 ].c_str() );

}

//-----------------------------------------------------------------------------
void vtkAMREnzoReader::ParseConversionFactors( )
{
  assert( "pre: FileName should not be NULL" && (this->FileName != NULL) );

  // STEP 0: Extract the parameters file from the user-supplied filename
  std::string baseDir =
153
      vtksys::SystemTools::GetFilenamePath( std::string( this->FileName ) );
George Zagaris's avatar
George Zagaris committed
154
155

  std::string paramsFile = baseDir + "/" +
156
157
   vtksys::SystemTools::GetFilenameWithoutExtension(
    std::string( this->FileName ) );
George Zagaris's avatar
George Zagaris committed
158
159
160
161

  // STEP 1: Open Parameters file
  std::ifstream ifs;
  ifs.open( paramsFile.c_str() );
162
163
  if( !ifs.is_open( ) )
    {
164
165
    vtkWarningMacro( "Cannot open ENZO parameters file!\n" );
    return;
166
    }
George Zagaris's avatar
George Zagaris committed
167
168
169
170
171
172
173
174

  // STEP 2: Parsing parameters file
  std::string line;  // temp string to store a line read from the params file
  std::string label; // stores the attribute name
  double cf;         // stores the conversion factor
  int    idx;        // stores the attribute label index
  while( getline(ifs, line) )
    {
175
176
177
178
179
180
181
182
183
184
185
186
    if( vtksys::SystemTools::StringStartsWith(
        line.c_str(), "DataLabel" ) )
      {
      this->ParseLabel( line, idx, label );
      this->label2idx[ label ] = idx;
      }
    else if( vtksys::SystemTools::StringStartsWith(
             line.c_str(), "#DataCGSConversionFactor" ) )
      {
      this->ParseCFactor( line, idx, cf );
      this->conversionFactors[ idx ] = cf;
      }
George Zagaris's avatar
George Zagaris committed
187
188
189
190
191
192
    }

  // STEP 3: Close parameters file
  ifs.close();
}

193
194
//-----------------------------------------------------------------------------
void vtkAMREnzoReader::SetFileName( const char* fileName )
195
{
George Zagaris's avatar
George Zagaris committed
196
197
  assert("pre: Internal Enzo AMR Reader is NULL" && (this->Internal != NULL ));

198
199
200
201
202
  int isValid=0;

  if( fileName && strcmp( fileName, "" ) &&
    ( (this->FileName==NULL) || (strcmp(fileName,this->FileName ) ) ) )
    {
203
204
205
    std::string  tempName( fileName );
    std::string  bExtName( ".boundary" );
    std::string  hExtName( ".hierarchy" );
206

207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
    if( tempName.length() > hExtName.length() &&
        tempName.substr(tempName.length()-hExtName.length() )== hExtName )
      {
      this->Internal->MajorFileName =
          tempName.substr( 0, tempName.length() - hExtName.length() );
      this->Internal->HierarchyFileName = tempName;
      this->Internal->BoundaryFileName  = this->Internal->MajorFileName+bExtName;
      }
   else if( tempName.length() > bExtName.length() &&
       tempName.substr( tempName.length() - bExtName.length() )==bExtName )
     {
     this->Internal->MajorFileName =
        tempName.substr( 0, tempName.length() - bExtName.length() );
     this->Internal->BoundaryFileName  = tempName;
     this->Internal->HierarchyFileName = this->Internal->MajorFileName+hExtName;
     }
    else
      {
      vtkErrorMacro( "Enzo file has invalid extension!");
      return;
      }
    isValid = 1;
    this->Internal->DirectoryName =
     GetEnzoDirectory(this->Internal->MajorFileName.c_str());
231
232
233
234
    }

  if( isValid )
    {
235
236
237
238
    this->BlockMap.clear();
    this->Internal->Blocks.clear();
    this->Internal->NumberOfBlocks = 0;
    this->LoadedMetaData = false;
239

240
241
242
243
244
245
246
247
248
249
250
    if ( this->FileName != NULL)
    {
    delete [] this->FileName;
    this->FileName = NULL;
    this->Internal->SetFileName( NULL );
    }
    this->FileName = new char[  strlen( fileName ) + 1  ];
    strcpy( this->FileName, fileName );
    this->FileName[ strlen( fileName ) ] = '\0';
    this->Internal->SetFileName( this->FileName );
    this->ParseConversionFactors( );
251
252
253
    }

  this->Internal->ReadMetaData();
254
255
  this->SetUpDataArraySelections();
  this->InitializeArraySelections();
256
  this->Modified();
George Zagaris's avatar
George Zagaris committed
257
  return;
258
259
}

260
//-----------------------------------------------------------------------------
261
262
void vtkAMREnzoReader::ReadMetaData()
{
263
264
  assert( "pre: Internal Enzo Reader is NULL" && (this->Internal != NULL) );
  this->Internal->ReadMetaData();
265
266
}

267
//-----------------------------------------------------------------------------
268
269
void vtkAMREnzoReader::GenerateBlockMap()
{
270
271
272
273
  assert( "pre: Internal Enzo Reader is NULL" && (this->Internal != NULL) );

  this->BlockMap.clear();
  this->Internal->ReadMetaData();
George Zagaris's avatar
George Zagaris committed
274

275
276
  for( int i=0; i < this->Internal->NumberOfBlocks; ++i )
    {
277
278
279
280
    if( this->GetBlockLevel( i ) <= this->MaxLevel )
      {
      this->BlockMap.push_back( i );
      }
281
    } // END for all blocks
George Zagaris's avatar
George Zagaris committed
282
283
}

284
//-----------------------------------------------------------------------------
George Zagaris's avatar
George Zagaris committed
285
286
int vtkAMREnzoReader::GetBlockLevel( const int blockIdx )
{
287
288
289
  assert( "pre: Internal Enzo Reader is NULL" && (this->Internal != NULL) );

  this->Internal->ReadMetaData();
George Zagaris's avatar
George Zagaris committed
290

291
292
  if( blockIdx < 0 || blockIdx >= this->Internal->NumberOfBlocks )
    {
293
294
    vtkErrorMacro( "Block Index (" << blockIdx << ") is out-of-bounds!" );
    return( -1 );
295
296
    }
  return( this->Internal->Blocks[ blockIdx+1 ].Level );
297
298
}

299
//-----------------------------------------------------------------------------
300
301
int vtkAMREnzoReader::GetNumberOfBlocks()
{
302
303
304
  assert( "pre: Internal Enzo Reader is NULL" && (this->Internal != NULL) );
  this->Internal->ReadMetaData();
  return( this->Internal->NumberOfBlocks );
305
306
}

307
//-----------------------------------------------------------------------------
308
309
int vtkAMREnzoReader::GetNumberOfLevels()
{
310
311
312
  assert( "pre: Internal Enzo Reader is NULL" && (this->Internal != NULL) );
  this->Internal->ReadMetaData();
  return( this->Internal->NumberOfLevels );
313
314
}

315
//-----------------------------------------------------------------------------
316
int vtkAMREnzoReader::FillMetaData( )
317
318
{
  assert( "pre: Internal Enzo Reader is NULL" && (this->Internal != NULL) );
319
  assert( "pre: metadata object is NULL" && (this->Metadata != NULL) );
320
321
322

  this->Internal->ReadMetaData();
  std::vector< int > b2level;
323
  b2level.resize( this->Internal->NumberOfLevels+1, 0 );
324
325
326
327

  // this->Internal->Blocks includes a pseudo block -- the root as block #0
  for( int i=0; i < this->Internal->NumberOfBlocks; ++i )
    {
328
329
330
331
    vtkEnzoReaderBlock &theBlock = this->Internal->Blocks[ i+1 ];
    int level                    = theBlock.Level;
    int id                       = b2level[ level ];
    int internalIdx              = i;
332

333
334
335
    double blockMin[ 3 ];
    double blockMax[ 3 ];
    double spacings[ 3 ];
336

337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
    for( int j=0; j < 3; ++j )
      {
      blockMin[j] = theBlock.MinBounds[j];
      blockMax[j] = theBlock.MaxBounds[j];
      spacings[j] = (theBlock.BlockNodeDimensions[j] > 1)?
      (blockMax[j]-blockMin[j])/(theBlock.BlockNodeDimensions[j]-1.0):1.0;
      }

    vtkUniformGrid *ug = vtkUniformGrid::New();
    ug->SetDimensions( theBlock.BlockNodeDimensions );
    ug->SetOrigin( blockMin[0], blockMin[1], blockMin[2] );
    ug->SetSpacing( spacings[0], spacings[1], spacings[2] );

    this->Metadata->SetDataSet( level, id, ug );
    this->Metadata->SetCompositeIndex( level, id, internalIdx );
    ug->Delete();
    b2level[ level ]++;
354
355
356
    } // END for all blocks

  // NOTE: the controller here is null since each process loads its own metadata
357
  vtkAMRUtilities::GenerateMetaData( this->Metadata, NULL );
358
359
360
  return( 1 );
}

361
//-----------------------------------------------------------------------------
362
vtkUniformGrid* vtkAMREnzoReader::GetAMRGrid( const int blockIdx )
363
{
364
365
366
  assert( "pre: Internal Enzo Reader is NULL" && (this->Internal != NULL) );

  this->Internal->ReadMetaData();
George Zagaris's avatar
George Zagaris committed
367
368

  // this->Internal->Blocks includes a pseudo block --- the root as block #0
369
370
371
372
373
374
375
  vtkEnzoReaderBlock &theBlock = this->Internal->Blocks[ blockIdx+1 ];
  double blockMin[3];
  double blockMax[3];
  double spacings[3];

  for( int i=0; i < 3; ++i )
    {
376
377
378
379
    blockMin[ i ] = theBlock.MinBounds[ i ];
    blockMax[ i ] = theBlock.MaxBounds[ i ];
    spacings[ i ] = (theBlock.BlockNodeDimensions[i] > 1)?
        (blockMax[i]-blockMin[i])/(theBlock.BlockNodeDimensions[i]-1.0) : 1.0;
380
381
382
383
    }

  vtkUniformGrid *ug = vtkUniformGrid::New();
  ug->SetDimensions( theBlock.BlockNodeDimensions );
384
385
386
387
  ug->SetOrigin( blockMin[0], blockMin[1], blockMin[2] );
  ug->SetSpacing( spacings[0], spacings[1], spacings[2] );
  return( ug );
}
388

389
390
391
392
393
//-----------------------------------------------------------------------------
void vtkAMREnzoReader::GetAMRGridData(
    const int blockIdx, vtkUniformGrid *block, const char *field)
{
  assert( "pre: AMR block is NULL" && (block != NULL));
394

395
396
397
  this->Internal->GetBlockAttribute( field, blockIdx, block );
  if( this->ConvertToCGS == 1 )
    {
398
399
400
401
402
403
404
405
    double conversionFactor = this->GetConversionFactor(field);
    if( conversionFactor != 1.0 )
      {
      vtkDataArray *data = block->GetCellData()->GetArray( field );
      assert( "pre: data array is NULL!" && (data != NULL) );

      vtkIdType numTuples = data->GetNumberOfTuples();
      for( vtkIdType t=0; t < numTuples; ++t )
406
        {
407
408
409
410
411
412
413
414
        int numComp = data->GetNumberOfComponents();
        for( int c=0; c < numComp; ++c )
         {
         double f = data->GetComponent( t, c );
         data->SetComponent( t, c, f*conversionFactor );
         } // END for all components
        } // END for all tuples
      } // END if the conversion factor is not 1.0
415
    } // END if conversion to CGS units is requested
416
417
}

418

419
420
//-----------------------------------------------------------------------------
void vtkAMREnzoReader::SetUpDataArraySelections()
George Zagaris's avatar
George Zagaris committed
421
{
422
423
424
  assert( "pre: Internal Enzo Reader is NULL" && (this->Internal != NULL) );
  this->Internal->ReadMetaData();
  this->Internal->GetAttributeNames();
George Zagaris's avatar
George Zagaris committed
425

426
  int numAttrs = static_cast< int >(this->Internal->BlockAttributeNames.size());
427
428
  for( int i=0; i < numAttrs; i++ )
    {
429
430
    this->CellDataArraySelection->AddArray(
        this->Internal->BlockAttributeNames[i].c_str() );
431
    } // END for all attributes
432

George Zagaris's avatar
George Zagaris committed
433
}