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

 Program:   Visualization Toolkit
 Module:    vtkAMRCutPlane.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 "vtkAMRCutPlane.h"
#include "vtkObjectFactory.h"
17
18
19
20
21
22
23
24
25
26
27
28
29
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkMultiBlockDataSet.h"
#include "vtkCompositeDataPipeline.h"
#include "vtkMultiProcessController.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkIndent.h"
#include "vtkPlane.h"
#include "vtkAMRBox.h"
#include "vtkAMRUtilities.h"
#include "vtkUniformGrid.h"
#include "vtkCutter.h"
30
#include "vtkUnstructuredGrid.h"
31
32
33
34
35
36
37
#include "vtkIdList.h"
#include "vtkDoubleArray.h"
#include "vtkPointData.h"
#include "vtkCellData.h"
#include "vtkCell.h"
#include "vtkCellArray.h"
#include "vtkPoints.h"
38
#include "vtkOverlappingAMR.h"
39
40
#include "vtkPointData.h"
#include "vtkCellData.h"
41
42
43

#include <cassert>
#include <algorithm>
44

45
46
47
48
49
50

vtkStandardNewMacro(vtkAMRCutPlane);

//------------------------------------------------------------------------------
vtkAMRCutPlane::vtkAMRCutPlane()
{
51
52
  this->SetNumberOfInputPorts( 1 );
  this->SetNumberOfOutputPorts( 1 );
53
  this->LevelOfResolution = 0;
54
  this->initialRequest    = true;
55
56
  for( int i=0; i < 3; ++i )
    {
George Zagaris's avatar
George Zagaris committed
57
58
    this->Center[i] = 0.0;
    this->Normal[i] = 0.0;
59
    }
60
61
  this->Controller       = vtkMultiProcessController::GetGlobalController();
  this->UseNativeCutter  = 1;
62
63
64
65
66
}

//------------------------------------------------------------------------------
vtkAMRCutPlane::~vtkAMRCutPlane()
{
67
  this->blocksToLoad.clear();
68
69
70
71
72
73
}

//------------------------------------------------------------------------------
void vtkAMRCutPlane::PrintSelf( std::ostream &oss, vtkIndent indent )
{
  this->Superclass::PrintSelf( oss, indent );
74
  oss << indent << "LevelOfResolution: "
75
      << this->LevelOfResolution << endl;
76
  oss << indent << "UseNativeCutter: "
77
      << this->UseNativeCutter << endl;
78
  oss << indent << "Controller: "
79
      << this->Controller << endl;
80
81
82
83
84
  oss << indent << "Center: ";
  for( int i=0; i < 3; ++i )
    {
    oss << this->Center[i ] << " ";
    }
85
  oss << endl;
86
87
88
89
90
  oss << indent << "Normal: ";
  for( int i=0; i < 3; ++i )
    {
    oss << this->Normal[i] << " ";
    }
91
  oss << endl;
92
93
94
95
96
97
98
}

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

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

//------------------------------------------------------------------------------
int vtkAMRCutPlane::RequestInformation(
    vtkInformation* vtkNotUsed(rqst), vtkInformationVector** inputVector,
George Zagaris's avatar
George Zagaris committed
115
    vtkInformationVector* vtkNotUsed(outputVector) )
116
{
117
118
119
120
121
122
123
  this->blocksToLoad.clear();

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

  if( input->Has(vtkCompositeDataPipeline::COMPOSITE_DATA_META_DATA() ) )
    {
124
125
126
    vtkOverlappingAMR *metadata =
        vtkOverlappingAMR::SafeDownCast(
          input->Get(vtkCompositeDataPipeline::COMPOSITE_DATA_META_DATA()));
127

128
129
    vtkPlane *cutPlane = this->GetCutPlane( metadata );
    assert( "Cut plane is NULL" && (cutPlane != NULL) );
130

131
132
    this->ComputeAMRBlocksToLoad(cutPlane, metadata);
    cutPlane->Delete();
133
134
135
    }

  this->Modified();
136
137
138
139
140
  return 1;
}

//------------------------------------------------------------------------------
int vtkAMRCutPlane::RequestUpdateExtent(
141
    vtkInformation* vtkNotUsed(rqst), vtkInformationVector** inputVector,
George Zagaris's avatar
George Zagaris committed
142
    vtkInformationVector* vtkNotUsed(outputVector) )
