vtkAMREnzoReader.cxx 13.8 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
#include "vtkUniformGrid.h"
#include "vtkDataArraySelection.h"
#include "vtkDataArray.h"
#include "vtkPolyData.h"
#include "vtkIndent.h"
24
#include "vtkInformation.h"
25 26
#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
#define H5_USE_16_API
40
#include "vtk_hdf5.h"
41

42
#include <sstream>
43 44
#include <vector>
#include <string>
45 46
#include <cassert>

47
#include "vtkAMREnzoReaderInternal.h"
48 49 50

vtkStandardNewMacro(vtkAMREnzoReader);

51 52 53 54 55 56 57 58 59
#include "vtkAMRInformation.h"
#include <limits>

void vtkAMREnzoReader::ComputeStats(vtkEnzoReaderInternal* internal, std::vector<int>& numBlocks, double min[3])
{
  min[0] = min[1] = min[2] = std::numeric_limits<double>::max();
  numBlocks.resize( this->Internal->NumberOfLevels, 0 );

  for( int i=0; i < internal->NumberOfBlocks; ++i )
60
  {
61 62 63
    vtkEnzoReaderBlock &theBlock = internal->Blocks[ i+1 ];
    double* gridMin = theBlock.MinBounds;
    if( gridMin[0] < min[0] )
64
    {
65
      min[0] = gridMin[0];
66
    }
67
    if( gridMin[1] < min[1] )
68
    {
69
      min[1] = gridMin[1];
70
    }
71
    if( gridMin[2] < min[2] )
72
    {
73 74
      min[2] = gridMin[2];
    }
75 76
    numBlocks[theBlock.Level]++;
  }
77 78
}

79
//-----------------------------------------------------------------------------
80 81
vtkAMREnzoReader::vtkAMREnzoReader()
{
82
  this->Internal = new vtkEnzoReaderInternal();
83
  this->IsReady  = false;
84
  this->Initialize();
85
  this->ConvertToCGS = 1;
86 87
}

88
//-----------------------------------------------------------------------------
89 90
vtkAMREnzoReader::~vtkAMREnzoReader()
{
91
  delete this->Internal;
92
  this->Internal=nullptr;
93 94 95

  this->BlockMap.clear();

96
  delete [] this->FileName;
97
  this->FileName = nullptr;
98 99 100 101 102 103
}

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

106 107 108
//-----------------------------------------------------------------------------
int vtkAMREnzoReader::GetIndexFromArrayName( std::string arrayName )
{
109 110 111 112
  char stringIdx[2];
  stringIdx[0] = arrayName.at( arrayName.size()-2 );
  stringIdx[1] = '\0';
  return( atoi( stringIdx ) );
113 114 115
}

//-----------------------------------------------------------------------------
116
double vtkAMREnzoReader::GetConversionFactor( const std::string &name )
117 118
{
  if( this->label2idx.find( name ) != this->label2idx.end() )
119
  {
120 121
    int idx = this->label2idx[ name ];
    if( this->conversionFactors.find( idx ) != this->conversionFactors.end() )
122
    {
123
      return( this->conversionFactors[ idx ] );
124
    }
125
    else
126
    {
127
      return( 1.0 );
128
    }
129
  }
130
  return( 1.0 );
131 132 133 134
}

//-----------------------------------------------------------------------------
void vtkAMREnzoReader::ParseLabel(
135
    const std::string &labelString, int &idx, std::string &label)
