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

 Program:   Visualization Toolkit
 Module:    vtkAMRFlashParticlesReader.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 "vtkAMRFlashParticlesReader.h"
#include "vtkObjectFactory.h"
#include "vtkPolyData.h"
18
#include "vtkDataArraySelection.h"
19 20 21 22 23
#include "vtkIdList.h"
#include "vtkCellArray.h"
#include "vtkDoubleArray.h"
#include "vtkIntArray.h"
#include "vtkPointData.h"
24

25 26
#include "vtkAMRFlashReaderInternal.h"

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

30
#include <vector>
31 32
#include <cassert>

33 34 35 36 37 38 39 40 41 42
#define  FLASH_READER_MAX_DIMS     3
#define  FLASH_READER_LEAF_BLOCK   1
#define  FLASH_READER_FLASH3_FFV8  8
#define  FLASH_READER_FLASH3_FFV9  9

//------------------------------------------------------------------------------
// Description:
// Helper function that reads the particle coordinates
// NOTE: it is assumed that H5DOpen has been called on the
// internal file index this->FileIndex.
43
static void GetParticleCoordinates( hid_t &dataIdx,
44 45
  std::vector< double > &xcoords, std::vector< double > &ycoords,
  std::vector< double > &zcoords,
46 47 48 49
  vtkFlashReaderInternal *iReader,
  int NumParticles )
{

50
  assert( "pre: internal reader should not be nullptr" && (iReader != nullptr) );
51 52

  hid_t theTypes[3];
Dave DeMarle's avatar
Dave DeMarle committed
53
  theTypes[0] = theTypes[1] = theTypes[2] = H5I_UNINIT;
54 55 56 57 58
  xcoords.resize( NumParticles );
  ycoords.resize( NumParticles );
  zcoords.resize( NumParticles );

  if( iReader->FileFormatVersion < FLASH_READER_FLASH3_FFV8  )
59
  {
60 61 62 63 64 65
    theTypes[0] = H5Tcreate( H5T_COMPOUND, sizeof(double) );
    theTypes[1] = H5Tcreate( H5T_COMPOUND, sizeof(double) );
    theTypes[2] = H5Tcreate( H5T_COMPOUND, sizeof(double) );
    H5Tinsert( theTypes[0], "particle_x", 0, H5T_NATIVE_DOUBLE );
    H5Tinsert( theTypes[1], "particle_y", 0, H5T_NATIVE_DOUBLE );
    H5Tinsert( theTypes[2], "particle_z", 0, H5T_NATIVE_DOUBLE );
66
  }
67 68 69

  // Read the coordinates from the file
  switch( iReader->NumberOfDimensions )
70
  {
71 72
    case 1:
      if( iReader->FileFormatVersion < FLASH_READER_FLASH3_FFV8 )
73
      {
74
        H5Dread(dataIdx,theTypes[0],H5S_ALL,H5S_ALL,H5P_DEFAULT,&xcoords[0]);
75
      }
76
      else
77
      {
78
        iReader->ReadParticlesComponent(dataIdx,"Particles/posx",&xcoords[0]);
79
      }
80 81 82
      break;
    case 2:
      if( iReader->FileFormatVersion < FLASH_READER_FLASH3_FFV8 )
83
      {
84 85
        H5Dread(dataIdx,theTypes[0],H5S_ALL,H5S_ALL,H5P_DEFAULT,&xcoords[0]);
        H5Dread(dataIdx,theTypes[1],H5S_ALL,H5S_ALL,H5P_DEFAULT,&ycoords[0]);
86
      }
87
      else
88
      {
89 90
        iReader->ReadParticlesComponent(dataIdx,"Particles/posx",&xcoords[0]);
        iReader->ReadParticlesComponent(dataIdx,"Particles/posy",&ycoords[0]);
91
      }
92 93 94
      break;
    case 3:
      if( iReader->FileFormatVersion < FLASH_READER_FLASH3_FFV8 )
95
      {
96 97 98
        H5Dread(dataIdx,theTypes[0],H5S_ALL,H5S_ALL,H5P_DEFAULT,&xcoords[0]);
        H5Dread(dataIdx,theTypes[1],H5S_ALL,H5S_ALL,H5P_DEFAULT,&ycoords[0]);
        H5Dread(dataIdx,theTypes[2],H5S_ALL,H5S_ALL,H5P_DEFAULT,&zcoords[0]);
99
      }
100
      else
101
      {
102 103 104
        iReader->ReadParticlesComponent(dataIdx,"Particles/posx",&xcoords[0]);
        iReader->ReadParticlesComponent(dataIdx,"Particles/posy",&ycoords[0]);
        iReader->ReadParticlesComponent(dataIdx,"Particles/posz",&zcoords[0]);
105
      }
106 107 108 109 110
      break;
    default:
      std::cerr << "ERROR: Undefined dimension!\n" << std::endl;
      std::cerr.flush();
      return;
111
  }
112 113 114 115 116 117 118 119 120

}


