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

 Program:   Visualization Toolkit
 Module:    vtkAMRBaseReader.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 "vtkAMRBaseReader.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkDataObject.h"
#include "vtkMultiProcessController.h"
20
#include "vtkOverlappingAMR.h"
21
22
#include "vtkDataArraySelection.h"
#include "vtkCallbackCommand.h"
23
#include "vtkIndent.h"
24
25
#include "vtkSmartPointer.h"
#include "vtkCompositeDataPipeline.h"
George Zagaris's avatar
George Zagaris committed
26
#include "vtkStreamingDemandDrivenPipeline.h"
27
28
29
30
31
#include "vtkAMRDataSetCache.h"
#include "vtkUniformGrid.h"
#include "vtkDataArray.h"
#include "vtkCellData.h"
#include "vtkPointData.h"
32
33
#include "vtkAMRUtilities.h"

34
35
#include "vtkTimerLog.h"

36

37
38
#include <cassert>

39

40
41
vtkAMRBaseReader::vtkAMRBaseReader()
{
George Zagaris's avatar
George Zagaris committed
42
  this->LoadedMetaData     = false;
43
44
  this->NumBlocksFromCache = 0;
  this->NumBlocksFromFile  = 0;
45
46
47
48
49
50
51
52
53
54
}

//------------------------------------------------------------------------------
vtkAMRBaseReader::~vtkAMRBaseReader()
{
  this->PointDataArraySelection->RemoveObserver( this->SelectionObserver );
  this->CellDataArraySelection->RemoveObserver( this->SelectionObserver );
  this->SelectionObserver->Delete( );
  this->CellDataArraySelection->Delete( );
  this->PointDataArraySelection->Delete( );
55

56
57
58
59
  if( this->Cache != NULL )
    {
    this->Cache->Delete();
    }
60

61
  if( this->Metadata != NULL )
62
    {
63
    this->Metadata->Delete();
64
    }
65
66
67
68
69
70
}

//------------------------------------------------------------------------------
int vtkAMRBaseReader::FillOutputPortInformation(
    int vtkNotUsed(port),vtkInformation *info )
{
71
  info->Set( vtkDataObject::DATA_TYPE_NAME(), "vtkOverlappingAMR" );
72
73
74
75
76
77
78
79
80
81
82
83
  return 1;
}

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

//------------------------------------------------------------------------------
void vtkAMRBaseReader::Initialize()
{
84
85
  vtkTimerLog::MarkStartEvent( "vtkAMRBaseReader::Initialize" );

86
  this->SetNumberOfInputPorts( 0 );
87
88
  this->FileName       = NULL;
  this->MaxLevel       = 0;
89
  this->Metadata       = NULL;
90
  this->Controller     = vtkMultiProcessController::GetGlobalController();
91
  this->InitialRequest = true;
92
  this->Cache          = vtkAMRDataSetCache::New();
93
94
95
96

  this->CellDataArraySelection  = vtkDataArraySelection::New();
  this->PointDataArraySelection = vtkDataArraySelection::New();
  this->SelectionObserver       = vtkCallbackCommand::New();
97
98
  this->SelectionObserver->SetCallback(
      &vtkAMRBaseReader::SelectionModifiedCallback);
99
100
101
102
103
104
  this->SelectionObserver->SetClientData( this );
  this->CellDataArraySelection->AddObserver(
     vtkCommand::ModifiedEvent,this->SelectionObserver );
  this->PointDataArraySelection->AddObserver(
     vtkCommand::ModifiedEvent, this->SelectionObserver );

105
  vtkTimerLog::MarkEndEvent( "vtkAMRBaseReader::Initialize" );
106
107
108
}

//----------------------------------------------------------------------------
109
110
111
112
113
void vtkAMRBaseReader::SelectionModifiedCallback(
    vtkObject*, unsigned long, void* clientdata, void*)
{
  static_cast<vtkAMRBaseReader*>(clientdata)->Modified();
}
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

//------------------------------------------------------------------------------
int vtkAMRBaseReader::GetNumberOfPointArrays()
{
  return( this->PointDataArraySelection->GetNumberOfArrays() );
}

