vtkAMREnzoParticlesReader.cxx 13.4 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
/*=========================================================================

 Program:   Visualization Toolkit
 Module:    vtkAMREnzoParticlesReader.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.

 =========================================================================*/
#include "vtkAMREnzoParticlesReader.h"
#include "vtkObjectFactory.h"
#include "vtkPolyData.h"
18
#include "vtkCellArray.h"
19 20
#include "vtkDataArraySelection.h"
#include "vtkPointData.h"
21 22
#include "vtkIntArray.h"
#include "vtkDataArray.h"
23 24

#include <cassert>
25
#include <vector>
26

27 28
#include "vtksys/SystemTools.hxx"

29
#define H5_USE_16_API
30
#include "vtk_hdf5.h"      // for the HDF data loading engine
31

32 33 34 35 36 37 38 39 40 41
#include "vtkAMREnzoReaderInternal.h"

//------------------------------------------------------------------------------
//            HDF5 Utility Routines
//------------------------------------------------------------------------------

//------------------------------------------------------------------------------
// Description:
// Finds the block index (blockIndx) within the HDF5 file associated with
// the given file index.
42
static bool FindBlockIndex( hid_t fileIndx, const int blockIdx, hid_t &rootIndx )
43 44 45 46 47 48
{
  // retrieve the contents of the root directory to look for a group
  // corresponding to the target block, if available, open that group
  hsize_t numbObjs;
  rootIndx = H5Gopen( fileIndx, "/" );
  if( rootIndx < 0 )
49
  {
50 51
    vtkGenericWarningMacro( "Failed to open root node of particles file" );
    return false;
52
  }
53 54 55 56

  bool found = false;
  H5Gget_num_objs( rootIndx, &numbObjs );
  for( int objIndex=0; objIndex < static_cast < int >(numbObjs); objIndex++  )
57
  {
58
    if( H5Gget_objtype_by_idx( rootIndx, objIndex ) == H5G_GROUP  )
59
    {
60 61 62 63 64 65 66
      int    blckIndx;
      char   blckName[65];
      H5Gget_objname_by_idx( rootIndx, objIndex, blckName, 64 );

      // Is this the target block?
      if( (sscanf( blckName, "Grid%d", &blckIndx )==1) &&
          (blckIndx  ==  blockIdx) )
67
      {
68 69 70
        // located the target block
        rootIndx = H5Gopen( rootIndx, blckName );
        if( rootIndx < 0 )
71
        {
72
          vtkGenericWarningMacro( "Could not locate target block!\n" );
73
        }
74 75
        found = true;
        break;
76 77 78
      }
    } // END if group
  } // END for all objects
79 80 81 82 83 84
  return( found );
}

//------------------------------------------------------------------------------
// Description:
// Returns the double array
85
static void GetDoubleArrayByName(
86
    const hid_t rootIdx, const char* name, std::vector<double> &array)
87 88
{
  // turn off warnings
89
  void       * pContext = nullptr;
90 91
  H5E_auto_t   erorFunc;
  H5Eget_auto( &erorFunc, &pContext );
92
  H5Eset_auto( nullptr, nullptr );
93 94 95

  hid_t arrayIdx = H5Dopen( rootIdx, name );
  if( arrayIdx < 0 )
96
  {
97 98
      vtkGenericWarningMacro( "Cannot open array: " << name << "\n");
      return;
99
  }
100 101 102

  // turn warnings back on
  H5Eset_auto( erorFunc, pContext );
103
  pContext = nullptr;
104 105 106 107

  // get the number of particles
  hsize_t      dimValus[3];
  hid_t        spaceIdx = H5Dget_space( arrayIdx );
108
  H5Sget_simple_extent_dims( spaceIdx, dimValus, nullptr );
109 110 111 112 113
  int          numbPnts = dimValus[0];

  array.resize( numbPnts );
  H5Dread( arrayIdx,H5T_NATIVE_DOUBLE,H5S_ALL,H5S_ALL,H5P_DEFAULT,&array[0] );

114 115
//  H5Dclose( spaceIdx );
//  H5Dclose( arrayIdx );
116 117 118 119 120
}