//------------------------------------------------------------------------------


//------------------------------------------------------------------------------

121 122 123 124
vtkStandardNewMacro( vtkAMRFlashParticlesReader );

vtkAMRFlashParticlesReader::vtkAMRFlashParticlesReader()
{
125
  this->Internal = new vtkFlashReaderInternal();
126
  this->Initialized = false;
127 128 129 130 131 132
  this->Initialize();
}

//------------------------------------------------------------------------------
vtkAMRFlashParticlesReader::~vtkAMRFlashParticlesReader()
{
133
  delete this->Internal;
134 135 136 137 138 139 140 141 142 143 144
}

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

//------------------------------------------------------------------------------
void vtkAMRFlashParticlesReader::ReadMetaData()
{
145
  if( this->Initialized )
146
  {
147
    return;
148
  }
149 150 151

  this->Internal->SetFileName( this->FileName );
  this->Internal->ReadMetaData();
152 153 154 155 156

  // In some cases, the FLASH file format may have no blocks but store
  // just particles in a single block. However, the AMRBaseParticles reader
  // would expect that in this case, the number of blocks is set to 1. The
  // following lines of code provide a simple workaround for that.
157
  this->NumberOfBlocks = this->Internal->NumberOfBlocks;
158
  if( this->NumberOfBlocks == 0 && this->Internal->NumberOfParticles > 0)
159
  {
160
    this->NumberOfBlocks=1;
161
  }
162 163
  this->Initialized    = true;
  this->SetupParticleDataSelections();
164 165
}

166 167 168
//------------------------------------------------------------------------------
int vtkAMRFlashParticlesReader::GetTotalNumberOfParticles()
{
169
  assert( "Internal reader is null" && (this->Internal!=nullptr) );
170 171 172
  return( this->Internal->NumberOfParticles );
}

173 174
//------------------------------------------------------------------------------
vtkPolyData* vtkAMRFlashParticlesReader::GetParticles(
George Zagaris's avatar
George Zagaris committed
175
    const char *file, const int vtkNotUsed(blkidx) )