//------------------------------------------------------------------------------
int vtkAMRBaseReader::GetNumberOfCellArrays()
{
  return( this->CellDataArraySelection->GetNumberOfArrays() );
}

//------------------------------------------------------------------------------
const char* vtkAMRBaseReader::GetPointArrayName(int index)
{
  return( this->PointDataArraySelection->GetArrayName( index ) );
}

//------------------------------------------------------------------------------
const char* vtkAMRBaseReader::GetCellArrayName(int index)
{
  return( this->CellDataArraySelection->GetArrayName( index )  );
}

//------------------------------------------------------------------------------
int vtkAMRBaseReader::GetPointArrayStatus(const char* name)
{
  return( this->PointDataArraySelection->ArrayIsEnabled( name ) );
}

//------------------------------------------------------------------------------
int vtkAMRBaseReader::GetCellArrayStatus(const char* name)
{
  return( this->CellDataArraySelection->ArrayIsEnabled( name ) );
}

//------------------------------------------------------------------------------
void vtkAMRBaseReader::SetPointArrayStatus(const char* name, int status)
{

  if( status )
    {
    this->PointDataArraySelection->EnableArray(name);
    }
  else
    {
    this->PointDataArraySelection->DisableArray(name);
    }

}

//------------------------------------------------------------------------------
void vtkAMRBaseReader::SetCellArrayStatus(const char* name, int status)
{
  if( status )
    {
    this->CellDataArraySelection->EnableArray( name );
    }
  else
    {
    this->CellDataArraySelection->DisableArray( name );
    }
}

//------------------------------------------------------------------------------
int vtkAMRBaseReader::GetBlockProcessId( const int blockIdx )
{
  // If this is reader instance is serial, return Process 0
  // as the Process ID for the corresponding block.
184
  if( !this->IsParallel() )
185
    {
186
    return 0;
187
    }
188
189
190
191
192
193
194
195
196
197

  int N = this->Controller->GetNumberOfProcesses();
  return( blockIdx%N );
}

//------------------------------------------------------------------------------
bool vtkAMRBaseReader::IsBlockMine( const int blockIdx )
{
  // If this reader instance does not run in parallel, then,
  // all blocks are owned by this reader.
198
  if( !this->IsParallel() )
199
    {
200
    return true;
201
    }
202
203
204

  int myRank = this->Controller->GetLocalProcessId();
  if( myRank == this->GetBlockProcessId( blockIdx ) )
205
    {
206
    return true;
207
    }
208
209
210
  return false;
}

211
212
213
214
215
//------------------------------------------------------------------------------
void vtkAMRBaseReader::InitializeArraySelections()
{
  if( this->InitialRequest )
    {
216
217
218
    this->PointDataArraySelection->DisableAllArrays();
    this->CellDataArraySelection->DisableAllArrays();
    this->InitialRequest=false;
219
220
221
    }
}

222
223
224
225
//------------------------------------------------------------------------------
bool vtkAMRBaseReader::IsParallel( )
{
  if( this->Controller == NULL )
226
    {
227
    return false;
228
    }
229
230

  if( this->Controller->GetNumberOfProcesses() > 1 )
231
    {
232
    return true;
233
    }
234
235
236
237

  return false;
}

238
239
240
241
242
243
244