//------------------------------------------------------------------------------
//          END of HDF5 Utility Routine definitions
//------------------------------------------------------------------------------
121 122 123 124 125 126

vtkStandardNewMacro(vtkAMREnzoParticlesReader);

//------------------------------------------------------------------------------
vtkAMREnzoParticlesReader::vtkAMREnzoParticlesReader()
{
127 128
  this->Internal     = new vtkEnzoReaderInternal();
  this->ParticleType = -1; /* undefined particle type */
129 130 131 132 133 134
  this->Initialize();
}

//------------------------------------------------------------------------------
vtkAMREnzoParticlesReader::~vtkAMREnzoParticlesReader()
{
135
  delete this->Internal;
136
  this->Internal = nullptr;
137 138 139 140 141 142 143 144 145 146 147 148
}

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

//------------------------------------------------------------------------------
void vtkAMREnzoParticlesReader::ReadMetaData()
{
149
  if( this->Initialized )
150
  {
151
    return;
152
  }
153

154
  if( !this->FileName )
155
  {
156 157
    vtkErrorMacro("No FileName set!");
    return;
158
  }
159

160
  this->Internal->SetFileName( this->FileName );
161 162 163
  std::string  tempName( this->FileName );
  std::string  bExtName( ".boundary" );
  std::string  hExtName( ".hierarchy" );
164 165 166

  if( tempName.length() > hExtName.length() &&
       tempName.substr(tempName.length()-hExtName.length() )== hExtName )
167
  {
168 169 170 171 172
     this->Internal->MajorFileName =
         tempName.substr( 0, tempName.length() - hExtName.length() );
     this->Internal->HierarchyFileName = tempName;
     this->Internal->BoundaryFileName  =
         this->Internal->MajorFileName + bExtName;
173
  }
174 175
  else if( tempName.length() > bExtName.length() &&
      tempName.substr( tempName.length() - bExtName.length() )==bExtName )
176
  {
177 178 179 180 181
    this->Internal->MajorFileName =
       tempName.substr( 0, tempName.length() - bExtName.length() );
    this->Internal->BoundaryFileName  = tempName;
    this->Internal->HierarchyFileName =
       this->Internal->MajorFileName + hExtName;
182
  }
183
  else
184
  {
185 186
   vtkErrorMacro( "Enzo file has invalid extension!");
   return;
187
  }
188 189 190 191

   this->Internal->DirectoryName =
       GetEnzoDirectory(this->Internal->MajorFileName.c_str());

192
  this->Internal->ReadMetaData();
193
  this->Internal->CheckAttributeNames();
194

195 196
  this->NumberOfBlocks = this->Internal->NumberOfBlocks;
  this->Initialized    = true;
197 198

  this->SetupParticleDataSelections();
199 200
}

201 202 203 204 205 206 207
//------------------------------------------------------------------------------
vtkDataArray* vtkAMREnzoParticlesReader::GetParticlesTypeArray(
    const int blockIdx )
{

  vtkIntArray *array = vtkIntArray::New();
  if( this->ParticleDataArraySelection->ArrayExists( "particle_type" ) )
208
  {
209 210
    this->Internal->LoadAttribute( "particle_type", blockIdx );
    array->DeepCopy( this->Internal->DataArray );
211
  }
212 213 214 215 216 217 218
  return( array );
}

//------------------------------------------------------------------------------
bool vtkAMREnzoParticlesReader::CheckParticleType(
    const int idx, vtkIntArray *ptypes )
{
219
  assert( "pre: particles type array should not be nullptr" && (ptypes != nullptr) );
220 221 222

  if( ptypes->GetNumberOfTuples() > 0 &&
      this->ParticleDataArraySelection->ArrayExists( "particle_type" ) )
223
  {
224 225
    int ptype = ptypes->GetValue( idx );
    if( (this->ParticleType==0) || (ptype==this->ParticleType) )
226
    {
227
      return true;
228
    }
229
    else
230
    {
231
      return false;
232
    }
233
  }
234
  else
235
  {
236
    return true;
237
  }
238 239
}