143
{
144
145
146
147
148
149
  vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
  assert( "pre: inInfo is NULL" && (inInfo != NULL) );

  inInfo->Set(
      vtkCompositeDataPipeline::UPDATE_COMPOSITE_INDICES(),
      &this->blocksToLoad[0], this->blocksToLoad.size() );
150
151
152
153
154
  return 1;
}

//------------------------------------------------------------------------------
int vtkAMRCutPlane::RequestData(
155
156
    vtkInformation* vtkNotUsed(rqst), vtkInformationVector** inputVector,
    vtkInformationVector* outputVector )
157
{
158
159
160
  // STEP 0: Get input object
  vtkInformation *input = inputVector[0]->GetInformationObject( 0 );
  assert( "pre: input information object is NULL" && (input != NULL)  );
161
162
  vtkOverlappingAMR *inputAMR=
      vtkOverlappingAMR::SafeDownCast(
163
164
165
166
167
168
169
170
171
172
173
          input->Get( vtkDataObject::DATA_OBJECT() ) );
  assert( "pre: input AMR dataset is NULL!" && (inputAMR != NULL) );

  // STEP 1: Get output object
  vtkInformation *output = outputVector->GetInformationObject( 0 );
  assert( "pre: output information is NULL" && (output != NULL) );
  vtkMultiBlockDataSet *mbds=
      vtkMultiBlockDataSet::SafeDownCast(
          output->Get( vtkDataObject::DATA_OBJECT() ) );
  assert( "pre: output multi-block dataset is NULL" && (mbds != NULL) );

George Zagaris's avatar
George Zagaris committed
174
175
  if( this->IsAMRData2D( inputAMR ) )
    {
176
177
    // Return an empty multi-block, we cannot cut a 2-D dataset
    return 1;
George Zagaris's avatar
George Zagaris committed
178
    }
179

180
181
  vtkPlane *cutPlane = this->GetCutPlane( inputAMR );
  assert("pre: cutPlane should not be NULL!" && (cutPlane != NULL) );
182

183
184
185
186
  unsigned int blockIdx = 0;
  unsigned int level    = 0;
  for( ; level < inputAMR->GetNumberOfLevels(); ++level )
    {
187
188
189
190
191
    unsigned int dataIdx = 0;
    for( ; dataIdx < inputAMR->GetNumberOfDataSets( level ); ++dataIdx )
      {
      vtkUniformGrid *grid = inputAMR->GetDataSet( level, dataIdx );
      if( this->UseNativeCutter == 1 )
192
        {
193
194
195
196
        if( grid != NULL )
          {
          vtkCutter *myCutter = vtkCutter::New();
          myCutter->SetInputData( grid );
197
          myCutter->SetCutFunction( cutPlane );
198
199
200
201
202
203
204
205
206
207
208
209
210
          myCutter->Update();
          mbds->SetBlock( blockIdx, myCutter->GetOutput( ) );
          ++blockIdx;
          myCutter->Delete();
          }
        else
          {
          mbds->SetBlock(blockIdx,NULL);
          ++blockIdx;
          }
        }
      else
        {
211
212
        if( grid != NULL )
          {
213
          this->CutAMRBlock( cutPlane, blockIdx, grid, mbds );
214
          ++blockIdx;
215
216
217
218
219
220
          }
        else
          {
          mbds->SetBlock(blockIdx,NULL);
          ++blockIdx;
          }
221
222
        }
      } // END for all data
223
224
    } // END for all levels

225
  cutPlane->Delete();
226
227
  return 1;
}
228

229
230
//------------------------------------------------------------------------------
void vtkAMRCutPlane::CutAMRBlock(
231
    vtkPlane *cutPlane,
232
    unsigned int blockIdx, vtkUniformGrid *grid, vtkMultiBlockDataSet *output )