//------------------------------------------------------------------------------
int vtkAMRBaseReader::RequestInformation(
    vtkInformation *rqst,
    vtkInformationVector **inputVector,
    vtkInformationVector *outputVector )
{
245
  if( this->LoadedMetaData )
246
    {
247
    return( 1 );
248
    }
249

250
251
  vtkTimerLog::MarkStartEvent( "vtkAMRBaseReader::GenerateMetadata" );
  this->Superclass::RequestInformation( rqst, inputVector, outputVector );
252
  if( this->Metadata == NULL )
253
    {
254
    this->Metadata = vtkOverlappingAMR::New();
255
256
257
258
    vtkInformation* info = outputVector->GetInformationObject(0);
    assert( "pre: output information object is NULL" && (info != NULL) );
    this->FillMetaData( );
    info->Set( vtkCompositeDataPipeline::COMPOSITE_DATA_META_DATA(),
259
        this->Metadata );
260
    }
261
262
263
264
  vtkTimerLog::MarkStartEvent("vtkAMRBaseReader::GenerateParentChildInformation");
  this->Metadata->GenerateParentChildInformation();
  vtkTimerLog::MarkEndEvent("vtkAMRBaseReader::GenerateParentChildInformation");

265
  vtkTimerLog::MarkEndEvent( "vtkAMRBaseReader::GenerateMetadata" );
266

267
268
269
270
271
272
273
274
275
276
  std::cout << "TOTAL NUMBER OF LEVELS: " << this->Metadata->GetNumberOfLevels()
            << "\n";
  std::cout.flush();
  unsigned int levelIdx = 0;
  for( ; levelIdx < this->Metadata->GetNumberOfLevels(); ++levelIdx )
    {
    std::cout << " \tL(" << levelIdx << ") = "
              << this->Metadata->GetNumberOfDataSets( levelIdx ) << "\n";
    std::cout.flush();
    } // END for levels
277
278
279
  return 1;
}

280
//------------------------------------------------------------------------------
281
void vtkAMRBaseReader::SetupBlockRequest( vtkInformation *outInf )
282
{
283
  assert( "pre: output information is NULL" && (outInf != NULL) );
284

285
  if( outInf->Has(
George Zagaris's avatar
George Zagaris committed
286
      vtkCompositeDataPipeline::UPDATE_COMPOSITE_INDICES() ) )
287
    {
288
    assert( "Metadata should not be null" && (this->Metadata!=NULL) );
289
    this->ReadMetaData();
290

291
292
293
294
    int size =
     outInf->Length(vtkCompositeDataPipeline::UPDATE_COMPOSITE_INDICES() );
    int *indices =
      outInf->Get(vtkCompositeDataPipeline::UPDATE_COMPOSITE_INDICES() );
George Zagaris's avatar
George Zagaris committed
295

296
297
    this->BlockMap.clear();
    this->BlockMap.resize( size );
George Zagaris's avatar
George Zagaris committed
298

299
300
301
302
    for( int i=0; i < size; ++i )
      {
      this->BlockMap[ i ] = indices[ i ];
      }
303
    }
George Zagaris's avatar
George Zagaris committed
304
305
  else
    {
306
307
    this->ReadMetaData();
    this->GenerateBlockMap();
George Zagaris's avatar
George Zagaris committed
308
309
    }

310
311
312
313
314
315
316
317
318
319
320
321
}

//------------------------------------------------------------------------------
void vtkAMRBaseReader::GetAMRData(
    const int blockIdx, vtkUniformGrid *block, const char *fieldName )
{
  assert( "pre: AMR block is NULL" && (block != NULL) );
  assert( "pre: field name is NULL" && (fieldName != NULL) );

  // If caching is disabled load the data from file
  if( !this->IsCachingEnabled() )
    {
322
323
324
325
    vtkTimerLog::MarkStartEvent( "GetAMRGridDataFromFile" );
    this->GetAMRGridData( blockIdx, block, fieldName );
    vtkTimerLog::MarkEndEvent( "GetAMRGridDataFromFile" );
    return;
326
327
328
329
330
    }

  // Caching is enabled.
  // Check the cache to see if the data has already been read.
  // Otherwise, read it and cache it.
331
  if( this->Cache->HasAMRBlockCellData( blockIdx, fieldName ) )
332
    {
333
334
335
336
337
    vtkTimerLog::MarkStartEvent( "GetAMRGridDataFromCache" );
    vtkDataArray *data =
       this->Cache->GetAMRBlockCellData( blockIdx, fieldName );
    assert( "pre: cached data is NULL!" && (data != NULL) );
    vtkTimerLog::MarkEndEvent( "GetAMRGridDataFromCache" );
338

339
    block->GetCellData()->AddArray( data );
340
341
342
    }
  else
    {
343
344
345
346
347
348
349
350
    vtkTimerLog::MarkStartEvent( "GetAMRGridDataFromFile" );
    this->GetAMRGridData( blockIdx, block, fieldName );
    vtkTimerLog::MarkEndEvent( "GetAMRGridDataFromFile" );

    vtkTimerLog::MarkStartEvent( "CacheAMRData" );
    this->Cache->InsertAMRBlockCellData(
        blockIdx, block->GetCellData()->GetArray( fieldName ) );
    vtkTimerLog::MarkEndEvent( "CacheAMRData" );
351
352
353
354
355
356
357
358
359
360
361
    }

}