136 137 138 139 140 141 142
{

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

  std::istringstream iss( labelString );
  std::string word;
  while ( iss >> word )
143
  {
144
    if( !vtksys::SystemTools::StringStartsWith( word.c_str(), "=") )
145
    {
146
      strings.push_back( word );
147
    }
148
  }
149 150 151 152 153 154 155

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

//-----------------------------------------------------------------------------
void vtkAMREnzoReader::ParseCFactor(
156
    const std::string &labelString, int &idx, double &factor )
157 158 159 160 161 162
{
  std::vector< std::string > strings;

  std::istringstream iss( labelString );
  std::string word;
  while( iss >> word )
163
  {
164
    if( ! vtksys::SystemTools::StringStartsWith( word.c_str(), "=") )
165
    {
166
      strings.push_back( word );
167
    }
168
  }
169 170 171 172 173 174 175 176 177

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

}

//-----------------------------------------------------------------------------
void vtkAMREnzoReader::ParseConversionFactors( )
{
178
  assert( "pre: FileName should not be nullptr" && (this->FileName != nullptr) );
179 180 181

  // STEP 0: Extract the parameters file from the user-supplied filename
  std::string baseDir =
182
      vtksys::SystemTools::GetFilenamePath( std::string( this->FileName ) );
183 184

  std::string paramsFile = baseDir + "/" +
185 186
   vtksys::SystemTools::GetFilenameWithoutExtension(
    std::string( this->FileName ) );
187 188 189 190

  // STEP 1: Open Parameters file
  std::ifstream ifs;
  ifs.open( paramsFile.c_str() );
191
  if( !ifs.is_open( ) )
192
  {
193 194
    vtkWarningMacro( "Cannot open ENZO parameters file!\n" );
    return;
195
  }
196 197 198 199 200 201 202

  // 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) )
203
  {
204 205
    if( vtksys::SystemTools::StringStartsWith(
        line.c_str(), "DataLabel" ) )
206
    {
207 208
      this->ParseLabel( line, idx, label );
      this->label2idx[ label ] = idx;
209
    }
210 211
    else if( vtksys::SystemTools::StringStartsWith(
             line.c_str(), "#DataCGSConversionFactor" ) )
212
    {
213 214
      this->ParseCFactor( line, idx, cf );
      this->conversionFactors[ idx ] = cf;
215
    }
216
  }
217 218 219 220 221

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

222 223
//-----------------------------------------------------------------------------
void vtkAMREnzoReader::SetFileName( const char* fileName )
224
{
225
  assert("pre: Internal Enzo AMR Reader is nullptr" && (this->Internal != nullptr ));
226

227
  if( fileName && strcmp( fileName, "" ) &&
228
    ( (this->FileName==nullptr) || (strcmp(fileName,this->FileName ) ) ) )
229
  {
230 231 232
    std::string  tempName( fileName );
    std::string  bExtName( ".boundary" );
    std::string  hExtName( ".hierarchy" );
233

234 235
    if( tempName.length() > hExtName.length() &&
        tempName.substr(tempName.length()-hExtName.length() )== hExtName )
236
    {
237 238 239 240
      this->Internal->MajorFileName =
          tempName.substr( 0, tempName.length() - hExtName.length() );
      this->Internal->HierarchyFileName = tempName;
      this->Internal->BoundaryFileName  = this->Internal->MajorFileName+bExtName;
241
    }
242 243 244 245 246 247 248 249
    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;
    }
250
    else
251
    {
252 253
      vtkErrorMacro( "Enzo file has invalid extension!");
      return;
254
    }
255 256

    this->IsReady = true;
257
    this->Internal->DirectoryName =
258
        GetEnzoDirectory(this->Internal->MajorFileName.c_str());
259
  }
260

261
  if( this->IsReady )
262
  {
263 264 265 266
    this->BlockMap.clear();
    this->Internal->Blocks.clear();
    this->Internal->NumberOfBlocks = 0;
    this->LoadedMetaData = false;
267

268
    if ( this->FileName != nullptr)
269 270
    {
    delete [] this->FileName;
271 272
    this->FileName = nullptr;
    this->Internal->SetFileName( nullptr );
273 274 275 276 277 278
    }
    this->FileName = new char[  strlen( fileName ) + 1  ];
    strcpy( this->FileName, fileName );
    this->FileName[ strlen( fileName ) ] = '\0';
    this->Internal->SetFileName( this->FileName );
    this->ParseConversionFactors( );
279 280 281 282

    this->Internal->ReadMetaData();
    this->SetUpDataArraySelections();
    this->InitializeArraySelections();
283
  }
284

285

286 287 288
  this->Modified();
}

289
//-----------------------------------------------------------------------------
290 291
void vtkAMREnzoReader::ReadMetaData()
{
292
  assert( "pre: Internal Enzo Reader is nullptr" && (this->Internal != nullptr) );
293 294

  if( !this->IsReady )
295
  {
296
    return;
297
  }
298

299
  this->Internal->ReadMetaData();
300 301
}

302

303
//-----------------------------------------------------------------------------
304 305
int vtkAMREnzoReader::GetBlockLevel( const int blockIdx )
{
306
  assert( "pre: Internal Enzo Reader is nullptr" && (this->Internal != nullptr) );
307

308
  if( !this->IsReady )
309
  {
310
    return( -1 );
311
  }
312

313
  this->Internal->ReadMetaData();
314

315
  if( blockIdx < 0 || blockIdx >= this->Internal->NumberOfBlocks )
316
  {
317 318
    vtkErrorMacro( "Block Index (" << blockIdx << ") is out-of-bounds!" );
    return( -1 );
319
  }
320
  return( this->Internal->Blocks[ blockIdx+1 ].Level );
321 322
}

323
//-----------------------------------------------------------------------------
324 325
int vtkAMREnzoReader::GetNumberOfBlocks()
{
326
  assert( "pre: Internal Enzo Reader is nullptr" && (this->Internal != nullptr) );
327
  if( !this->IsReady )
328
  {
329
    return 0;
330
  }
331

332 333
  this->Internal->ReadMetaData();
  return( this->Internal->NumberOfBlocks );
334 335
}

336
//-----------------------------------------------------------------------------
337 338
int vtkAMREnzoReader::GetNumberOfLevels()
{
339
  assert( "pre: Internal Enzo Reader is nullptr" && (this->Internal != nullptr) );
340
  if( !this->IsReady )
341
  {
342
    return 0;
343
  }
344

345 346
  this->Internal->ReadMetaData();
  return( this->Internal->NumberOfLevels );
347 348
}

349
//-----------------------------------------------------------------------------
350
int vtkAMREnzoReader::FillMetaData( )
351
{
352 353
  assert( "pre: Internal Enzo Reader is nullptr" && (this->Internal != nullptr) );
  assert( "pre: metadata object is nullptr" && (this->Metadata != nullptr) );
354
  if( !this->IsReady )
355
  {
356
    return 0;
357
  }
358

359 360
  this->Internal->ReadMetaData();

361 362 363 364
  double origin[3];
  std::vector<int> blocksPerLevel;
  this->ComputeStats(this->Internal, blocksPerLevel, origin);

365 366 367
  this->Metadata->Initialize( static_cast<int>(blocksPerLevel.size()), &blocksPerLevel[0]);
  this->Metadata->SetGridDescription(VTK_XYZ_GRID);
  this->Metadata->SetOrigin(origin);
368 369 370

  std::vector< int > b2level( this->Internal->NumberOfLevels+1, 0 );
  for( int block=0; block < this->Internal->NumberOfBlocks; ++block )
371
  {
372
    vtkEnzoReaderBlock &theBlock = this->Internal->Blocks[ block+1 ];
373
    int level                    = theBlock.Level;
374
    int internalIdx              = block;
375
    int id                       = b2level[ level ];
376 377 378 379

    //compute spacing
    double spacing[3];
    for(int d=0; d<3; ++d)
380
    {
381 382
      spacing[d] = (theBlock.BlockNodeDimensions[d] > 1)?
        (theBlock.MaxBounds[d]-theBlock.MinBounds[d])/(theBlock.BlockNodeDimensions[d]-1.0):1.0;
383
    }
384 385 386 387 388 389 390
    //compute AMRBox
    vtkAMRBox box(theBlock.MinBounds, theBlock.BlockNodeDimensions, spacing, origin, VTK_XYZ_GRID);

    //set meta data
    this->Metadata->SetSpacing(level,spacing);
    this->Metadata->SetAMRBox(level, id, box);
    this->Metadata->SetAMRBlockSourceIndex(level,id, internalIdx);
391
    b2level[ level ]++;
392
  }
393
  this->Metadata->GenerateParentChildInformation();
394
  this->Metadata->GetInformation()->Set(vtkDataObject::DATA_TIME_STEP(),this->Internal->DataTime);
395 396 397
  return( 1 );
}

398
//-----------------------------------------------------------------------------
399
vtkUniformGrid* vtkAMREnzoReader::GetAMRGrid( const int blockIdx )
400
{
401
  assert( "pre: Internal Enzo Reader is nullptr" && (this->Internal != nullptr) );
402

403
  if( !this->IsReady )
404
  {
405
    return nullptr;
406
  }
407

408
  this->Internal->ReadMetaData();
409 410

  // this->Internal->Blocks includes a pseudo block --- the root as block #0
411 412 413 414 415 416
  vtkEnzoReaderBlock &theBlock = this->Internal->Blocks[ blockIdx+1 ];
  double blockMin[3];
  double blockMax[3];
  double spacings[3];

  for( int i=0; i < 3; ++i )
417
  {
418 419 420 421
    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;
422
  }
423 424 425

  vtkUniformGrid *ug = vtkUniformGrid::New();
  ug->SetDimensions( theBlock.BlockNodeDimensions );
426 427 428 429
  ug->SetOrigin( blockMin[0], blockMin[1], blockMin[2] );
  ug->SetSpacing( spacings[0], spacings[1], spacings[2] );
  return( ug );
}
430

431 432 433 434
//-----------------------------------------------------------------------------
void vtkAMREnzoReader::GetAMRGridData(
    const int blockIdx, vtkUniformGrid *block, const char *field)
{
435
  assert( "pre: AMR block is nullptr" && (block != nullptr));
436

437 438
  this->Internal->GetBlockAttribute( field, blockIdx, block );
  if( this->ConvertToCGS == 1 )
439
  {
440 441
    double conversionFactor = this->GetConversionFactor(field);
    if( conversionFactor != 1.0 )
442
    {
443
      vtkDataArray *data = block->GetCellData()->GetArray( field );
444
      assert( "pre: data array is nullptr!" && (data != nullptr) );
445 446 447

      vtkIdType numTuples = data->GetNumberOfTuples();
      for( vtkIdType t=0; t < numTuples; ++t )
448
      {
449 450
        int numComp = data->GetNumberOfComponents();
        for( int c=0; c < numComp; ++c )
451
        {
452 453
         double f = data->GetComponent( t, c );
         data->SetComponent( t, c, f*conversionFactor );
454 455 456 457
        } // END for all components
      } // END for all tuples
    } // END if the conversion factor is not 1.0
  } // END if conversion to CGS units is requested
458 459
}

460

461 462
//-----------------------------------------------------------------------------
void vtkAMREnzoReader::SetUpDataArraySelections()
463
{
464
  assert( "pre: Internal Enzo Reader is nullptr" && (this->Internal != nullptr) );
465 466
  this->Internal->ReadMetaData();
  this->Internal->GetAttributeNames();
467

468
  int numAttrs = static_cast< int >(this->Internal->BlockAttributeNames.size());
469
  for( int i=0; i < numAttrs; i++ )
470
  {
471 472
    this->CellDataArraySelection->AddArray(
        this->Internal->BlockAttributeNames[i].c_str() );
473
  } // END for all attributes
474

475
}