vtkAMRResampleFilter.cxx 42.2 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
/*=========================================================================

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

 =========================================================================*/

16
#include "vtkAMRResampleFilter.h"
George Zagaris's avatar
George Zagaris committed
17
#include "vtkAMRBox.h"
18
#include "vtkObjectFactory.h"
19
#include "vtkStreamingDemandDrivenPipeline.h"
20
21
22
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkMultiBlockDataSet.h"
23
#include "vtkOverlappingAMR.h"
24
25
26
#include "vtkMultiProcessController.h"
#include "vtkUniformGrid.h"
#include "vtkIndent.h"
27
28
29
30
#include "vtkAMRUtilities.h"
#include "vtkBoundingBox.h"
#include "vtkMath.h"
#include "vtkCompositeDataPipeline.h"
George Zagaris's avatar
George Zagaris committed
31
32
33
#include "vtkFieldData.h"
#include "vtkCellData.h"
#include "vtkPointData.h"
34
#include "vtkCell.h"
George Zagaris's avatar
George Zagaris committed
35
#include "vtkXMLImageDataWriter.h"
36
#include "vtkExtentRCBPartitioner.h"
37
#include "vtkUniformGridPartitioner.h"
38
#include "vtkDataArray.h"
39

40
41
#include "vtkTimerLog.h"

42
#include <cassert>
43
#include <cmath>
George Zagaris's avatar
George Zagaris committed
44
#include <sstream>
George Zagaris's avatar
George Zagaris committed
45
#include <cmath>
46
#include <algorithm>
47

48
49


50
vtkStandardNewMacro( vtkAMRResampleFilter );
51
52

//-----------------------------------------------------------------------------
53
vtkAMRResampleFilter::vtkAMRResampleFilter()
54
{
55
  this->TransferToNodes      = 1;
56
57
  this->DemandDrivenMode     = 0;
  this->NumberOfPartitions   = 1;
George Zagaris's avatar
George Zagaris committed
58
  this->LevelOfResolution    = 0;
59
  this->AMRMetaData          = NULL;
George Zagaris's avatar
George Zagaris committed
60
  this->NumberOfSamples[0]   = this->NumberOfSamples[1] = this->NumberOfSamples[2] = 10;
61
  this->Controller           = vtkMultiProcessController::GetGlobalController();
62
  this->ROI                  = vtkMultiBlockDataSet::New();
63
64
65
66
67
68

  for( int i=0; i < 3; ++i )
    {
    this->Min[i] = 0.;
    this->Max[i] = 1.;
    }
69
70
  this->SetNumberOfInputPorts( 1 );
  this->SetNumberOfOutputPorts( 1 );
71
72
  this->UseBiasVector = false;
  this->BiasVector[0] = this->BiasVector[1] = this->BiasVector[2] = 0.0;
73
74
75
}

//-----------------------------------------------------------------------------
76
vtkAMRResampleFilter::~vtkAMRResampleFilter()
77
{
78
  this->BlocksToLoad.clear();
79
80

  if( this->ROI != NULL )
81
    {
82
    this->ROI->Delete();
83
    }
84
  this->ROI = NULL;
George Zagaris's avatar
George Zagaris committed
85

86
87
88
89
90
//  if( this->AMRMetaData != NULL )
//    {
//    this->AMRMetaData->Delete();
//    }
//  this->AMRMetaData = NULL;
91
92
93
}

//-----------------------------------------------------------------------------
94
void vtkAMRResampleFilter::PrintSelf( std::ostream &oss, vtkIndent indent )
95
96
97
98
99
{
  this->Superclass::PrintSelf( oss, indent );
}