//------------------------------------------------------------------------------
vtkUniformGrid* vtkAMRBaseReader::GetAMRBlock( const int blockIdx )
{

  // If caching is disabled load the data from file
  if( !this->IsCachingEnabled() )
    {
362
    ++this->NumBlocksFromFile;
363
364
365
366
367
    vtkTimerLog::MarkStartEvent( "ReadAMRBlockFromFile" );
    vtkUniformGrid *gridPtr = this->GetAMRGrid( blockIdx );
    vtkTimerLog::MarkEndEvent( "ReadAMRBlockFromFile" );
    assert( "pre: grid pointer is NULL" && (gridPtr != NULL) );
    return( gridPtr );
368
369
370
371
372
    }

  // Caching is enabled.
  // Check the cache to see if the block has already been read.
  // Otherwise, read it and cache it.
373
  if( this->Cache->HasAMRBlock( blockIdx ) )
374
    {
375
    ++this->NumBlocksFromCache;
376
377
378
379
380
381
    vtkTimerLog::MarkStartEvent("ReadAMRBlockFromCache");
    vtkUniformGrid *gridPtr    = vtkUniformGrid::New();
    vtkUniformGrid *cachedGrid = this->Cache->GetAMRBlock( blockIdx );
    gridPtr->CopyStructure( cachedGrid );
    vtkTimerLog::MarkEndEvent( "ReadAMRBlockFromCache" );
    return( gridPtr );
382
383
384
    }
  else
    {
385
    ++this->NumBlocksFromFile;
386
387
388
389
390
391
392
393
394
395
396
397
    vtkTimerLog::MarkStartEvent( "ReadAMRBlockFromFile" );
    vtkUniformGrid *cachedGrid = vtkUniformGrid::New();
    vtkUniformGrid *gridPtr    = this->GetAMRGrid( blockIdx );
    assert( "pre: grid pointer is NULL" && (gridPtr != NULL) );
    vtkTimerLog::MarkEndEvent( "ReadAMRBlockFromFile" );

    vtkTimerLog::MarkStartEvent( "CacheAMRBlock" );
    cachedGrid->CopyStructure( gridPtr );
    this->Cache->InsertAMRBlock( blockIdx, cachedGrid );
    vtkTimerLog::MarkEndEvent( "CacheAMRBlock" );

    return( gridPtr );
398
399
    }

400
401
//  assert( "Code should never reach here!" && (false) );
//  return NULL;
402
403
}

404
405
//------------------------------------------------------------------------------
void vtkAMRBaseReader::LoadPointData(
George Zagaris's avatar
George Zagaris committed
406
    const int vtkNotUsed(blockIdx), vtkUniformGrid* vtkNotUsed(block) )
407
408
{
  // TODO: implement this
George Zagaris's avatar
George Zagaris committed
409
  //vtkErrorMacro( "Node-centered AMR data are not currently supported" );
410
411
412
413
414
415
416
417
418
419
420
}

//------------------------------------------------------------------------------
void vtkAMRBaseReader::LoadCellData(
    const int blockIdx, vtkUniformGrid *block )
{
  // Sanity check!
  assert( "pre: AMR block should not be NULL" && (block != NULL) );

  for( int i=0; i < this->GetNumberOfCellArrays(); ++i )
    {
421
422
423
424
425
    if( this->GetCellArrayStatus( this->GetCellArrayName( i ) ) )
      {
      this->GetAMRData( blockIdx, block, this->GetCellArrayName( i ) );
      }
    } // END for all cell arrays
426
427
}