233
{
234
235
  assert("pre: multiblock output object is NULL!" && (output != NULL));
  assert("pre: grid is NULL" && (grid != NULL) );
236

237
  vtkUnstructuredGrid *mesh       = vtkUnstructuredGrid::New();
238
  vtkPoints *meshPts      = vtkPoints::New();
239
  meshPts->SetDataTypeToDouble();
240
  vtkCellArray *cells     = vtkCellArray::New();
241

242
243
244
  // Maps points from the input grid to the output grid
  std::map< vtkIdType, vtkIdType > grdPntMapping;
  std::vector< vtkIdType > extractedCells;
245
246
247
248

  vtkIdType cellIdx = 0;
  for( ; cellIdx < grid->GetNumberOfCells(); ++cellIdx )
    {
249
    if( grid->IsCellVisible( cellIdx ) &&
250
        this->PlaneIntersectsCell( cutPlane, grid->GetCell(cellIdx) ) )
251
      {
252
      extractedCells.push_back( cellIdx );
253
      this->ExtractCellFromGrid(
254
          grid,grid->GetCell(cellIdx),grdPntMapping,meshPts,cells );
255
      } // END if
256
257
    } // END for all cells

258
259
260
261
262
263
  // Sanity checks
  assert("post: Number of mesh points should match map size!" &&
       (static_cast<int>(grdPntMapping.size())==meshPts->GetNumberOfPoints()));
  assert("post: Number of cells mismatch"&&
       (cells->GetNumberOfCells()==static_cast<int>(extractedCells.size())));

264
265
266
267
  // Insert the points
  mesh->SetPoints( meshPts );
  meshPts->Delete();

268
269
270
271
272
273
  std::vector<int> types;
  if( grid->GetDataDimension() == 3 )
    {
    types.resize( cells->GetNumberOfCells(), VTK_VOXEL );
    }
  else
274
    {
275
276
277
278
    vtkErrorMacro("Cannot cut a grid of dimension=" << grid->GetDataDimension());
    output->SetBlock( blockIdx, NULL );
    return;
    }
279
280

  // Insert the cells
281
  mesh->SetCells( &types[0], cells );
282
283
  cells->Delete();

284
285
286
287
288
289
  // Extract fields
  this->ExtractPointDataFromGrid(
      grid,grdPntMapping,mesh->GetNumberOfPoints(),mesh->GetPointData());
  this->ExtractCellDataFromGrid(
      grid,extractedCells,mesh->GetCellData());

290
  output->SetBlock( blockIdx, mesh );
291
  mesh->Delete();
292
293
  grdPntMapping.clear();
  extractedCells.clear();
294
295
296
297
298
}

//------------------------------------------------------------------------------
void vtkAMRCutPlane::ExtractCellFromGrid(
            vtkUniformGrid *grid,
299
300
            vtkCell* cell, std::map< vtkIdType, vtkIdType >& grdPntMapping,
            vtkPoints *nodes,
301
            vtkCellArray *cells)
302
303
304
305
306
{
  assert( "pre: grid is NULL"  && (grid != NULL)  );
  assert( "pre: cell is NULL"  && (cell != NULL)  );
  assert( "pre: cells is NULL" && (cells != NULL) );

307
  cells->InsertNextCell( cell->GetNumberOfPoints() );
308
309
310
  vtkIdType nodeIdx = 0;
  for( ; nodeIdx < cell->GetNumberOfPoints(); ++nodeIdx )
    {
311
    // Get the point ID w.r.t. the grid
312
313
314
    vtkIdType meshPntIdx = cell->GetPointId( nodeIdx );
    assert( "pre: mesh point ID should within grid range point ID" &&
            (meshPntIdx < grid->GetNumberOfPoints()));
315

316
317
318
319
320
321
322
323
324
    if( grdPntMapping.find( meshPntIdx ) != grdPntMapping.end() )
      {
      // Point already exists in nodes
      cells->InsertCellPoint( grdPntMapping[meshPntIdx] );
      }
    else
      {
      // Push point to the end of the list
      vtkIdType nidx = nodes->GetNumberOfPoints();
325
326
      double *pnt    = grid->GetPoint(meshPntIdx);
      nodes->InsertPoint( nidx, pnt );
327
328
329
      assert("post: number of points should be increased by 1" &&
             (nodes->GetNumberOfPoints()==(nidx+1)));
      grdPntMapping[ meshPntIdx ] = nidx;
330
      cells->InsertCellPoint( nidx );
331
      }
332
333
334
335
    } // END for all nodes

}

