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

  Program:   Visualization Toolkit
  Module:    vtkPExtractRectilinearGrid.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 "vtkPExtractRectilinearGrid.h"

// VTK includes
18
#include "vtkStructuredExtent.h"
19
#include "vtkExtractStructuredGridHelper.h"
20 21 22 23 24 25 26 27
#include "vtkRectilinearGrid.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkMPIController.h"
#include "vtkMPIUtilities.h"
#include "vtkMultiProcessController.h"
#include "vtkObjectFactory.h"
#include "vtkStreamingDemandDrivenPipeline.h"
28
#include "vtkStructuredExtent.h"
29 30 31 32 33
#include "vtkStructuredImplicitConnectivity.h"

#include <cassert>
#include <sstream>

34 35 36 37
// Some useful extent macros
#define EMIN(ext, dim) (ext[2*dim])
#define EMAX(ext, dim) (ext[2*dim+1])

38 39
// #define DEBUG

40 41 42
#ifdef DEBUG
#define DEBUG_EXTENT(label, extent) \
  if (this->Controller) \
43
  { \
44 45 46
    vtkMPIUtilities::SynchronizedPrintf( \
        this->Controller, #label "=[%d,%d,%d,%d,%d,%d]\n", \
        extent[0], extent[1], extent[2], extent[3], extent[4], extent[5]); \
47
  } \
48
  else \
49
  { \
50 51 52 53
    std::cout << label << "=[" \
              << extent[0] << "," << extent[1] << "," \
              << extent[2] << "," << extent[3] << "," \
              << extent[4] << "," << extent[5] << "]\n"; \
54
  }
55 56 57

#define DEBUG_OUT(out) \
  if (this->Controller) \
58
  { \
59 60 61 62
    std::ostringstream tmpStreamOut; \
    tmpStreamOut << out; \
    vtkMPIUtilities::SynchronizedPrintf(this->Controller, \
                                        tmpStreamOut.str().c_str()); \
63
  } \
64
  else \
65
  { \
66
    std::cout << out; \
67
  }
68 69 70 71 72
#else // DEBUG
#define DEBUG_EXTENT(label, extent)
#define DEBUG_OUT(out)
#endif // DEBUG

73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95
vtkStandardNewMacro(vtkPExtractRectilinearGrid);

//------------------------------------------------------------------------------
vtkPExtractRectilinearGrid::vtkPExtractRectilinearGrid()
{
  this->Controller = vtkMPIController::SafeDownCast(
           vtkMultiProcessController::GetGlobalController());
}

//------------------------------------------------------------------------------
vtkPExtractRectilinearGrid::~vtkPExtractRectilinearGrid()
{

}

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

//------------------------------------------------------------------------------
int vtkPExtractRectilinearGrid::RequestData(
96 97 98
        vtkInformation *request,
        vtkInformationVector **inputVector,
        vtkInformationVector *outputVector)
99
{
100 101
  DEBUG_OUT("########### RequestData\n");

102 103 104 105 106 107
  bool isSubSampling = this->SampleRate[0] != 1 ||
                       this->SampleRate[1] != 1 ||
                       this->SampleRate[2] != 1;

  // No MPI, or no subsampling? Just run the serial implementation.
  if (!this->Controller || !isSubSampling)
108
  {
109
    return this->Superclass::RequestData(request, inputVector, outputVector);
110
  }
111

112
    if (!this->Internal->IsValid())
113
    {
114
      return 0;
115
    }
116

117 118 119
  // Collect information:
  vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
  vtkInformation *outInfo = outputVector->GetInformationObject(0);
120

121 122 123 124 125 126
  int inputWholeExtent[6];
  inInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),
               inputWholeExtent);
  int outputWholeExtent[6];
  outInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),
               outputWholeExtent);
127

128 129 130 131
  vtkRectilinearGrid *input = vtkRectilinearGrid::SafeDownCast(
          inInfo->Get(vtkDataObject::DATA_OBJECT()));
  vtkRectilinearGrid *output = vtkRectilinearGrid::SafeDownCast(
        outInfo->Get(vtkDataObject::DATA_OBJECT()));
132

133 134
  int inputExtent[6];
  input->GetExtent(inputExtent);
135

136 137 138 139 140
  // Clamp the global VOI to the whole extent:
  int globalVOI[6];
  std::copy(this->VOI, this->VOI + 6, globalVOI);
  vtkStructuredExtent::Clamp(globalVOI, inputWholeExtent);

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
  // 1D Example:
  //   InputWholeExtent = [0, 20]
  //   GlobalVOI = [3, 17]
  //   SampleRate = 2
  //   OutputWholeExtent = [0, 7]
  //   Processes = 2
  //
  // Process 0:
  //   PartitionedInputExtent = [0, 10]
  //   PartitionedVOI = [3, 9] (due to sampling)
  //   OutputExtent = [0, 3]
  //   SerialOutputExtent = [0, 3]
  //   FinalOutputExtent = [0, 4] (after gap closing)
  //
  // Process 1:
  //   PartitionedInputExtent = [10, 20]
  //   PartitionedVOI = [11, 17] (offset due to sampling)
  //   OutputExtent = [4, 7]
  //   SerialOutputExtent = [0, 3]
  //   FinalOutputExtent = [4, 7]
  //
  // This filter should:
  // 1) Compute ParititonedVOI that will allow the base class to produce as much
  //    of the output data set as possible from the partitioned piece.
  //