428
//------------------------------------------------------------------------------
429
void vtkAMRBaseReader::LoadRequestedBlocks( vtkOverlappingAMR *output )
430
{
431
  assert( "pre: AMR data-structure is NULL" && (output != NULL) );
432

433
434
435
436
  //STEP 1: Gather all blocks loaded by each process
  int numBlocks = static_cast< int >( this->BlockMap.size() );
  for( int block=0; block < numBlocks; ++block )
    {
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
    int blockIdx = this->BlockMap[ block ];
    int level    = this->GetBlockLevel( blockIdx );

     // STEP 0: Get the AMR block
    vtkTimerLog::MarkStartEvent( "GetAMRBlock" );
    vtkUniformGrid *amrBlock = this->GetAMRBlock( blockIdx );
    vtkTimerLog::MarkEndEvent( "GetAMRBlock" );
    assert( "pre: AMR block is NULL" && (amrBlock != NULL) );

    // STEP 2: Load any point-data
    vtkTimerLog::MarkStartEvent( "vtkARMBaseReader::LoadPointData" );
    this->LoadPointData( blockIdx, amrBlock );
    vtkTimerLog::MarkEndEvent( "vtkAMRBaseReader::LoadPointData" );

    // STEP 3: Load any cell data
    vtkTimerLog::MarkStartEvent( "vtkAMRBaseReader::LoadCellData" );
    this->LoadCellData( blockIdx, amrBlock );
    vtkTimerLog::MarkEndEvent( "vtkAMRBaseReader::LoadCellData" );

    // STEP 4: Add dataset
457
458
459
460
461
    unsigned int metaLevel = 0;
    unsigned int metaIdx   = 0;
    this->Metadata->GetLevelAndIndex( blockIdx, metaLevel, metaIdx );

    output->SetDataSet( level,metaIdx,amrBlock );
462
    amrBlock->Delete();
463
464
    } // END for all blocks

465
466
467
468
  assert( "pre: Metadata should not be NULL" && (this->Metadata != NULL) );
  assert( "pre: NumLevels in output should be <= to NumLevels in metadata" &&
           output->GetNumberOfLevels() <= this->Metadata->GetNumberOfLevels() );

469
#ifndef NDEBUG
470
471
472
473
474
475
476
477
  unsigned int levelIdx = 0;
  for( ; levelIdx < output->GetNumberOfLevels(); ++levelIdx )
    {
    unsigned int N         = output->GetNumberOfDataSets( levelIdx );
    unsigned int Nexpected = this->Metadata->GetNumberOfDataSets( levelIdx );
    assert( "pre: NumData at level must correspond to the metadata" &&
        N <= Nexpected );
    }
478
#endif
479

480
481
482
}