240
//------------------------------------------------------------------------------
241 242
vtkPolyData* vtkAMREnzoParticlesReader::GetParticles(
    const char* file, const int blockIdx )
243 244
{
  vtkPolyData *particles = vtkPolyData::New();
245 246 247 248 249 250
  vtkPoints   *positions = vtkPoints::New();
  positions->SetDataTypeToDouble();
  vtkPointData *pdata    = particles->GetPointData();

  hid_t fileIndx = H5Fopen( file, H5F_ACC_RDONLY, H5P_DEFAULT );
  if( fileIndx < 0 )
251
  {
252
    vtkErrorMacro( "Failed opening particles file!" );
253
    return nullptr;
254
  }
255 256

  hid_t rootIndx;
257
  if( ! FindBlockIndex( fileIndx, blockIdx+1,rootIndx ) )
258
  {
259
    vtkErrorMacro( "Could not locate target block!" );
260
    return nullptr;
261
  }
262 263 264 265 266 267 268 269

  //
  // Load the particles position arrays by name.
  // In Enzo the following arrays are available:
  //  ( 1 ) particle_position_i
  //  ( 2 ) tracer_particle_position_i
  //
  // where i \in {x,y,z}.
270 271 272
  std::vector< double > xcoords;
  std::vector< double > ycoords;
  std::vector< double > zcoords;
273

274
  // TODO: should we handle 2-D particle datasets?
275 276 277
  GetDoubleArrayByName( rootIndx, "particle_position_x", xcoords );
  GetDoubleArrayByName( rootIndx, "particle_position_y", ycoords );
  GetDoubleArrayByName( rootIndx, "particle_position_z", zcoords );
278

279
  vtkIntArray *particleTypes = vtkArrayDownCast<vtkIntArray>(
280 281
      this->GetParticlesTypeArray( blockIdx ) );

282 283 284 285 286
  assert( "Coordinate arrays must have the same size: " &&
           (xcoords.size()==ycoords.size() ) );
  assert( "Coordinate arrays must have the same size: " &&
           (ycoords.size()==zcoords.size() ) );

287
  int TotalNumberOfParticles = static_cast<int>(xcoords.size());
288
  positions->SetNumberOfPoints( TotalNumberOfParticles );
289 290 291 292

  vtkIdList *ids = vtkIdList::New();
  ids->SetNumberOfIds( TotalNumberOfParticles );

293
  vtkIdType NumberOfParticlesLoaded = 0;
294
  for( int i=0; i < TotalNumberOfParticles; ++i )
295
  {
296
    if( (i%this->Frequency) == 0  )
297
    {
298 299
      if( this->CheckLocation( xcoords[i], ycoords[i],zcoords[i] ) &&
          this->CheckParticleType( i, particleTypes) )
300
      {
301 302 303 304
        int pidx = NumberOfParticlesLoaded;
        ids->InsertId( pidx, i );
        positions->SetPoint( pidx, xcoords[i],ycoords[i],zcoords[i] );
        ++NumberOfParticlesLoaded;
305 306 307
      } // END if within requested region
    } // END if within requested interval
  } // END for all particles
308 309 310
  H5Gclose( rootIndx );
  H5Fclose( fileIndx );

311 312 313
  ids->SetNumberOfIds( NumberOfParticlesLoaded );
  ids->Squeeze();

314 315
  positions->SetNumberOfPoints( NumberOfParticlesLoaded );
  positions->Squeeze();
316

317 318
  particles->SetPoints( positions );
  positions->Delete();
319 320 321 322 323 324 325 326 327 328 329

  // Create CellArray consisting of a single polyvertex cell
  vtkCellArray *polyVertex     = vtkCellArray::New();

  polyVertex->InsertNextCell( NumberOfParticlesLoaded  );
  for( vtkIdType idx=0; idx < NumberOfParticlesLoaded; ++idx )
    polyVertex->InsertCellPoint( idx );

  particles->SetVerts( polyVertex );
  polyVertex->Delete();

330 331
  // Release the particle types array
  particleTypes->Delete();
332 333 334

  int numArrays = this->ParticleDataArraySelection->GetNumberOfArrays();
  for( int i=0; i < numArrays; ++i )
335
  {
336 337
    const char* name = this->ParticleDataArraySelection->GetArrayName( i );
    if( this->ParticleDataArraySelection->ArrayIsEnabled( name ) )
338
    {
339 340 341 342 343 344 345 346 347 348 349 350 351 352
      // Note: 0-based indexing is used for loading particles
      this->Internal->LoadAttribute( name, blockIdx );
      assert( "pre: particle attribute size mismatch" &&
         (this->Internal->DataArray->GetNumberOfTuples()==
          TotalNumberOfParticles) );

      vtkDataArray *array = this->Internal->DataArray->NewInstance();
      array->SetName(this->Internal->DataArray->GetName( )  );
      array->SetNumberOfTuples( NumberOfParticlesLoaded );
      array->SetNumberOfComponents(
          this->Internal->DataArray->GetNumberOfComponents() );

      vtkIdType numIds = ids->GetNumberOfIds();
      for( vtkIdType pidx=0; pidx < numIds; ++pidx )
353
      {
354 355 356
          vtkIdType particleIdx = ids->GetId( pidx );
          int numComponents = array->GetNumberOfComponents();
          for( int k=0; k < numComponents; ++k )
357
          {
358 359 360
              array->SetComponent( pidx, k,
                  this->Internal->DataArray->GetComponent(
                      particleIdx,k ) );
361 362
          } // END for all components
      } // END for all ids of loaded particles
363 364
      pdata->AddArray( array );
      array->Delete();
365 366
    } // END if the array is supposed to be loaded
  } // END for all particle arrays
367 368

  ids->Delete();
369 370 371
  return( particles );
}

