File: vtkAMRFlashParticlesReader.cxx

package info (click to toggle)
vtk6 6.3.0%2Bdfsg2-8.1
  • links: PTS, VCS
  • area: main
  • in suites: bullseye
  • size: 118,972 kB
  • sloc: cpp: 1,442,790; ansic: 113,395; python: 72,383; tcl: 46,998; xml: 8,119; yacc: 4,525; java: 4,239; perl: 3,108; lex: 1,694; sh: 1,093; asm: 154; makefile: 68; objc: 17
file content (353 lines) | stat: -rw-r--r-- 12,047 bytes parent folder | download | duplicates (5)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
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
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
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
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
/*=========================================================================

 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"
#include "vtkDataArraySelection.h"
#include "vtkIdList.h"
#include "vtkCellArray.h"
#include "vtkDoubleArray.h"
#include "vtkIntArray.h"
#include "vtkPointData.h"

#include "vtkAMRFlashReaderInternal.h"

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

#include <vector>
#include <cassert>

#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.
static void GetParticleCoordinates( hid_t &dataIdx,
  std::vector< double > &xcoords, std::vector< double > &ycoords,
  std::vector< double > &zcoords,
  vtkFlashReaderInternal *iReader,
  int NumParticles )
{

  assert( "pre: internal reader should not be NULL" && (iReader != NULL) );

  hid_t theTypes[3];
  theTypes[0] = theTypes[1] = theTypes[2] = H5I_UNINIT;
  xcoords.resize( NumParticles );
  ycoords.resize( NumParticles );
  zcoords.resize( NumParticles );

  if( iReader->FileFormatVersion < FLASH_READER_FLASH3_FFV8  )
    {
    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 );
    }

  // Read the coordinates from the file
  switch( iReader->NumberOfDimensions )
    {
    case 1:
      if( iReader->FileFormatVersion < FLASH_READER_FLASH3_FFV8 )
        {
        H5Dread(dataIdx,theTypes[0],H5S_ALL,H5S_ALL,H5P_DEFAULT,&xcoords[0]);
        }
      else
        {
        iReader->ReadParticlesComponent(dataIdx,"Particles/posx",&xcoords[0]);
        }
      break;
    case 2:
      if( iReader->FileFormatVersion < FLASH_READER_FLASH3_FFV8 )
        {
        H5Dread(dataIdx,theTypes[0],H5S_ALL,H5S_ALL,H5P_DEFAULT,&xcoords[0]);
        H5Dread(dataIdx,theTypes[1],H5S_ALL,H5S_ALL,H5P_DEFAULT,&ycoords[0]);
        }
      else
        {
        iReader->ReadParticlesComponent(dataIdx,"Particles/posx",&xcoords[0]);
        iReader->ReadParticlesComponent(dataIdx,"Particles/posy",&ycoords[0]);
        }
      break;
    case 3:
      if( iReader->FileFormatVersion < FLASH_READER_FLASH3_FFV8 )
        {
        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]);
        }
      else
        {
        iReader->ReadParticlesComponent(dataIdx,"Particles/posx",&xcoords[0]);
        iReader->ReadParticlesComponent(dataIdx,"Particles/posy",&ycoords[0]);
        iReader->ReadParticlesComponent(dataIdx,"Particles/posz",&zcoords[0]);
        }
      break;
    default:
      std::cerr << "ERROR: Undefined dimension!\n" << std::endl;
      std::cerr.flush();
      return;
    }

}


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


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

vtkStandardNewMacro( vtkAMRFlashParticlesReader );

vtkAMRFlashParticlesReader::vtkAMRFlashParticlesReader()
{
  this->Internal = new vtkFlashReaderInternal();
  this->Initialized = false;
  this->Initialize();
}

//------------------------------------------------------------------------------
vtkAMRFlashParticlesReader::~vtkAMRFlashParticlesReader()
{
  delete this->Internal;
}

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

//------------------------------------------------------------------------------
void vtkAMRFlashParticlesReader::ReadMetaData()
{
  if( this->Initialized )
    {
    return;
    }

  this->Internal->SetFileName( this->FileName );
  this->Internal->ReadMetaData();

  // 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.
  this->NumberOfBlocks = this->Internal->NumberOfBlocks;
  if( this->NumberOfBlocks == 0 && this->Internal->NumberOfParticles > 0)
    {
    this->NumberOfBlocks=1;
    }
  this->Initialized    = true;
  this->SetupParticleDataSelections();
}

//------------------------------------------------------------------------------
int vtkAMRFlashParticlesReader::GetTotalNumberOfParticles()
{
  assert( "Internal reader is null" && (this->Internal!=NULL) );
  return( this->Internal->NumberOfParticles );
}

//------------------------------------------------------------------------------
vtkPolyData* vtkAMRFlashParticlesReader::GetParticles(
    const char *file, const int vtkNotUsed(blkidx) )
{
  hid_t dataIdx = H5Dopen( this->Internal->FileIndex, file );
  if( dataIdx < 0 )
    {
    vtkErrorMacro( "Could not open particles file!" );
    return NULL;
    }

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

  vtkPointData *pdata = particles->GetPointData();
  assert( "pre: PointData is NULL" && (pdata != NULL) );

  // Load the particle position arrays by name
  std::vector< double > xcoords;
  std::vector< double > ycoords;
  std::vector< double > zcoords;
  GetParticleCoordinates(
      dataIdx, xcoords, ycoords, zcoords,
      this->Internal, this->Internal->NumberOfParticles );

  // Sub-sample particles
  int TotalNumberOfParticles        = static_cast<int>(xcoords.size());
  vtkIdList *ids                    = vtkIdList::New();
  ids->SetNumberOfIds( TotalNumberOfParticles );

  vtkIdType NumberOfParticlesLoaded = 0;
  for( int i=0; i < TotalNumberOfParticles; ++i )
    {
    if( i%this->Frequency == 0 )
      {
      if( this->CheckLocation(xcoords[i],ycoords[i],zcoords[i] ) )
        {
        int pidx = NumberOfParticlesLoaded;
        ids->InsertId( pidx, i );
        positions->SetPoint( pidx, xcoords[i], ycoords[i], zcoords[i] );
        ++NumberOfParticlesLoaded;
        } // END if within requested region
      } // END if within requested interval
    } // END for all particles

  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 )
    {
    polyVertex->InsertCellPoint( idx );
    }
  particles->SetVerts( polyVertex );
  polyVertex->Delete();

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

      if( attrType == H5T_NATIVE_DOUBLE )
        {
        double *data = new double[ this->Internal->NumberOfParticles ];
        assert( data != NULL );

        if( this->Internal->FileFormatVersion < FLASH_READER_FLASH3_FFV8 )
          {
          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 );
          }
        else
          {
          this->Internal->ReadParticlesComponent( dataIdx, name, data );
          }

        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 )
          {
          vtkIdType particleIdx = ids->GetId( pidx );
          array->SetComponent( pidx, 0, data[ particleIdx ] );
          } // END for all ids of loaded particles
        pdata->AddArray( array );
        delete [] data;
        }
      else if( attrType == H5T_NATIVE_INT )
        {
        hid_t dataType = H5Tcreate( H5T_COMPOUND, sizeof(int) );
        H5Tinsert( dataType, name, 0, H5T_NATIVE_INT );

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

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

        vtkIdType numIds = ids->GetNumberOfIds( );
        for( vtkIdType pidx=0; pidx < numIds; ++pidx )
          {
          vtkIdType particleIdx = ids->GetId( pidx );
          array->SetComponent( pidx, 0, data[ particleIdx ] );
          } // END for all ids of loaded particles
        pdata->AddArray( array );
        delete [] data;
        }
      else
        {
        vtkErrorMacro( "Unsupport array type in HDF5 file!" );
        return NULL;
        }
      } // END if the array is supposed to be loaded
    } // END for all arrays

  H5Dclose(dataIdx);

  return( particles );
}

//------------------------------------------------------------------------------
vtkPolyData* vtkAMRFlashParticlesReader::ReadParticles( const int blkidx )
{
  assert( "pre: Internal reader is NULL" && (this->Internal != NULL) );
  assert( "pre: Not initialized " && (this->Initialized) );

  int NumberOfParticles = this->Internal->NumberOfParticles;
  if( NumberOfParticles <= 0 )
    {
    vtkPolyData *emptyParticles = vtkPolyData::New();
    assert( "Cannot create particle dataset" && (emptyParticles != NULL)  );
    return( emptyParticles );
    }

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

  return( particles );
}

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

  unsigned int N =
      static_cast<unsigned int>(this->Internal->ParticleAttributeNames.size());
  for( unsigned int i=0; i < N; ++i )
    {
    this->ParticleDataArraySelection->AddArray(
      this->Internal->ParticleAttributeNames[ i ].c_str( ) );
    } // END for all particles attributes

  this->InitializeParticleDataSelections();
}