//------------------------------------------------------------------------------
483
void vtkAMRBaseReader::AssignAndLoadBlocks( vtkOverlappingAMR *output )
484
485
{
  assert( "pre: AMR data-structure is NULL" && (output != NULL) );
486

487
488
489
490
491
  // Initialize counter of the number of blocks at each level.
  // This counter is used to compute the block index w.r.t. the
  // hierarchical box data-structure. Note that then number of blocks
  // can change based on user constraints, e.g., the number of levels
  // visible.
492
  std::vector< int > idxcounter;
493
  idxcounter.resize(this->GetNumberOfLevels()+1, 0);
494
495
496
497
498
499
500

  // Find the number of blocks to be processed. BlockMap.size()
  // has all the blocks that are to be processesed and may be
  // less than or equal to this->GetNumberOfBlocks(), i.e., the
  // total number of blocks.
  int numBlocks = static_cast< int >( this->BlockMap.size() );
  for( int block=0; block < numBlocks; ++block )
501
    {
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
    int blockIdx = this->BlockMap[ block ];
    int level    = this->GetBlockLevel( blockIdx );

    if( this->IsBlockMine(block) )
      {
      // STEP 0: Get the AMR block
      vtkTimerLog::MarkStartEvent( "GetAMRBlock" );
      vtkUniformGrid *amrBlock = this->GetAMRBlock( blockIdx );
      vtkTimerLog::MarkEndEvent( "GetAMRBlock" );
      assert( "pre: AMR block is NULL" && (amrBlock != NULL) );

      // STEP 2: Load any point-data
      vtkTimerLog::MarkStartEvent( "vtkARMBaseReader::LoadPointData" );
      this->LoadPointData( blockIdx, amrBlock );
      vtkTimerLog::MarkEndEvent( "vtkAMRBaseReader::LoadPointData" );
517

518
519
520
521
      // STEP 3: Load any cell data
      vtkTimerLog::MarkStartEvent( "vtkAMRBaseReader::LoadCellData" );
      this->LoadCellData( blockIdx, amrBlock );
      vtkTimerLog::MarkEndEvent( "vtkAMRBaseReader::LoadCellData" );
522

523
524
525
526
527
528
529
530
531
532
      // STEP 4: Add dataset
      output->SetDataSet( level,idxcounter[level],amrBlock );
      amrBlock->Delete();
      idxcounter[level]++;
      } // END if the block belongs to this process
    else
      {
      output->SetDataSet( level, idxcounter[level], NULL );
      idxcounter[level]++;
      }
533
    } // END for all blocks
534
535
536
537
538
539
540
541
542
}

//------------------------------------------------------------------------------
int vtkAMRBaseReader::RequestData(
        vtkInformation* vtkNotUsed(request),
        vtkInformationVector** vtkNotUsed(inputVector),
        vtkInformationVector* outputVector )
{
  vtkTimerLog::MarkStartEvent( "vtkAMRBaseReader::RqstData" );
543
544
  this->NumBlocksFromCache = 0;
  this->NumBlocksFromFile  = 0;
545

George Zagaris's avatar
George Zagaris committed
546
  vtkInformation    *outInf = outputVector->GetInformationObject( 0 );
547
548
  vtkOverlappingAMR *output =
    vtkOverlappingAMR::SafeDownCast(
549
550
     outInf->Get( vtkDataObject::DATA_OBJECT() ) );
  assert( "pre: output AMR dataset is NULL" && ( output != NULL ) );
551

552
553
554
555
  // Setup the block request
  vtkTimerLog::MarkStartEvent( "vtkAMRBaseReader::SetupBlockRequest" );
  this->SetupBlockRequest( outInf );
  vtkTimerLog::MarkEndEvent( "vtkAMRBaseReader::SetupBlockRequest" );
556

557
  if( outInf->Has( vtkCompositeDataPipeline::LOAD_REQUESTED_BLOCKS() ) )
558
    {
559
    this->LoadRequestedBlocks( output );
560
561
562
    }
  else
    {
563
564
565
566
567
568
569
570
571
572
573
574
575
576
    this->AssignAndLoadBlocks( output );

    // Generate all the AMR metadata & the visibility arrays
    vtkTimerLog::MarkStartEvent( "AMRUtilities::GenerateMetaData" );
    vtkAMRUtilities::GenerateMetaData( output, this->Controller );
    vtkTimerLog::MarkEndEvent( "AMRUtilities::GenerateMetaData" );

    //If there is a downstream module, do not generate visibility arrays here.
    if(!outInf->Has( vtkCompositeDataPipeline::UPDATE_COMPOSITE_INDICES() ) )
      {
      vtkTimerLog::MarkStartEvent( "AMR::GenerateVisibilityArrays" );
      output->GenerateVisibilityArrays();
      vtkTimerLog::MarkEndEvent( "AMR::GenerateVisibilityArrays" );
      }
577
    }
578
579
580

  // If this instance of the reader is not parallel, block until all processes
  // read their blocks.
581
  if( this->IsParallel() )
582
    {
583
    this->Controller->Barrier();
584
    }
585

586
587
  outInf = NULL;
  output = NULL;
588
589

  vtkTimerLog::MarkEndEvent( "vtkAMRBaseReader::RqstData" );
George Zagaris's avatar
George Zagaris committed
590

591
592
  return 1;
}