336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
//------------------------------------------------------------------------------
void vtkAMRCutPlane::ExtractPointDataFromGrid(
    vtkUniformGrid *grid,
    std::map<vtkIdType,vtkIdType>& gridPntMapping,
    vtkIdType NumNodes,
    vtkPointData *PD)
{
  assert("pre: grid is NULL!" && (grid != NULL) );
  assert("pre: target point data is NULL!" && (PD != NULL) );

  if( (grid->GetPointData()->GetNumberOfArrays()==0) ||
      (gridPntMapping.size() == 0))
    {
    // Nothing to extract short-circuit here
    return;
    }

  vtkPointData *GPD = grid->GetPointData();
  for( int fieldArray=0; fieldArray < GPD->GetNumberOfArrays(); ++fieldArray)
    {
    vtkDataArray *sourceArray = GPD->GetArray(fieldArray);
    int dataType = sourceArray->GetDataType();
    vtkDataArray *array = vtkDataArray::CreateDataArray( dataType );
    assert( "pre: failed to create array!" && (array != NULL) );

    array->SetName(sourceArray->GetName());
    array->SetNumberOfComponents(sourceArray->GetNumberOfComponents());
    array->SetNumberOfTuples(NumNodes);

    // Copy tuples from source array
    std::map<vtkIdType,vtkIdType>::iterator iter = gridPntMapping.begin();
    for( ; iter != gridPntMapping.end(); ++iter )
      {
      vtkIdType srcIdx    = iter->first;
      vtkIdType targetIdx = iter->second;
      assert( "pre: source node index is out-of-bounds" &&
              (srcIdx >= 0) && (srcIdx < grid->GetNumberOfPoints()) );
      assert( "pre: target node index is out-of-bounds" &&
              (targetIdx >=0) && (targetIdx < NumNodes)  );
      array->SetTuple( targetIdx, srcIdx, sourceArray );
      } // END for all extracted nodes

    PD->AddArray( array );
    array->Delete();
    } // END for all arrays
}

//------------------------------------------------------------------------------
void vtkAMRCutPlane::ExtractCellDataFromGrid(
    vtkUniformGrid *grid,
    std::vector<vtkIdType>& cellIdxList,
    vtkCellData *CD)
{
  assert("pre: grid is NULL!" && (grid != NULL) );
  assert("pre: target cell data is NULL!" && (CD != NULL) );

  if( (grid->GetCellData()->GetNumberOfArrays()==0) ||
      (cellIdxList.size()==0) )
    {
    // Nothing to extract short-circuit here
    return;
    }

  int NumCells = static_cast<int>(cellIdxList.size());
  vtkCellData *GCD = grid->GetCellData();
  for( int fieldArray=0; fieldArray < GCD->GetNumberOfArrays(); ++fieldArray)
  {
  vtkDataArray *sourceArray = GCD->GetArray(fieldArray);
  int dataType = sourceArray->GetDataType();
  vtkDataArray *array = vtkDataArray::CreateDataArray(dataType);
  assert( "pre: failed to create array!" && (array != NULL) );

  array->SetName(sourceArray->GetName());
  array->SetNumberOfComponents(sourceArray->GetNumberOfComponents());
  array->SetNumberOfTuples(NumCells);

  // Copy tuples from source array
  for( int i=0; i < NumCells; ++i )
    {
    vtkIdType cellIdx = cellIdxList[ i ];
    assert( "pre: cell index is out-of-bounds!" &&
          (cellIdx >= 0) && (cellIdx < grid->GetNumberOfCells()) );
    array->SetTuple(i,cellIdx,sourceArray);
    } // END for all extracted cells

  CD->AddArray( array );
  array->Delete();
  } // END for all arrays

}

427
//------------------------------------------------------------------------------
428
vtkPlane* vtkAMRCutPlane::GetCutPlane( vtkOverlappingAMR *metadata )
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
{
  assert( "pre: metadata is NULL" && (metadata != NULL) );

  vtkPlane *pl = vtkPlane::New();

  // Get global bounds
  double minBounds[3];
  double maxBounds[3];
  vtkAMRBox root;
  metadata->GetRootAMRBox( root );
  root.GetMinBounds( minBounds );
  root.GetMaxBounds( maxBounds );

  this->InitializeCenter( minBounds, maxBounds );

  pl->SetNormal( this->Normal );
  pl->SetOrigin( this->Center );
  return( pl );
}