//-----------------------------------------------------------------------------
100
int vtkAMRResampleFilter::FillInputPortInformation(
101
102
103
104
    int vtkNotUsed(port), vtkInformation *info )
{
  assert( "pre: information object is NULL" && (info != NULL) );
  info->Set(
105
   vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkOverlappingAMR" );
106
107
108
109
  return 1;
}

//-----------------------------------------------------------------------------
110
int vtkAMRResampleFilter::FillOutputPortInformation(
111
112
113
    int vtkNotUsed(port), vtkInformation *info )
{
  assert( "pre: information object is NULL" && (info != NULL) );
114
  info->Set(vtkDataObject::DATA_TYPE_NAME(),"vtkMultiBlockDataSet");
115
116
117
118
  return 1;
}

//-----------------------------------------------------------------------------
119
int vtkAMRResampleFilter::RequestUpdateExtent(
120
    vtkInformation*, vtkInformationVector **inputVector,
George Zagaris's avatar
George Zagaris committed
121
    vtkInformationVector* vtkNotUsed(outputVector) )
122
{
George Zagaris's avatar
George Zagaris committed
123
  assert( "pre: inputVector is NULL" && (inputVector != NULL) );
124
125
126
  vtkInformation *info = inputVector[0]->GetInformationObject(0);
  assert( "pre: info is NULL" && (info != NULL) );

127
128
129
130
  if( this->DemandDrivenMode == 1 )
    {
    // Tell reader to load all requested blocks.
    info->Set( vtkCompositeDataPipeline::LOAD_REQUESTED_BLOCKS(), 1 );
131

132
133
134
    // Tell reader which blocks this process requires
    info->Set(
        vtkCompositeDataPipeline::UPDATE_COMPOSITE_INDICES(),
135
        &this->BlocksToLoad[0], static_cast<int>(this->BlocksToLoad.size()));
136
    }
137
138
139
140
  return 1;
}

//-----------------------------------------------------------------------------
141
int vtkAMRResampleFilter::RequestInformation(
142
    vtkInformation* vtkNotUsed(rqst),
143
    vtkInformationVector** inputVector,
George Zagaris's avatar
George Zagaris committed
144
    vtkInformationVector* vtkNotUsed(outputVector) )
145
{
146

George Zagaris's avatar
George Zagaris committed
147
148
  assert( "pre: inputVector is NULL" && (inputVector != NULL) );

149
150
151
  vtkInformation *input = inputVector[0]->GetInformationObject( 0 );
  assert( "pre: input is NULL" && (input != NULL)  );

152
153
  if( this->DemandDrivenMode == 1 &&
      input->Has(vtkCompositeDataPipeline::COMPOSITE_DATA_META_DATA() ) )
154
    {
155
156
//    this->AMRMetaData = vtkOverlappingAMR::New();
    this->AMRMetaData =
157
    vtkOverlappingAMR::SafeDownCast(
158
      input->Get( vtkCompositeDataPipeline::COMPOSITE_DATA_META_DATA() ) );
159

160

161
162
    // Get Region
    double h[3];
163
    this->ComputeAndAdjustRegionParameters( this->AMRMetaData, h );
164
    this->GetRegion( h );
165

166
    // Compute which blocks to load
167
    this->ComputeAMRBlocksToLoad( this->AMRMetaData );
168
    }
169

170
171
// Don't we need to call this->Modified() here?
// this->Modified();
172
 return 1;
173
174
175
}

//-----------------------------------------------------------------------------
176
int vtkAMRResampleFilter::RequestData(
George Zagaris's avatar
George Zagaris committed
177
    vtkInformation* vtkNotUsed(rqst), vtkInformationVector** inputVector,
178
    vtkInformationVector* outputVector )
179
{
180
181
  cerr << "Running Resampler\n";

182
183
184
  // STEP 0: Get input object
  vtkInformation *input = inputVector[0]->GetInformationObject( 0 );
  assert( "pre: Null information object!" && (input != NULL) );
185
186
  vtkOverlappingAMR *amrds=
     vtkOverlappingAMR::SafeDownCast(
187
188
189
      input->Get(vtkDataObject::DATA_OBJECT()));
  assert( "pre: input AMR dataset is NULL" && (amrds != NULL) );

190
191
192
193
194
195
196
197
198
199
  // STEP 1: Get output object
   vtkInformation *output = outputVector->GetInformationObject( 0 );
   assert( "pre: Null output information object!" && (output != NULL) );

   vtkMultiBlockDataSet *mbds =
      vtkMultiBlockDataSet::SafeDownCast(
          output->Get(vtkDataObject::DATA_OBJECT() ) );
   assert( "pre: ouput grid is NULL" && (mbds != NULL) );

  // STEP 2: Get Metadata
200
  if( this->DemandDrivenMode == 1 )
201
    {
202
203
204
    assert( "pre: Metadata must have been populated in RqstInfo" &&
            (this->AMRMetaData != NULL) );
    this->ExtractRegion( amrds, mbds, this->AMRMetaData );
205
206
207
    }
  else
    {
208
209
210
211
212
    // GetRegion
    double h[3];
    this->ComputeAndAdjustRegionParameters(amrds, h );
    this->GetRegion( h );
    this->ExtractRegion( amrds, mbds, amrds);
213
    }
214

215
216
  return 1;
}
217

George Zagaris's avatar
George Zagaris committed
218
//-----------------------------------------------------------------------------
219
bool vtkAMRResampleFilter::FoundDonor(
220
    double q[3],vtkUniformGrid *&donorGrid,int &cellIdx)
George Zagaris's avatar
George Zagaris committed
221
222
{
  assert( "pre: donor grid is NULL" && (donorGrid != NULL) );
223
224
  double gbounds[6];
  // Lets do a trival spatial check
225
  this->NumberOfBlocksTested++;
226
227
228
229
230
231
232
  donorGrid->GetBounds(gbounds);
  if ((q[0] < gbounds[0]) || (q[0] > gbounds[1]) ||
      (q[1] < gbounds[2]) || (q[1] > gbounds[3]) || 
      (q[2] < gbounds[4]) || (q[2] > gbounds[5]))
    {
    return false;
    }
George Zagaris's avatar
George Zagaris committed
233
234
235
236
237
  int ijk[3];
  double pcoords[3];
  int status = donorGrid->ComputeStructuredCoordinates( q, ijk, pcoords );
  if( status == 1 )
    {
238
239
    cellIdx=vtkStructuredData::ComputeCellId(donorGrid->GetDimensions(),ijk);
    return true;
George Zagaris's avatar
George Zagaris committed
240
241
242
243
244
    }
  return false;
}

//-----------------------------------------------------------------------------
245
void vtkAMRResampleFilter::InitializeFields(
George Zagaris's avatar
George Zagaris committed
246
247
248
249
250
251
252
      vtkFieldData *f, vtkIdType size, vtkCellData *src )
{
  assert( "pre: field data is NULL!" && (f != NULL) );
  assert( "pre: source cell data is NULL" && (src != NULL) );

  for( int arrayIdx=0; arrayIdx < src->GetNumberOfArrays(); ++arrayIdx )
    {
253
254
255
    int dataType        = src->GetArray( arrayIdx )->GetDataType();
    vtkDataArray *array = vtkDataArray::CreateDataArray( dataType );
    assert( "pre: failed to create array!" && (array != NULL) );
George Zagaris's avatar
George Zagaris committed
256

257
258
259
260
261
262
    array->SetName( src->GetArray(arrayIdx)->GetName() );
    array->SetNumberOfComponents(
             src->GetArray(arrayIdx)->GetNumberOfComponents() );
    array->SetNumberOfTuples( size );
    assert( "post: array size mismatch" &&
            (array->GetNumberOfTuples() == size) );
George Zagaris's avatar
George Zagaris committed
263

264
265
    f->AddArray( array );
    array->Delete();
266

267
268
269
270
271
272
273
//      int myIndex = -1;
//      vtkDataArray *arrayPtr = f->GetArray(
//          src->GetArray(arrayIdx)->GetName(), myIndex );
//      assert( "post: array index mismatch" && (myIndex == arrayIdx) );
//      assert( "post: array size mismatch" &&
//              (arrayPtr->GetNumberOfTuples()==size) );

274
275
    assert( "post: array size mismatch" &&
            (f->GetArray( arrayIdx)->GetNumberOfTuples() == size) );
George Zagaris's avatar
George Zagaris committed
276
277
278
279
280
    } // END for all arrays

}

//-----------------------------------------------------------------------------
281
void vtkAMRResampleFilter::CopyData(
George Zagaris's avatar
George Zagaris committed
282
283
284
285
286
287
288
289
290
291
292
    vtkFieldData *target, vtkIdType targetIdx,
    vtkCellData *src, vtkIdType srcIdx )
{
  assert( "pre: target field data is NULL" && (target != NULL) );
  assert( "pre: source field data is NULL" && (src != NULL) );
  assert( "pre: number of arrays does not match" &&
           (target->GetNumberOfArrays() == src->GetNumberOfArrays() ) );

  int arrayIdx = 0;
  for( ; arrayIdx < src->GetNumberOfArrays(); ++arrayIdx )
    {
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
    vtkDataArray *targetArray = target->GetArray( arrayIdx );
    vtkDataArray *srcArray    = src->GetArray( arrayIdx );
    assert( "pre: target array is NULL!" && (targetArray != NULL) );
    assert( "pre: source array is NULL!" && (srcArray != NULL) );
    assert( "pre: targer/source array number of components mismatch!" &&
            (targetArray->GetNumberOfComponents()==
             srcArray->GetNumberOfComponents() ) );
    assert( "pre: target/source array names mismatch!" &&
            (strcmp(targetArray->GetName(),srcArray->GetName()) == 0) );
    assert( "pre: source index is out-of-bounds" &&
            (srcIdx >=0) &&
            (srcIdx < srcArray->GetNumberOfTuples() ) );
    assert( "pre: target index is out-of-bounds" &&
            (targetIdx >= 0) &&
            (targetIdx < targetArray->GetNumberOfTuples() ) );

    int c=0;
    for( ; c < srcArray->GetNumberOfComponents(); ++c )
      {
      double f = srcArray->GetComponent( srcIdx, c );
      targetArray->SetComponent( targetIdx, c, f );
      } // END for all componenents
George Zagaris's avatar
George Zagaris committed
315
316
317
318

    } // END for all arrays
}

319
//-----------------------------------------------------------------------------
320
void vtkAMRResampleFilter::ComputeCellCentroid(
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
    vtkUniformGrid *g, const vtkIdType cellIdx, double c[3] )
{
  assert( "pre: uniform grid is NULL" && (g != NULL) );
  assert( "pre: centroid is NULL" && (c != NULL) );
  assert( "pre: cell index out-of-bounds" &&
           (cellIdx >= 0) && (cellIdx < g->GetNumberOfCells()) );


  vtkCell *myCell = g->GetCell( cellIdx );
  assert( "post: cell is NULL!" && (myCell != NULL) );

  double pc[3]; // the parametric center
  double *weights = new double[ myCell->GetNumberOfPoints() ];
  assert( "post: weights vector is NULL" && (weights != NULL) );

  int subId = myCell->GetParametricCenter( pc );
  myCell->EvaluateLocation( subId, pc, c, weights );
  delete [] weights;
}

George Zagaris's avatar
George Zagaris committed
341
//-----------------------------------------------------------------------------
342
void vtkAMRResampleFilter::TransferToCellCenters(
343
        vtkUniformGrid *g, vtkOverlappingAMR *amrds )
George Zagaris's avatar
George Zagaris committed
344
345
346
347
{
  assert( "pre: uniform grid is NULL" && (g != NULL) );
  assert( "pre: AMR data-strucutre is NULL" && (amrds != NULL) );

348
  // STEP 0: Get the first block so that we know the arrays
349
350
351
//  vtkUniformGrid *refGrid = amrds->GetDataSet(0,0);
//  assert( "pre: Block(0,0) is NULL!" && (refGrid != NULL) );
  vtkUniformGrid *refGrid = this->GetReferenceGrid( amrds );
352
353
354
355
356
357
358
359
360
361
362
363
364

  // STEP 1: Get the cell-data of the reference grid
  vtkCellData *CD = refGrid->GetCellData();
  assert( "pre: Donor CellData is NULL!" && (CD != NULL)  );

  // STEP 2: Get the cell data of the resampled grid
  vtkCellData *fieldData = g->GetCellData();
  assert( "pre: Target PointData is NULL!" && (fieldData != NULL) );

  // STEP 3: Initialize the fields on the resampled grid
  this->InitializeFields( fieldData, g->GetNumberOfCells(), CD );

  if(fieldData->GetNumberOfArrays() == 0)
365
366
367
    {
    return;
    }
368

369
370
371
  // TODO: this is a very naive implementation and should be optimized. However,
  // mostly this filter is used to transfer the solution to the grid nodes and
  // not on the cell nodes.
372
373
374
  vtkIdType cellIdx = 0;
  for( ; cellIdx < g->GetNumberOfCells(); ++cellIdx )
    {
375
376
    double qPoint[3];
    this->ComputeCellCentroid( g, cellIdx, qPoint );
377

378
379
380
381
382
    unsigned int level=0;
    for( ; level < amrds->GetNumberOfDataSets( level ); ++level )
      {
      unsigned int dataIdx = 0;
      for( ; dataIdx < amrds->GetNumberOfDataSets( level ); ++dataIdx )
383
        {
384
385
386
387
388
          int donorCellIdx = -1;
          vtkUniformGrid *donorGrid = amrds->GetDataSet(level,dataIdx);
          if( (donorGrid!=NULL) &&
              this->FoundDonor(qPoint,donorGrid,donorCellIdx) )
            {
389
390
391
392
393
            assert( "pre: donorCellIdx is invalid" &&
                    (donorCellIdx >= 0) &&
                    (donorCellIdx < donorGrid->GetNumberOfCells()) );
            CD = donorGrid->GetCellData();
            this->CopyData( fieldData, cellIdx, CD, donorCellIdx );
394
            } // END if
395
396
        } // END for all datasets
      } // END for all levels
397
    } // END for all cells
George Zagaris's avatar
George Zagaris committed
398
399
}

400
401
//-----------------------------------------------------------------------------
void vtkAMRResampleFilter::SearchForDonorGridAtLevel(
402
    double q[3], vtkOverlappingAMR *amrds,
403
    unsigned int level, unsigned int donorGridId, vtkUniformGrid *&donorGrid,
404
    int &donorCellIdx )
405
406
{
  assert( "pre: AMR dataset is NULL" && (amrds != NULL) );
407
  this->NumberOfBlocksTestedForLevel = 0;
408
409
410
411
412
  std::ostringstream oss;
  oss << "SearchLevel-" << level;

  vtkTimerLog::MarkStartEvent( oss.str().c_str() );

413
  for(donorGridId = 0; donorGridId < amrds->GetNumberOfDataSets(level); ++donorGridId )
414
415
    {
    donorCellIdx = -1;
416
    donorGrid    = amrds->GetDataSet(level,donorGridId);
417
    this->NumberOfBlocksTestedForLevel++;
418
419
420
421
422
423
    if( (donorGrid!=NULL) &&
       this->FoundDonor(q,donorGrid,donorCellIdx) )
     {
     assert( "pre: donorCellIdx is invalid" &&
             (donorCellIdx >= 0) &&
             (donorCellIdx < donorGrid->GetNumberOfCells()) );
424
     vtkTimerLog::MarkEndEvent( oss.str().c_str() );
425
426
427
428
429
     return;
     } // END if

    } // END for all data at level

430

431
432
433
  // No suitable grid is found at the requested level, set donorGrid to NULL
  // to indicate that to the caller.
  donorGrid = NULL;
434
  vtkTimerLog::MarkEndEvent( oss.str().c_str() );
435
436
}

437
438
439
//-----------------------------------------------------------------------------
int vtkAMRResampleFilter::ProbeGridPointInAMR(
    double q[3], vtkUniformGrid *&donorGrid, unsigned int &donorLevel,
440
    vtkOverlappingAMR *amrds, unsigned int maxLevel )
441
442
443
444
445
446
{
  assert( "pre: AMR dataset is NULL" && amrds != NULL );

  vtkUniformGrid *currentGrid = NULL;
  int currentCellIdx          = -1;
  int donorCellIdx            = -1;
George Zagaris's avatar
George Zagaris committed
447
  unsigned int currentLevel = 0;
Bob Obara's avatar
Bob Obara committed
448
  bool hadDonorGrid = false;
George Zagaris's avatar
George Zagaris committed
449
450
  unsigned int donorGridId = 0;
  unsigned int currentGridId = 0;
451
452
453
454

  // STEP 0: Check the previously cached donor-grid
  if( donorGrid != NULL )
    {
Bob Obara's avatar
Bob Obara committed
455
    hadDonorGrid = true;
456
    this->NumberOfBlocksTested++;
457
458
459
    if(!this->FoundDonor( q, donorGrid, donorCellIdx ) )
      {
      // Lets see if the point is contained by a grid at the same donar level
460
461
      this->SearchForDonorGridAtLevel(q,amrds,donorLevel,donorGridId,
                                      donorGrid,donorCellIdx);
462
      this->NumberOfBlocksTested += this->NumberOfBlocksTestedForLevel;
463
464
465
466
      }
    // If donorGrid is still not NULL then we found the grid and potential starting
    // level
    if (donorGrid != NULL)
467
468
      {
      assert( "pre: donorCellIdx is invalid" &&
469
470
              (donorCellIdx >= 0) && (donorCellIdx < donorGrid->GetNumberOfCells()) );
      
471
      this->NumberOfTimesFoundOnDonorLevel++;
472
473
474
475
      // Lets see if the cell is blanked - if it isn't then we found the highest
      // resolution grid that contains the point
      if (donorGrid->IsCellVisible(donorCellIdx))
        {
Bob Obara's avatar
Bob Obara committed
476
        //return donorCellIdx;
477
        }
478
479
480
481

      // Initialize values for step 1 s.t. that the search will start from the
      // current donorLevel
      currentGrid    = donorGrid;
482
      currentGridId = donorGridId;
483
      currentCellIdx = donorCellIdx;
Bob Obara's avatar
Bob Obara committed
484
      currentLevel = donorLevel+1;
485
      }
Bob Obara's avatar
Bob Obara committed
486
487
488
489
490
491
492
493
494
    else if (donorLevel == 0)
      {
      //if we are here then the point is not contained in any of the level 0
      // blocks!
      this->NumberOfFailedPoints++;
      donorGrid    = NULL;
      donorLevel = 0;
      return -1;
      }
495
496
    else
      {
497
498
499
500
      // If we are here then we know the point is not on the donor level
      // and therefore not contained in any of the more refined levels -
      // Base on the assumption of overlapping AMR
      assert("pre:Donor Level is 0" && donorLevel != 0);
501
502
      // Initialize values for step 1 s.t. the search will start from level 0.
      donorGrid  = NULL;
503
504
      maxLevel = donorLevel;
      donorLevel = 0;
Bob Obara's avatar
Bob Obara committed
505
      currentLevel = 0;
506
507
508
      }
    }

Bob Obara's avatar
Bob Obara committed
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
  // If we didn't have an initial donor grid or if we still have one
  // we need to test higher res grids
  int startLevel, endLevel;
  int incLevel;
  if (!((donorGrid == NULL) && hadDonorGrid))
    {
    startLevel = currentLevel;
    endLevel = maxLevel;
    incLevel = 1;
    }
  else
    {
    startLevel = maxLevel-1;
    endLevel = -1;
    incLevel = -1;
    }
525
  // STEP 1: Search in the AMR hierarchy for the donor-grid
Bob Obara's avatar
Bob Obara committed
526
  for( int level=startLevel; level != endLevel; level += incLevel )
527
    {
528
529
530
531
532
533
534
535
    if (incLevel == 1)
      {
      this->NumberOfTimesLevelUp++;
      }
    else
      {
      this->NumberOfTimesLevelDown++;
      }
536
    this->SearchForDonorGridAtLevel(q,amrds,level,donorGridId,donorGrid,donorCellIdx);
537
    this->NumberOfBlocksTested += this->NumberOfBlocksTestedForLevel;
538
539
    if( donorGrid != NULL )
      {
540
541
542
543
544
545
546
      donorLevel = level;
      // if we are going from fine to coarse then we can stop the search
      if (incLevel == -1)
        {
        return donorCellIdx;
        }

547
548
      // Lets see if this is the highest resolution grid that contains the 
      // point
549
550
      if (donorGrid->IsCellVisible(donorCellIdx))
        {
Bob Obara's avatar
Bob Obara committed
551
        //return donorCellIdx;
552
        }
553
554
555
556
557
      // we found a grid that contains the point at level l, let's store it
      // here temporatily in case there is a grid at a higher resolution that
      // we need to use.
      currentGrid    = donorGrid;
      currentCellIdx = donorCellIdx;
Bob Obara's avatar
Bob Obara committed
558
      currentLevel = level;
559
      currentGridId = donorGridId;
560
561
562
563
564
      }
    else if( currentGrid != NULL )
      {
      // we did not find the point at a higher res, but, we did find at a lower
      // resolution, so we will use the solution we found previously
565
      // THIS SHOULD NOW NOT HAPPEN!!
Bob Obara's avatar
Bob Obara committed
566
      //vtkErrorMacro("Could not find point in an unblanked cell.");
567
      this->NumberOfBlocksVisSkipped +=  this->NumberOfBlocksTestedForLevel;
568
569
      donorGrid    = currentGrid;
      donorCellIdx = currentCellIdx;
Bob Obara's avatar
Bob Obara committed
570
      donorLevel = currentLevel;
571
      donorGridId = currentGridId;
572
573
574
575
576
577
      break;
      }
    else
      {
      // we are not able to find a grid/cell that contains the query point, in
      // this case we will just return.
578
      this->NumberOfFailedPoints++;
579
580
      donorCellIdx = -1;
      donorGrid    = NULL;
Bob Obara's avatar
Bob Obara committed
581
      donorLevel = 0;
582
583
584
585
586
587
      break;
      }
    } // END for all levels
  return( donorCellIdx );
}

588
589
//-----------------------------------------------------------------------------
void vtkAMRResampleFilter::SearchGridAncestors(double q[3], 
590
                                               vtkOverlappingAMR *amrds,
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
                                               unsigned int &level,
                                               unsigned int &gridId,
                                               vtkUniformGrid *&grid,
                                               int &cellId )
{
  assert( "pre: AMR dataset is NULL" && (amrds != NULL) );
  unsigned int *parents, plevel;
  for (; level > 0; --level)
    {
    ++this->NumberOfTimesLevelUp;
    // Get the parents of the grid
    parents = amrds->GetParents(level, gridId);
    plevel = level - 1;
    // There should be at least 1 parent
    assert( "Found non-level 0 grid with no parents" && (parents != NULL) && (parents[0] > 0) );
    if (parents[0] > 1)
      {
      vtkDebugMacro( "Number of parents: " << parents[0] << " - Only processing 1 route");
      }
    gridId = parents[1];
    grid  = amrds->GetDataSet(plevel,gridId);
    if (this->FoundDonor( q, grid, cellId ))
      {
      level = plevel;
      return;
      }
    }
  // If we are here then we could not find an ancestor
  grid = NULL;
  cellId = -1;
}
//-----------------------------------------------------------------------------
void vtkAMRResampleFilter::SearchGridDecendants(double q[3], 
624
                                                vtkOverlappingAMR *amrds,
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
                                                unsigned int maxLevel,
                                                unsigned int &level,
                                                unsigned int &gridId,
                                                vtkUniformGrid *&grid,
                                                int &cellId)
{
  assert( "pre: AMR dataset is NULL" && (amrds != NULL) );
  unsigned int *children, clevel, n, i;
  vtkUniformGrid *childGrid;
  for (; level < maxLevel-1; ++level)
    {
    // Get the children of the grid
    children = amrds->GetChildren(level, gridId);
    clevel = level + 1;
    // If there are no children then we found the grid!
    if (children == NULL)
      {
      return;
      }
    n = children[0];
    for (i = 1; i <= n; ++i)
      {
      childGrid  = amrds->GetDataSet(clevel,children[i]);
      if (this->FoundDonor( q, childGrid, cellId ))
        {
        // We found a decendant so stop searching the
        // children and can instead search that grid's
        // children
        gridId = children[i];
        grid  = childGrid;
        ++this->NumberOfTimesLevelDown;
        break;
        }
      }
    if (i > n)
      {
      // We tested some children that we didn't need to if
      // we had visibility info
      this->NumberOfBlocksVisSkipped += n;
      // If we are here then no child contains the point
      // so don't search any further
      return;
      }
    }
}
//-----------------------------------------------------------------------------
int vtkAMRResampleFilter::
ProbeGridPointInAMRGraph(double q[3], vtkUniformGrid *&donorGrid,
                         unsigned int &donorLevel,  unsigned int &donorGridId, 
674
                         vtkOverlappingAMR *amrds, unsigned int maxLevel)
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
{
  assert( "pre: AMR dataset is NULL" && amrds != NULL );

  int donorCellIdx = -1;

  // STEP 0: Check the previously cached donor-grid
  if( donorGrid != NULL )
    {
    if(!this->FoundDonor( q, donorGrid, donorCellIdx ) )
      {
      // Lets find the grid's ancestor that contains the point
      this->SearchGridAncestors(q,amrds,donorLevel,donorGridId,donorGrid,donorCellIdx);
      }
    else
      {
      ++this->NumberOfTimesFoundOnDonorLevel;
      }
    // if the point is not contained in an ancestor then lets just assume its on level
    // 0 which is the default
    }

  // If there is no initial donor grid then search level 0
  if (donorGrid == NULL)
    {
    this->SearchForDonorGridAtLevel(q,amrds,0,donorGridId,donorGrid,donorCellIdx);
    // If we still can't find a grid then the point is not contained in the
    // AMR Data
    if (donorGrid == NULL)
      {
      this->NumberOfFailedPoints++;
      donorLevel = 0;
      return -1;
      }
    }

  // Now search the decendants of the donor grid
  this->SearchGridDecendants(q,amrds,maxLevel,donorLevel,donorGridId,donorGrid,donorCellIdx);
  return( donorCellIdx );
}

George Zagaris's avatar
George Zagaris committed
715
//-----------------------------------------------------------------------------
716
void vtkAMRResampleFilter::TransferToGridNodes(
717
    vtkUniformGrid *g, vtkOverlappingAMR *amrds )
George Zagaris's avatar
George Zagaris committed
718
{
719
720
721
722
723
724
725
  this->NumberOfBlocksTested = 0;
  this->NumberOfBlocksVisSkipped = 0;
  this->NumberOfTimesFoundOnDonorLevel = 0;
  this->NumberOfTimesLevelUp = 0;
  this->NumberOfTimesLevelDown = 0;
  this->NumberOfFailedPoints = 0;
  this->AverageLevel = 0.0;
George Zagaris's avatar
George Zagaris committed
726
727
728
  assert( "pre: uniform grid is NULL" && (g != NULL) );
  assert( "pre: AMR data-structure is NULL" && (amrds != NULL) );

729
  // STEP 0: Initialize the fields on the grid
730
  vtkUniformGrid *refGrid = this->GetReferenceGrid( amrds );
731
732

  vtkCellData *CD = refGrid->GetCellData();
George Zagaris's avatar
George Zagaris committed
733
734
735
736
737
  assert( "pre: Donor CellData is NULL!" && (CD != NULL)  );

  vtkPointData *PD = g->GetPointData();
  assert( "pre: Target PointData is NULL!" && (PD != NULL) );

738
  // STEP 0: Initialize the fields on the grid
George Zagaris's avatar
George Zagaris committed
739
740
  this->InitializeFields( PD, g->GetNumberOfPoints(), CD );

741
742
  // STEP 1: If no arrays are selected, there is no need to interpolate
  // anything on the grid, just return
743
  if(PD->GetNumberOfArrays() == 0)
744
    {
745
    return;
746
    }
747

748
749
750
751
752
753
754
755
756
757
  // STEP 2: Fix the maximum level at which the search algorithm will operate
  unsigned int maxLevelToLoad = 0;
  if( this->LevelOfResolution < static_cast<int>(amrds->GetNumberOfLevels()))
    {
    maxLevelToLoad = this->LevelOfResolution+1;
    }
  else
    {
    maxLevelToLoad = amrds->GetNumberOfLevels();
    }
758

759
760
761
762
  // STEP 3: Loop through all the points and find the donors.
  int numPoints           = 0;
  unsigned int donorLevel = 0;
  vtkUniformGrid *donorGrid = NULL;
763
764
765
766
767
768
  unsigned int donorGridId = 0;
  double qPoint[3];
  vtkIdType pIdx;
  int donorCellIdx;
  // Do we have parent/child meta information
  if (this->AMRMetaData)
George Zagaris's avatar
George Zagaris committed
769
    {
770
    for(pIdx = 0; pIdx < g->GetNumberOfPoints(); ++pIdx )
771
      {
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
      g->GetPoint( pIdx, qPoint );
      donorCellIdx =
        this->ProbeGridPointInAMRGraph(qPoint,donorGrid, 
                                       donorLevel, donorGridId,
                                       this->AMRMetaData,maxLevelToLoad );
      if( donorCellIdx != -1 )
        {
        vtkUniformGrid *amrGrid = amrds->GetDataSet(donorLevel,donorGridId);
        this->AverageLevel += donorLevel;
        assert(donorGrid != NULL);
        CD = amrGrid->GetCellData();
        this->CopyData( PD, pIdx, CD, donorCellIdx );
        }
      else
        {
        // Point is outside the domain, blank it
        ++ numPoints;
        g->BlankPoint( pIdx );
        }
      } // END for all grid nodes
    }
  else
    {
    for(pIdx = 0; pIdx < g->GetNumberOfPoints(); ++pIdx )
796
      {
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
      g->GetPoint( pIdx, qPoint );
      
      donorCellIdx =
        this->ProbeGridPointInAMR(qPoint,donorGrid, donorLevel, 
                                  amrds,maxLevelToLoad );
      
      if( donorCellIdx != -1 )
        {
        this->AverageLevel += donorLevel;
        assert(donorGrid != NULL);
        CD = donorGrid->GetCellData();
        this->CopyData( PD, pIdx, CD, donorCellIdx );
        }
      else
        {
        // Point is outside the domain, blank it
        ++ numPoints;
        g->BlankPoint( pIdx );
        }
      } // END for all grid nodes
    }
  std::cerr << "********* Resample Stats *************\n";
  double c = this->NumberOfSamples[0] * this->NumberOfSamples[1] * this->NumberOfSamples[2];
  double b = g->GetNumberOfPoints();
  std::cerr << "Number of Requested Points: " << c << " Number of Actual Points: " << b << "\n";
  std::cerr << " Percentage of Requested Points in Grid: " << 100.0 * b / c << "\n";
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
  std::cerr << "Total Number of Blocks Tested: " << this->NumberOfBlocksTested << "\n";
  std::cerr << " Number of Blocks that could be skipped by Visibility: " << this->NumberOfBlocksVisSkipped << "\n";
  double a = 100.0 * (double)(this->NumberOfBlocksVisSkipped) / (double) this->NumberOfBlocksTested;
  std::cerr << "Percentage of Blocks skipped via Visibility: " << a << "\n";
  a = (double) this->NumberOfBlocksTested / b;
  std::cerr << "Ave Number of Blocks Tested per Point: " << a << "\n";
  a = 100.0 * (double) this->NumberOfTimesFoundOnDonorLevel / b;
  std::cerr << "Percentage of Times we found point on Previous Level: " << a << "\n";
  a = 100.0 * (double) this->NumberOfTimesLevelUp / b;
  std::cerr << "Percentage of Times went to finer level: " << a << "\n";
  a = 100.0 * (double) this->NumberOfTimesLevelDown / b;
  std::cerr << "Percentage of Times went to coarser level: " << a << "\n";
  a = this->AverageLevel / b;
  std::cerr << "Average Level: " << a << "\n";
  std::cerr << "Number Of Failed Points: " << this->NumberOfFailedPoints << "\n";

George Zagaris's avatar
George Zagaris committed
839
840
841
}

//-----------------------------------------------------------------------------
842
void vtkAMRResampleFilter::TransferSolution(
843
    vtkUniformGrid *g, vtkOverlappingAMR *amrds)
George Zagaris's avatar
George Zagaris committed
844
845
846
847
848
{
  assert( "pre: uniform grid is NULL" && (g != NULL) );
  assert( "pre: AMR data-strucutre is NULL" && (amrds != NULL) );

  if( this->TransferToNodes == 1 )
849
    {
George Zagaris's avatar
George Zagaris committed
850
    this->TransferToGridNodes( g, amrds );
851
    }
George Zagaris's avatar
George Zagaris committed
852
  else
853
    {
George Zagaris's avatar
George Zagaris committed
854
    this->TransferToCellCenters( g, amrds );
855
    }
George Zagaris's avatar
George Zagaris committed
856
857
}

858
//-----------------------------------------------------------------------------
859
void vtkAMRResampleFilter::ExtractRegion(
860
    vtkOverlappingAMR *amrds, vtkMultiBlockDataSet *mbds,
861
    vtkOverlappingAMR * vtkNotUsed(metadata) )
862
{
863

864
  assert( "pre: input AMR data-structure is NULL" && (amrds != NULL) );
865
  assert( "pre: resampled grid should not be NULL" && (mbds != NULL) );
866

867
868
869
//  std::cout << "NumBlocks: " << this->ROI->GetNumberOfBlocks() << std::endl;
//  std::cout << "NumProcs: "  << this->Controller->GetNumberOfProcesses() << std::endl;
//  std::cout.flush();
870

Bob Obara's avatar
Bob Obara committed
871
872
  assert( "pre: NumProcs must be less than or equal to NumBlocks" &&
   ( static_cast<int>(this->ROI->GetNumberOfBlocks()) <= this->Controller->GetNumberOfProcesses()));
873

874
  mbds->SetNumberOfBlocks( this->ROI->GetNumberOfBlocks( ) );
875
  for( unsigned int block=0; block < this->ROI->GetNumberOfBlocks(); ++block )
876
    {
877
878
    if( this->IsRegionMine( block ) )
      {
879
880
881
882
883
884
885
886
887
      vtkUniformGrid *grid = vtkUniformGrid::New();
      grid->ShallowCopy( this->ROI->GetBlock( block ) );
      this->TransferSolution( grid, amrds );
      mbds->SetBlock( block, grid );
      grid->Delete();
      }
    else
      {
      mbds->SetBlock( block, NULL );
888
      }
889
    } // END for all blocks
890
891
892
893

}

//-----------------------------------------------------------------------------
894
void vtkAMRResampleFilter::ComputeAMRBlocksToLoad(
895
    vtkOverlappingAMR *metadata )
896
897
898
{
  assert( "pre: metadata is NULL" && (metadata != NULL) );

899
  this->BlocksToLoad.clear();
George Zagaris's avatar
George Zagaris committed
900

901
  unsigned int maxLevelToLoad = 0;
902
  if( this->LevelOfResolution < static_cast<int>(metadata->GetNumberOfLevels()))
903
    {
904
    maxLevelToLoad = this->LevelOfResolution+1;
905
    }
906
  else
907
    {
908
    maxLevelToLoad = metadata->GetNumberOfLevels();
909
    }
910
911
912
913

  unsigned int level=0;
  for( ;level < maxLevelToLoad; ++level )
    {
914
915
916
    unsigned int dataIdx = 0;
    for( ; dataIdx < metadata->GetNumberOfDataSets( level ); ++dataIdx )
      {
917
918
      vtkUniformGrid *grd = metadata->GetDataSet( level, dataIdx );
      assert( "pre: Metadata grid is NULL" && (grd != NULL) );
919

920
921
922
      if( this->IsBlockWithinBounds( grd ) )
        {
        this->BlocksToLoad.push_back(
923
        metadata->GetCompositeIndex(level,dataIdx) );
924
        } // END check if the block is within the bounds of the ROI
925
      } // END for all data
926
927
    } // END for all levels

928
   std::sort( this->BlocksToLoad.begin(), this->BlocksToLoad.end() );
Bob Obara's avatar
Bob Obara committed
929
   cerr << "Number Levels Loaded = " << maxLevelToLoad << " Number of Blocks = " << BlocksToLoad.size() << "\n";
930
931
932
}

//-----------------------------------------------------------------------------
George Zagaris's avatar
George Zagaris committed
933
void vtkAMRResampleFilter::GetDomainParameters(
934
    vtkOverlappingAMR *amr,
George Zagaris's avatar
George Zagaris committed
935
936
    double domainMin[3], double domainMax[3], double h[3],
    int dims[3], double &rf )
937
{
George Zagaris's avatar
George Zagaris committed
938
  assert( "pre: AMR dataset is NULL!" && (amr != NULL) );
George Zagaris's avatar
George Zagaris committed
939
940

  vtkAMRBox amrBox;
George Zagaris's avatar
George Zagaris committed
941
  amr->GetRootAMRBox( amrBox );
George Zagaris's avatar
George Zagaris committed
942

George Zagaris's avatar
George Zagaris committed
943
  rf = amr->GetRefinementRatio(1);
944

George Zagaris's avatar
George Zagaris committed
945
946
947
  amrBox.GetNumberOfNodes( dims );
  amrBox.GetMinBounds( domainMin );
  amrBox.GetMaxBounds( domainMax );
George Zagaris's avatar
George Zagaris committed
948
949
  amrBox.GetGridSpacing( h );
}
George Zagaris's avatar
George Zagaris committed
950

George Zagaris's avatar
George Zagaris committed
951
952
//-----------------------------------------------------------------------------
void vtkAMRResampleFilter::SnapBounds(
953
954
955
956
    const double* vtkNotUsed(h0[3]),
    const double domainMin[3],
    const double domainMax[3],
    const int* vtkNotUsed(dims[3]), bool outside[6] )
George Zagaris's avatar
George Zagaris committed
957
{
958
959
  int i, j;
  for(i=0, j=0; i < 3; ++i )
George Zagaris's avatar
George Zagaris committed
960
    {
961
962
    // Snap the parts of the bounds that lie outside of the AMR data
    if (this->Min[i] < domainMin[i])
George Zagaris's avatar
George Zagaris committed
963
      {
964
965
      outside[j++] = true;
      this->GridMin[i] = domainMin[i];
George Zagaris's avatar
George Zagaris committed
966
967
968
      }
    else
      {
969
      outside[j++] = false;
Bob Obara's avatar
Bob Obara committed
970
      this->GridMin[i] = this->Min[i];
George Zagaris's avatar
George Zagaris committed
971
      }
George Zagaris's avatar
George Zagaris committed
972

973
    if (this->Max[i] > domainMax[i])
George Zagaris's avatar
George Zagaris committed
974
      {
975
976
      outside[j++] = true;
      this->GridMax[i] = domainMax[i];
George Zagaris's avatar
George Zagaris committed
977
978
979
      }
    else
      {
980
981
      outside[j++] = false;
      this->GridMax[i] = this->Max[i];
George Zagaris's avatar
George Zagaris committed
982
983
      }
    }
George Zagaris's avatar
George Zagaris committed
984
}
George Zagaris's avatar
George Zagaris committed
985

George Zagaris's avatar
George Zagaris committed
986
987
988
989
990
//-----------------------------------------------------------------------------
void vtkAMRResampleFilter::ComputeLevelOfResolution(
    const int N[3], const double h0[3], const double L[3], const double rf )
{
  this->LevelOfResolution = 0;
991
992
  for( int i=0; i < 3; ++i )
    {
George Zagaris's avatar
George Zagaris committed
993
994
995
    double c1        = ( ( N[i]*h0[i] )/L[i] );
    int currentLevel = vtkMath::Floor( 0.5+(log(c1)/log(rf)) );
    if( currentLevel > this->LevelOfResolution )
996
      {
George Zagaris's avatar
George Zagaris committed
997
      this->LevelOfResolution = currentLevel;
998
      }
999
    } // END for all i
Bob Obara's avatar
Bob Obara committed
1000
  std::cerr << "Requested Max Level = " << this->LevelOfResolution << "\n";
George Zagaris's avatar
George Zagaris committed
1001
1002
}

George Zagaris's avatar
George Zagaris committed
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
//-----------------------------------------------------------------------------
bool vtkAMRResampleFilter::RegionIntersectsWithAMR(
    double domainMin[3], double domainMax[3],
    double regionMin[3], double regionMax[3])
{
  vtkBoundingBox domain;
  domain.SetMinPoint( domainMin );
  domain.SetMaxPoint( domainMax );

  vtkBoundingBox region;
  region.SetMinPoint( regionMin );
  region.SetMaxPoint( regionMax );

  if( domain.Intersects(region) )
1017
    {
George Zagaris's avatar
George Zagaris committed
1018
    return true;
1019
1020
    }

George Zagaris's avatar
George Zagaris committed
1021
1022
  return false;
}
1023

George Zagaris's avatar
George Zagaris committed
1024
//-----------------------------------------------------------------------------
George Zagaris's avatar
George Zagaris committed
1025
void vtkAMRResampleFilter::AdjustNumberOfSamplesInRegion(
1026
    const double Rh[3],
George Zagaris's avatar
George Zagaris committed
1027
    const bool outside[6], int N[3] )
George Zagaris's avatar
George Zagaris committed
1028
{
1029
1030
  int startIndex;
  int endIndex;
George Zagaris's avatar
George Zagaris committed
1031
1032
1033
1034
1035
1036
1037
1038
  double dx = 0.0;
  for( int i=0; i < 3; ++i )
    {
    N[i] = this->NumberOfSamples[i];

    // Get ijk of the snapped bounding wrt the requested virtual grid.
    if( outside[i*2] || outside[i*2+1] )
      {
1039
      dx = this->GridMin[i]-this->Min[i];
1040
1041
      if (dx > 0.0)
        {
1042
        startIndex = static_cast<int>( dx/Rh[i]+1 );
1043
1044
1045
        }
      else
        {
1046
        startIndex = 0;
1047
1048
        }

1049
1050
      dx  = this->GridMax[i]-this->Min[i];
      endIndex = static_cast<int>( dx/Rh[i]+1 );
George Zagaris's avatar
George Zagaris committed
1051

1052
      if (endIndex > N[i])
1053
        {
1054
        endIndex = N[i];
1055
        }
1056
1057
      int newN = endIndex - startIndex +1;
      if (newN <= N[i])
1058
        {
1059
        N[i] = newN;
1060
1061
1062
        }
      else
        {
George Zagaris's avatar
George Zagaris committed
1063
        assert("ERROR: code should not reach here!" && false );
1064
        }
George Zagaris's avatar
George Zagaris committed
1065
      }