166
  // 2) Update the output dataset's extents to match PartitionedOutputExtent (it
167 168
  //    will be [0, L] in each dimension by default).
  //
169 170
  // 3) Extract PartitionedVOI using the base class's implementation.
  //
171 172
  // 4) Close gaps using vtkStructuredImplicitConnectivity (e.g. [3, 4] in the
  //    above example).
173

174 175
  bool partitionContainsVOI = true;
  for (int dim = 0; partitionContainsVOI && dim < 3; ++dim)
176
  {
177 178
    partitionContainsVOI = EMAX(inputExtent, dim) >= EMIN(globalVOI, dim) &&
                           EMIN(inputExtent, dim) <= EMAX(globalVOI, dim);
179
  }
180 181

  DEBUG_EXTENT("InputWholeExtent", inputWholeExtent);
182
  DEBUG_EXTENT("OutputWholeExtent", outputWholeExtent);
183
  DEBUG_EXTENT("GlobalVOI", globalVOI);
184
  DEBUG_EXTENT("InputPartitionedExtent", inputExtent);
185

186 187 188 189
  int partitionedVOI[6] = {0, -1, 0, -1, 0, -1};
  int partitionedOutputExtent[6] = {0, -1, 0, -1, 0, -1};

  if (partitionContainsVOI)
190
  {
191 192 193
    ////////////////////////////////////////////////////////////////
    // 1) Compute actual VOI for aligning the partitions outputs: //
    ////////////////////////////////////////////////////////////////
194 195 196
    vtkExtractStructuredGridHelper::GetPartitionedVOI(
          globalVOI, inputExtent, this->SampleRate, this->IncludeBoundary != 0,
          partitionedVOI);
197
  }
198
  DEBUG_EXTENT("PartitionedVOI", partitionedVOI);
199

200
  if (partitionContainsVOI)
201
  {
202
    ////////////////////////////////////////////////////////////////
203
    // 2) Compute and update the output dataset's actual extents. //
204 205
    ////////////////////////////////////////////////////////////////
    vtkExtractStructuredGridHelper::GetPartitionedOutputExtent(
206
          globalVOI, partitionedVOI, outputWholeExtent, this->SampleRate,
207
          this->IncludeBoundary != 0, partitionedOutputExtent);
208
    output->SetExtent(partitionedOutputExtent);
209
  }
210
  DEBUG_EXTENT("PartitionedOutputExtent", partitionedOutputExtent);
211 212

  if (partitionContainsVOI)
213
  {
214 215 216 217
    ////////////////////////////////////////////////////////////
    // 3) Extract actual VOI using superclass implementation: //
    ////////////////////////////////////////////////////////////
    if (!this->Superclass::RequestDataImpl(inputVector, outputVector))
218
    {
219 220
      return 0;
    }
221
  }
222 223 224 225

  //////////////////////////////
  // 4: Detect & resolve gaps //
  //////////////////////////////
226 227
  vtkStructuredImplicitConnectivity* gridConnectivity =
      vtkStructuredImplicitConnectivity::New();
228
  gridConnectivity->SetWholeExtent(outputWholeExtent);
229 230 231 232

  // Register the grid, grid ID is the same as the process ID
  gridConnectivity->RegisterRectilinearGrid(
    this->Controller->GetLocalProcessId(),
233 234 235 236 237
    output->GetExtent(),
    output->GetXCoordinates(),
    output->GetYCoordinates(),
    output->GetZCoordinates(),
    output->GetPointData()
238 239 240 241 242 243 244
    );

  // Establish neighbor connectivity & detect any gaps
  gridConnectivity->EstablishConnectivity();

  // Check if there are any gaps, if any close them now
  if( gridConnectivity->HasImplicitConnectivity() )
245
  {
246
    DEBUG_OUT("Closing gaps...\n");
247 248
    // there are gaps, grow the grid to the right
    gridConnectivity->ExchangeData();
249

250
    gridConnectivity->GetOutputRectilinearGrid(
251
      this->Controller->GetLocalProcessId(),output);
252
  }
253

254 255
  gridConnectivity->Delete();

256
#ifdef DEBUG
257 258 259
  int finalOutputExtent[6];
  output->GetExtent(finalOutputExtent);
  DEBUG_EXTENT("FinalOutputExtent", finalOutputExtent);
260 261
#endif

262
  return 1;
263 264 265 266
}

//------------------------------------------------------------------------------
int vtkPExtractRectilinearGrid::RequestInformation(
267 268 269
        vtkInformation *request,
        vtkInformationVector **inputVector,
        vtkInformationVector *outputVector)
270
{
271
  DEBUG_OUT("########### RequestInformation\n");
272

273 274
  return this->Superclass::RequestInformation(request, inputVector,
                                              outputVector);
275 276 277 278 279 280 281 282
}

//------------------------------------------------------------------------------
int vtkPExtractRectilinearGrid::RequestUpdateExtent(
        vtkInformation* request,
        vtkInformationVector** inputVector,
        vtkInformationVector* outputVector)
{
283 284 285 286
  DEBUG_OUT("########### RequestUpdateExtent\n");

  return this->Superclass::RequestUpdateExtent(request, inputVector,
                                               outputVector);
287
}