176 177 178
{
  hid_t dataIdx = H5Dopen( this->Internal->FileIndex, file );
  if( dataIdx < 0 )
179
  {
180
    vtkErrorMacro( "Could not open particles file!" );
181
    return nullptr;
182
  }
183 184 185 186 187 188 189

  vtkPolyData *particles = vtkPolyData::New();
  vtkPoints   *positions = vtkPoints::New();
  positions->SetDataTypeToDouble();
  positions->SetNumberOfPoints( this->Internal->NumberOfParticles );

  vtkPointData *pdata = particles->GetPointData();
190
  assert( "pre: PointData is nullptr" && (pdata != nullptr) );
191 192

  // Load the particle position arrays by name
193 194 195
  std::vector< double > xcoords;
  std::vector< double > ycoords;
  std::vector< double > zcoords;
196 197 198 199 200
  GetParticleCoordinates(
      dataIdx, xcoords, ycoords, zcoords,
      this->Internal, this->Internal->NumberOfParticles );

  // Sub-sample particles
201
  int TotalNumberOfParticles        = static_cast<int>(xcoords.size());
202 203 204 205 206
  vtkIdList *ids                    = vtkIdList::New();
  ids->SetNumberOfIds( TotalNumberOfParticles );

  vtkIdType NumberOfParticlesLoaded = 0;
  for( int i=0; i < TotalNumberOfParticles; ++i )
207
  {
208
    if( i%this->Frequency == 0 )
209
    {
210
      if( this->CheckLocation(xcoords[i],ycoords[i],zcoords[i] ) )
211
      {
212 213 214 215
        int pidx = NumberOfParticlesLoaded;
        ids->InsertId( pidx, i );
        positions->SetPoint( pidx, xcoords[i], ycoords[i], zcoords[i] );
        ++NumberOfParticlesLoaded;
216 217 218
      } // END if within requested region
    } // END if within requested interval
  } // END for all particles
219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234

  xcoords.clear(); ycoords.clear(); zcoords.clear();

  ids->SetNumberOfIds( NumberOfParticlesLoaded );
  ids->Squeeze();

  positions->SetNumberOfPoints( NumberOfParticlesLoaded );
  positions->Squeeze();

  particles->SetPoints( positions );
  positions->Squeeze( );

  // Create CellArray consisting of a single polyvertex
  vtkCellArray *polyVertex = vtkCellArray::New();
  polyVertex ->InsertNextCell( NumberOfParticlesLoaded );
  for( vtkIdType idx=0; idx < NumberOfParticlesLoaded; ++idx )
235
  {
236
    polyVertex->InsertCellPoint( idx );
237
  }
238 239 240 241 242 243
  particles->SetVerts( polyVertex );
  polyVertex->Delete();

  // Load particle data arrays
  int numArrays = this->ParticleDataArraySelection->GetNumberOfArrays();
  for( int i=0; i < numArrays; ++i )
244
  {
245 246
    const char *name = this->ParticleDataArraySelection->GetArrayName( i );
    if( this->ParticleDataArraySelection->ArrayIsEnabled( name ) )
247
    {
248 249 250 251
      int attrIdx     = this->Internal->ParticleAttributeNamesToIds[ name ];
      hid_t attrType  = this->Internal->ParticleAttributeTypes[ attrIdx ];

      if( attrType == H5T_NATIVE_DOUBLE )
252
      {
253
        double *data = new double[ this->Internal->NumberOfParticles ];
254
        assert( data != nullptr );
255 256

        if( this->Internal->FileFormatVersion < FLASH_READER_FLASH3_FFV8 )
257
        {
258 259 260 261
          hid_t dataType = H5Tcreate( H5T_COMPOUND, sizeof(double) );
          H5Tinsert( dataType, name, 0, H5T_NATIVE_DOUBLE );
          H5Dread(dataIdx, dataType, H5S_ALL, H5S_ALL, H5P_DEFAULT,data);
          H5Tclose( dataType );
262
        }
263
        else
264
        {
265
          this->Internal->ReadParticlesComponent( dataIdx, name, data );
266
        }
267 268 269 270 271 272 273 274

        vtkDataArray *array = vtkDoubleArray::New();
        array->SetName( name );
        array->SetNumberOfTuples( ids->GetNumberOfIds() );
        array->SetNumberOfComponents( 1 );

        vtkIdType numIds = ids->GetNumberOfIds();
        for( vtkIdType pidx=0; pidx < numIds; ++pidx )
275
        {
276 277
          vtkIdType particleIdx = ids->GetId( pidx );
          array->SetComponent( pidx, 0, data[ particleIdx ] );
278
        } // END for all ids of loaded particles
279 280
        pdata->AddArray( array );
        delete [] data;
281
      }
282
      else if( attrType == H5T_NATIVE_INT )
283
      {
284 285 286 287
        hid_t dataType = H5Tcreate( H5T_COMPOUND, sizeof(int) );
        H5Tinsert( dataType, name, 0, H5T_NATIVE_INT );

        int *data = new int[ this->Internal->NumberOfParticles ];
288
        assert( data != nullptr );
289
        H5Dread( dataIdx, dataType, H5S_ALL, H5S_ALL, H5P_DEFAULT, data );
290

291 292 293 294
        vtkDataArray *array = vtkIntArray::New();
        array->SetName( name );
        array->SetNumberOfTuples( ids->GetNumberOfIds() );
        array->SetNumberOfComponents( 1 );
295

296 297
        vtkIdType numIds = ids->GetNumberOfIds( );
        for( vtkIdType pidx=0; pidx < numIds; ++pidx )
298
        {
299 300
          vtkIdType particleIdx = ids->GetId( pidx );
          array->SetComponent( pidx, 0, data[ particleIdx ] );
301
        } // END for all ids of loaded particles
302 303
        pdata->AddArray( array );
        delete [] data;
304
      }
305
      else
306
      {
luz.paz's avatar
luz.paz committed
307
        vtkErrorMacro( "Unsupported array type in HDF5 file!" );
308
        return nullptr;
309 310 311
      }
    } // END if the array is supposed to be loaded
  } // END for all arrays
312 313 314 315 316 317

  H5Dclose(dataIdx);

  return( particles );
}

318 319 320
//------------------------------------------------------------------------------
vtkPolyData* vtkAMRFlashParticlesReader::ReadParticles( const int blkidx )
{
321
  assert( "pre: Internal reader is nullptr" && (this->Internal != nullptr) );
322 323 324 325
  assert( "pre: Not initialized " && (this->Initialized) );

  int NumberOfParticles = this->Internal->NumberOfParticles;
  if( NumberOfParticles <= 0 )
326
  {
327
    vtkPolyData *emptyParticles = vtkPolyData::New();
328
    assert( "Cannot create particle dataset" && (emptyParticles != nullptr)  );
329
    return( emptyParticles );
330
  }
331 332 333

  vtkPolyData* particles = this->GetParticles(
     this->Internal->ParticleName.c_str(), blkidx );
334
  assert( "partciles should not be nullptr " && (particles != nullptr) );
335

336 337
  return( particles );
}
338 339 340 341

//------------------------------------------------------------------------------
void vtkAMRFlashParticlesReader::SetupParticleDataSelections()
{
342
  assert( "pre: Internal reader is nullptr" && (this->Internal != nullptr) );
343

344 345
  unsigned int N =
      static_cast<unsigned int>(this->Internal->ParticleAttributeNames.size());
346
  for( unsigned int i=0; i < N; ++i )
347
  {
348 349
    this->ParticleDataArraySelection->AddArray(
      this->Internal->ParticleAttributeNames[ i ].c_str( ) );
350
  } // END for all particles attributes
351 352

  this->InitializeParticleDataSelections();
353
}