//------------------------------------------------------------------------------
void vtkAMRCutPlane::ComputeAMRBlocksToLoad(
451
      vtkPlane* p, vtkOverlappingAMR *m)
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
{
  assert( "pre: Plane object is NULL" && (p != NULL) );
  assert( "pre: metadata is NULL" && (m != NULL) );

  // Store A,B,C,D from the plane equation
  double plane[4];
  plane[0] = p->GetNormal()[0];
  plane[1] = p->GetNormal()[1];
  plane[2] = p->GetNormal()[2];
  plane[3] = p->GetNormal()[0]*p->GetOrigin()[0] +
             p->GetNormal()[1]*p->GetOrigin()[1] +
             p->GetNormal()[2]*p->GetOrigin()[2];

  double bounds[6];

George Zagaris's avatar
George Zagaris committed
467
  int NumLevels = m->GetNumberOfLevels();
468
  int maxLevelToLoad =
George Zagaris's avatar
George Zagaris committed
469
470
     (this->LevelOfResolution < NumLevels )?
         this->LevelOfResolution : NumLevels;
471
472

  unsigned int level = 0;
George Zagaris's avatar
George Zagaris committed
473
  for( ; level <= static_cast<unsigned int>(maxLevelToLoad); ++level )
474
    {
475
476
477
478
479
480
481
482
483
484
485
486
487
    unsigned int dataIdx = 0;
    for( ; dataIdx < m->GetNumberOfDataSets( level ); ++dataIdx )
      {
      vtkAMRBox box;
      m->GetMetaData( level, dataIdx, box  );
      bounds[0] = box.GetMinX();
      bounds[1] = box.GetMaxX();
      bounds[2] = box.GetMinY();
      bounds[3] = box.GetMaxY();
      bounds[4] = box.GetMinZ();
      bounds[5] = box.GetMaxZ();

      if( this->PlaneIntersectsAMRBox( plane, bounds ) )
488
        {
489
490
491
492
        unsigned int amrGridIdx = m->GetCompositeIndex(level,dataIdx);
        this->blocksToLoad.push_back( amrGridIdx );
        }
      } // END for all data
493
494
    } // END for all levels

495
  std::sort( this->blocksToLoad.begin(), this->blocksToLoad.end() );
496
497
498
499
500
501
}

//------------------------------------------------------------------------------
void vtkAMRCutPlane::InitializeCenter( double min[3], double max[3] )
{
  if( !this->initialRequest )
502
    {
503
    return;
504
    }
505
506
507
508
509
510
511

  this->Center[0] = 0.5*( max[0]-min[0] );
  this->Center[1] = 0.5*( max[1]-min[1] );
  this->Center[2] = 0.5*( max[2]-min[2] );
  this->initialRequest = false;
}

512
//------------------------------------------------------------------------------
513
bool vtkAMRCutPlane::PlaneIntersectsCell( vtkPlane *pl, vtkCell *cell )
514
{
515
  assert( "pre: plane is NULL" && (pl != NULL) );
516
  assert( "pre: cell is NULL!" && (cell != NULL) );
517
  return( this->PlaneIntersectsAMRBox( pl, cell->GetBounds() ) );
518
519
}
//------------------------------------------------------------------------------
520
bool vtkAMRCutPlane::PlaneIntersectsAMRBox(vtkPlane *pl, double bounds[6] )
521
{
522
  assert( "pre: plane is NULL" && (pl != NULL) );
523

524
525
  // Store A,B,C,D from the plane equation
  double plane[4];
526
527
528
529
530
531
  plane[0] = pl->GetNormal()[0];
  plane[1] = pl->GetNormal()[1];
  plane[2] = pl->GetNormal()[2];
  plane[3] = pl->GetNormal()[0]*pl->GetOrigin()[0] +
             pl->GetNormal()[1]*pl->GetOrigin()[1] +
             pl->GetNormal()[2]*pl->GetOrigin()[2];
532
533
534
535

 return( this->PlaneIntersectsAMRBox( plane,bounds) );
}

536
537
538
539
540
541
542
543
//------------------------------------------------------------------------------
bool vtkAMRCutPlane::PlaneIntersectsAMRBox( double plane[4], double bounds[6] )
{
  bool lowPnt  = false;
  bool highPnt = false;

  for( int i=0; i < 8; ++i )
    {
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
    // Get box coordinates
    double x = ( i&1 ? bounds[1] : bounds[0] );
    double y = ( i&2 ? bounds[3] : bounds[2] );
    double z = ( i&3 ? bounds[5] : bounds[4] );

    // Plug-in coordinates to the plane equation
    double v = plane[3] - plane[0]*x - plane[1]*y - plane[2]*z;

    if( v == 0.0 ) // Point is on a plane
      {
      return true;
      }

    if( v < 0.0 )
      {
      lowPnt = true;
      }
    else
      {
      highPnt = true;
      }

    if( lowPnt && highPnt )
      {
      return true;
      }
570
571
    }

572
  return false;
573
574
575
}

//------------------------------------------------------------------------------
576
bool vtkAMRCutPlane::IsAMRData2D( vtkOverlappingAMR *input )
577
578
579
580
581
582
583
{
  assert( "pre: Input AMR dataset is NULL" && (input != NULL)  );

  vtkAMRBox box;
  input->GetMetaData( 0, 0, box );

  if( box.GetDimensionality() == 2 )
584
585
586
    {
    return true;
    }
587
588
589

 return false;
}