372 373 374
//------------------------------------------------------------------------------
void vtkAMREnzoParticlesReader::SetupParticleDataSelections()
{
375
  assert( "pre: Intenal reader is nullptr" && (this->Internal != nullptr) );
376

377 378
  unsigned int N =
    static_cast<unsigned int>( this->Internal->ParticleAttributeNames.size() );
379
  for( unsigned int i=0; i < N; ++i )
380
  {
381 382 383 384 385 386
    bool isParticleAttribute =
     vtksys::SystemTools::StringStartsWith(
       this->Internal->ParticleAttributeNames[ i ].c_str( ),
        "particle_" );

    if( isParticleAttribute )
387
    {
388 389
      this->ParticleDataArraySelection->AddArray(
          this->Internal->ParticleAttributeNames[ i ].c_str() );
390 391
    }
  } // END for all particles attributes
392 393 394
  this->InitializeParticleDataSelections();
}

395 396 397
//------------------------------------------------------------------------------
int vtkAMREnzoParticlesReader::GetTotalNumberOfParticles( )
{
398
  assert( "Internal reader is null" && (this->Internal!=nullptr) );
399 400
  int numParticles = 0;
  for( int blockIdx=0; blockIdx < this->NumberOfBlocks; ++blockIdx )
401
  {
402
    numParticles += this->Internal->Blocks[ blockIdx ].NumberOfParticles;
403
  }
404 405 406
  return( numParticles );
}

407 408 409 410 411 412 413 414
//------------------------------------------------------------------------------
vtkPolyData* vtkAMREnzoParticlesReader::ReadParticles(const int blkidx)
{
  // this->Internal->Blocks includes a pseudo block -- the roo as block #0
  int iBlockIdx    = blkidx+1;
  int NumParticles = this->Internal->Blocks[ iBlockIdx ].NumberOfParticles;

  if( NumParticles <= 0 )
415
  {
416
    vtkPolyData* emptyParticles = vtkPolyData::New();
417
    assert( "Cannot create particles dataset" && ( emptyParticles != nullptr ) );
418
    return( emptyParticles );
419
  }
420

421
  std::string pfile = this->Internal->Blocks[iBlockIdx].ParticleFileName;
422
  if( pfile.empty() )
423
  {
424
    vtkErrorMacro( "No particles file found, string is empty!" );
425
    return nullptr;
426
  }
427 428

  vtkPolyData* particles = this->GetParticles(pfile.c_str(), blkidx  );
429 430
  return( particles );
}