vtkGradientFilter.cxx 48.7 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28
/*=========================================================================

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

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

/*----------------------------------------------------------------------------
 Copyright (c) Sandia Corporation
 See Copyright.txt or http://www.paraview.org/HTML/Copyright.html for details.
----------------------------------------------------------------------------*/

#include "vtkGradientFilter.h"

#include "vtkCell.h"
#include "vtkCellData.h"
#include "vtkCellDataToPointData.h"
#include "vtkDataArray.h"
#include "vtkDataSet.h"
#include "vtkIdList.h"
29
#include "vtkImageData.h"
30 31 32
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkMath.h"
33
#include "vtkNew.h"
34 35 36
#include "vtkObjectFactory.h"
#include "vtkPointData.h"
#include "vtkPolyData.h"
37
#include "vtkRectilinearGrid.h"
38
#include "vtkSmartPointer.h"
39
#include "vtkStreamingDemandDrivenPipeline.h"
40
#include "vtkStructuredGrid.h"
41 42
#include "vtkUnstructuredGrid.h"

43
#include <limits>
44
#include <vector>
45

46 47 48 49
//-----------------------------------------------------------------------------

vtkStandardNewMacro(vtkGradientFilter);

50
namespace
51
{
52 53 54 55 56 57 58 59
// special template macro for only float and double types since we will never
// have other data types for output arrays and this helps with reducing template expansion
#define vtkFloatingPointTemplateMacro(call)         \
  vtkTemplateMacroCase(VTK_DOUBLE, double, call);   \
  vtkTemplateMacroCase(VTK_FLOAT, float, call);

// helper function to replace the gradient of a vector
// with the vorticity/curl of that vector
60 61
//-----------------------------------------------------------------------------
  template<class data_type>
62
  void ComputeVorticityFromGradient(data_type* gradients, data_type* vorticity)
63
  {
64 65 66
    vorticity[0] = gradients[7] - gradients[5];
    vorticity[1] = gradients[2] - gradients[6];
    vorticity[2] = gradients[3] - gradients[1];
67 68
  }

69 70 71 72 73 74
  template<class data_type>
  void ComputeDivergenceFromGradient(data_type* gradients, data_type* divergence)
  {
    divergence[0] = gradients[0]+gradients[4]+gradients[8];
  }

75 76 77
  template<class data_type>
  void ComputeQCriterionFromGradient(data_type* gradients, data_type* qCriterion)
  {
78 79 80 81 82
    // see http://public.kitware.com/pipermail/paraview/2015-May/034233.html for
    // paper citation and formula on Q-criterion.
    qCriterion[0] =
      - (gradients[0]*gradients[0]+gradients[4]*gradients[4]+gradients[8]*gradients[8])/2.
      - (gradients[1]*gradients[3]+gradients[2]*gradients[6]+gradients[5]*gradients[7]);
83 84
  }

85 86 87
  // Functions for unstructured grids and polydatas
  template<class data_type>
  void ComputePointGradientsUG(
88
    vtkDataSet *structure, vtkDataArray *array, data_type *gradients,
89
    int numberOfInputComponents, data_type* vorticity, data_type* qCriterion,
90
    data_type* divergence, int highestCellDimension, int contributingCellOption);
91 92

  int GetCellParametricData(
93
    vtkIdType pointId, double pointCoord[3], vtkCell *cell, int & subId,
94
    double parametricCoord[3]);
95

96 97
  template<class data_type>
  void ComputeCellGradientsUG(
98
    vtkDataSet *structure, vtkDataArray *array, data_type *gradients,
99 100
    int numberOfInputComponents, data_type* vorticity, data_type* qCriterion,
    data_type* divergence);
101 102 103

  // Functions for image data and structured grids
  template<class Grid, class data_type>
104
  void ComputeGradientsSG(Grid output, vtkDataArray* array, data_type* gradients,
105
                          int numberOfInputComponents, int fieldAssociation,
106 107
                          data_type* vorticity, data_type* qCriterion,
                          data_type* divergence);
108

109 110
  bool vtkGradientFilterHasArray(vtkFieldData *fieldData,
                                 vtkDataArray *array)
111 112 113
  {
    int numarrays = fieldData->GetNumberOfArrays();
    for (int i = 0; i < numarrays; i++)
114
    {
115
      if (array == fieldData->GetArray(i))
116
      {
117
        return true;
118
      }
119
    }
120
    return false;
121 122
  }

123
  // generic way to get the coordinate for either a cell (using
124
  // the parametric center) or a point
125
  void GetGridEntityCoordinate(vtkDataSet* grid, int fieldAssociation,
126
                               vtkIdType index, double coords[3])
127 128
  {
    if(fieldAssociation == vtkDataObject::FIELD_ASSOCIATION_POINTS)
129
    {
130
      grid->GetPoint(index, coords);
131
    }
132
    else
133
    {
134
      vtkCell* cell = grid->GetCell(index);
135
      double pcoords[3];
136
      int subId = cell->GetParametricCenter(pcoords);
137
      std::vector<double> weights(cell->GetNumberOfPoints()+1);
138
      cell->EvaluateLocation(subId, pcoords, coords, &weights[0]);
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

  template<class data_type>
  void Fill(vtkDataArray* array, data_type vtkNotUsed(data), int replacementValueOption)
  {
    switch (replacementValueOption)
    {
    case vtkGradientFilter::Zero:
      array->Fill(0);
      return;
    case vtkGradientFilter::NaN:
      array->Fill(vtkMath::Nan());
      return;
    case vtkGradientFilter::DataTypeMin:
      array->Fill(std::numeric_limits<data_type>::min());
      return;
    case vtkGradientFilter::DataTypeMax:
      array->Fill(std::numeric_limits<data_type>::max());
      return;
    }
  }
  template<class data_type>
  int GetOutputDataType(data_type vtkNotUsed(data))
  {
    if(sizeof(data_type)>=8)
    {
      return VTK_DOUBLE;
    }
    return VTK_FLOAT;
  }
170
} // end anonymous namespace
171 172 173 174

//-----------------------------------------------------------------------------
vtkGradientFilter::vtkGradientFilter()
{
175 176 177 178
  this->ResultArrayName = nullptr;
  this->DivergenceArrayName = nullptr;
  this->VorticityArrayName = nullptr;
  this->QCriterionArrayName = nullptr;
179
  this->FasterApproximation = 0;
180 181
  this->ComputeGradient = 1;
  this->ComputeDivergence = 0;
182
  this->ComputeVorticity = 0;
183
  this->ComputeQCriterion = 0;
184 185
  this->ContributingCellOption = vtkGradientFilter::All;
  this->ReplacementValueOption = vtkGradientFilter::Zero;
186 187 188 189
  this->SetInputScalars(vtkDataObject::FIELD_ASSOCIATION_POINTS_THEN_CELLS,
                        vtkDataSetAttributes::SCALARS);
}

190
//-----------------------------------------------------------------------------
191 192
vtkGradientFilter::~vtkGradientFilter()
{
193 194 195 196
  this->SetResultArrayName(nullptr);
  this->SetDivergenceArrayName(nullptr);
  this->SetVorticityArrayName(nullptr);
  this->SetQCriterionArrayName(nullptr);
197 198
}

199
//-----------------------------------------------------------------------------
200 201 202 203
void vtkGradientFilter::PrintSelf(ostream &os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os, indent);

204
  os << indent << "ResultArrayName:"
205
     << (this->ResultArrayName ? this->ResultArrayName : "Gradients") << endl;
206 207
  os << indent << "DivergenceArrayName:"
     << (this->DivergenceArrayName ? this->DivergenceArrayName : "Divergence") << endl;
208 209 210 211
  os << indent << "VorticityArrayName:"
     << (this->VorticityArrayName ? this->VorticityArrayName : "Vorticity") << endl;
  os << indent << "QCriterionArrayName:"
     << (this->QCriterionArrayName ? this->QCriterionArrayName : "Q-criterion") << endl;
212
  os << indent << "FasterApproximation:" << this->FasterApproximation << endl;
213 214
  os << indent << "ComputeGradient:"  << this->ComputeGradient << endl;
  os << indent << "ComputeDivergence:"  << this->ComputeDivergence << endl;
215
  os << indent << "ComputeVorticity:" << this->ComputeVorticity << endl;
216
  os << indent << "ComputeQCriterion:" << this->ComputeQCriterion << endl;
217 218
  os << indent << "ContributingCellOption:" << this->ContributingCellOption << endl;
  os << indent << "ReplacementValueOption:" << this->ReplacementValueOption << endl;
219 220 221 222 223 224 225 226
}

//-----------------------------------------------------------------------------
void vtkGradientFilter::SetInputScalars(int fieldAssociation, const char *name)
{
  if (   (fieldAssociation != vtkDataObject::FIELD_ASSOCIATION_POINTS)
      && (fieldAssociation != vtkDataObject::FIELD_ASSOCIATION_CELLS)
      && (fieldAssociation != vtkDataObject::FIELD_ASSOCIATION_POINTS_THEN_CELLS) )
227
  {
228
    vtkErrorMacro("Input Array must be associated with points or cells.");
229
    return;
230
  }
231 232 233 234

  this->SetInputArrayToProcess(0, 0, 0, fieldAssociation, name);
}

235
//-----------------------------------------------------------------------------
236 237 238 239 240 241
void vtkGradientFilter::SetInputScalars(int fieldAssociation,
                                        int fieldAttributeType)
{
  if (   (fieldAssociation != vtkDataObject::FIELD_ASSOCIATION_POINTS)
      && (fieldAssociation != vtkDataObject::FIELD_ASSOCIATION_CELLS)
      && (fieldAssociation != vtkDataObject::FIELD_ASSOCIATION_POINTS_THEN_CELLS) )
242
  {
243
    vtkErrorMacro("Input Array must be associated with points or cells.");
244
    return;
245
  }
246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262

  this->SetInputArrayToProcess(0, 0, 0, fieldAssociation, fieldAttributeType);
}

//-----------------------------------------------------------------------------
int vtkGradientFilter::RequestUpdateExtent(vtkInformation *vtkNotUsed(request),
                                           vtkInformationVector **inputVector,
                                           vtkInformationVector *outputVector)
{
  // get the info objects
  vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
  vtkInformation *outInfo = outputVector->GetInformationObject(0);

  // Technically, this code is only correct for pieces extent types.  However,
  // since this class is pretty inefficient for data types that use 3D extents,
  // we'll punt on the ghost levels for them, too.
  int piece, numPieces, ghostLevels;
263

264 265 266 267 268
  piece = outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER());
  numPieces = outInfo->Get(
                   vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_PIECES());
  ghostLevels = outInfo->Get(
             vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS());
269

270
  if (numPieces > 1)
271
  {
272
    ++ghostLevels;
273
  }
274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298

  inInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER(), piece);
  inInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_PIECES(),
              numPieces);
  inInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS(),
              ghostLevels);
  inInfo->Set(vtkStreamingDemandDrivenPipeline::EXACT_EXTENT(), 1);

  return 1;
}

//-----------------------------------------------------------------------------
int vtkGradientFilter::RequestData(vtkInformation *vtkNotUsed(request),
                                   vtkInformationVector **inputVector,
                                   vtkInformationVector *outputVector)
{
  vtkDebugMacro("RequestData");

  vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
  vtkInformation *outInfo = outputVector->GetInformationObject(0);

  vtkDataSet *input
    = vtkDataSet::SafeDownCast(inInfo->Get(vtkDataObject::DATA_OBJECT()));
  vtkDataSet *output
    = vtkDataSet::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()));
299
  vtkDataArray *array = this->GetInputArrayToProcess(0, inputVector);
300

301 302 303 304 305 306 307 308 309 310 311 312 313 314
  if (input->GetNumberOfCells() == 0)
  {
    // need cells to compute the gradient so if we don't have cells. we can't compute anything.
    // if we have points and an array with values provide a warning letting the user know that
    // no gradient will be computed because of the lack of cells. otherwise the dataset is
    // assumed empty and we can skip providing a warning message to the user.
    if (input->GetNumberOfPoints() && array && array->GetNumberOfTuples())
    {
      vtkWarningMacro("Cannot compute gradient for datasets without cells");
    }
    output->ShallowCopy(input);
    return 1;
  }

315
  if (array == nullptr)
316
  {
317 318
    vtkErrorMacro("No input array. If this dataset is part of a composite dataset"
                  << " check to make sure that all non-empty blocks have this array.");
319
    return 0;
320
  }
321
  if (array->GetNumberOfComponents() == 0)
322
  {
323
    vtkErrorMacro("Input array must have at least one component.");
324
    return 0;
325
  }
326

327 328 329 330
  // we can only compute vorticity and Q criterion if the input
  // array has 3 components. if we can't compute them because of
  // this we only mark internally the we aren't computing them
  // since we don't want to change the state of the filter.
331
  bool computeVorticity = this->ComputeVorticity != 0;
332
  bool computeDivergence = this->ComputeDivergence != 0;
333
  bool computeQCriterion = this->ComputeQCriterion != 0;
334 335
  if( (this->ComputeQCriterion || this->ComputeVorticity ||
       this->ComputeDivergence) && array->GetNumberOfComponents() != 3)
336
  {
337 338 339
    vtkWarningMacro("Input array must have exactly three components with "
                    << "ComputeDivergence, ComputeVorticity or ComputeQCriterion flag enabled."
                    << "Skipping divergence, vorticity and Q-criterion computation.");
340 341
    computeVorticity = false;
    computeQCriterion = false;
342
    computeDivergence = false;
343
  }
344

345
  int fieldAssociation;
346
  if (vtkGradientFilterHasArray(input->GetPointData(), array))
347
  {
348
    fieldAssociation = vtkDataObject::FIELD_ASSOCIATION_POINTS;
349
  }
350
  else if (vtkGradientFilterHasArray(input->GetCellData(), array))
351
  {
352
    fieldAssociation = vtkDataObject::FIELD_ASSOCIATION_CELLS;
353
  }
354
  else
355
  {
356
    vtkErrorMacro("Input arrays do not seem to be either point or cell arrays.");
357
    return 0;
358
  }
359 360 361 362 363

  output->CopyStructure(input);
  output->GetPointData()->PassData(input->GetPointData());
  output->GetCellData()->PassData(input->GetCellData());

364 365
  if(output->IsA("vtkImageData") || output->IsA("vtkStructuredGrid") ||
          output->IsA("vtkRectilinearGrid") )
366
  {
367
    this->ComputeRegularGridGradient(
368 369
      array, fieldAssociation, computeVorticity, computeQCriterion,
      computeDivergence, output);
370
  }
371
  else
372
  {
373
    this->ComputeUnstructuredGridGradient(
374 375
      array, fieldAssociation, input, computeVorticity, computeQCriterion,
      computeDivergence, output);
376
  }
377 378 379 380 381 382

  return 1;
}

//-----------------------------------------------------------------------------
int vtkGradientFilter::ComputeUnstructuredGridGradient(
383
  vtkDataArray* array, int fieldAssociation, vtkDataSet* input,
384 385
  bool computeVorticity, bool computeQCriterion, bool computeDivergence,
  vtkDataSet* output)
386
{
387
  int arrayType = this->GetOutputArrayType(array);
388
  int numberOfInputComponents = array->GetNumberOfComponents();
389
  vtkSmartPointer<vtkDataArray> gradients = nullptr;
390
  if(this->ComputeGradient)
391
  {
392
    gradients.TakeReference(vtkDataArray::CreateDataArray(arrayType));
393 394
    gradients->SetNumberOfComponents(3*numberOfInputComponents);
    gradients->SetNumberOfTuples(array->GetNumberOfTuples());
395 396 397 398
    switch (arrayType)
    {
      vtkFloatingPointTemplateMacro(Fill(gradients, static_cast<VTK_TT>(0), this->ReplacementValueOption));
    }
399
    if (this->ResultArrayName)
400
    {
401
      gradients->SetName(this->ResultArrayName);
402
    }
403
    else
404
    {
405
      gradients->SetName("Gradients");
406
    }
407
  }
408
  vtkSmartPointer<vtkDataArray> divergence = nullptr;
409
  if(computeDivergence)
410
  {
411
    divergence.TakeReference(vtkDataArray::CreateDataArray(arrayType));
412
    divergence->SetNumberOfTuples(array->GetNumberOfTuples());
413 414
    switch (arrayType)
    {
415
      vtkFloatingPointTemplateMacro(Fill(divergence, static_cast<VTK_TT>(0), this->ReplacementValueOption));
416
    }
417
    if (this->DivergenceArrayName)
418
    {
419
      divergence->SetName(this->DivergenceArrayName);
420
    }
421
    else
422
    {
423
      divergence->SetName("Divergence");
424
    }
425
  }
426 427
  vtkSmartPointer<vtkDataArray> vorticity;
  if(computeVorticity)
428
  {
429
    vorticity.TakeReference(vtkDataArray::CreateDataArray(arrayType));
430 431
    vorticity->SetNumberOfComponents(3);
    vorticity->SetNumberOfTuples(array->GetNumberOfTuples());
432 433
    switch (arrayType)
    {
434
      vtkFloatingPointTemplateMacro(Fill(vorticity, static_cast<VTK_TT>(0), this->ReplacementValueOption));
435
    }
436
    if (this->VorticityArrayName)
437
    {
438
      vorticity->SetName(this->VorticityArrayName);
439
    }
440
    else
441
    {
442 443
      vorticity->SetName("Vorticity");
    }
444
  }
445
  vtkSmartPointer<vtkDataArray> qCriterion;
446
  if(computeQCriterion)
447
  {
448
    qCriterion.TakeReference(vtkDataArray::CreateDataArray(arrayType));
449
    qCriterion->SetNumberOfTuples(array->GetNumberOfTuples());
450 451
    switch (arrayType)
    {
452
      vtkFloatingPointTemplateMacro(Fill(qCriterion, static_cast<VTK_TT>(0), this->ReplacementValueOption));
453
    }
454
    if (this->QCriterionArrayName)
455
    {
456
      qCriterion->SetName(this->QCriterionArrayName);
457
    }
458
    else
459
    {
460
      qCriterion->SetName("Q-criterion");
461
    }
462
  }
463

464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481
  int highestCellDimension = 0;
  if (this->ContributingCellOption == vtkGradientFilter::DataSetMax)
  {
    int maxDimension = input->IsA("vtkPolyData") == 1 ? 2 : 3;
    for (vtkIdType i=0;i<input->GetNumberOfCells();i++)
    {
      int dim = input->GetCell(i)->GetCellDimension();
      if (dim > highestCellDimension)
      {
        highestCellDimension = dim;
        if (highestCellDimension == maxDimension)
        {
          break;
        }
      }
    }
  }

482
  if (fieldAssociation == vtkDataObject::FIELD_ASSOCIATION_POINTS)
483
  {
484
    if (!this->FasterApproximation)
485
    {
486 487 488 489
      switch (arrayType)
      { // ok to use template macro here since we made the output arrays ourselves
        vtkFloatingPointTemplateMacro(ComputePointGradientsUG(
                           input, array,
490
                           (gradients == nullptr ? nullptr :
491
                            static_cast<VTK_TT *>(gradients->GetVoidPointer(0))),
492
                           numberOfInputComponents,
493
                           (vorticity == nullptr ? nullptr :
494
                            static_cast<VTK_TT *>(vorticity->GetVoidPointer(0))),
495
                           (qCriterion == nullptr ? nullptr :
496
                            static_cast<VTK_TT *>(qCriterion->GetVoidPointer(0))),
497
                           (divergence == nullptr ? nullptr :
498 499
                            static_cast<VTK_TT *>(divergence->GetVoidPointer(0))),
                           highestCellDimension, this->ContributingCellOption));
500
      }
501
      if(gradients)
502
      {
503
        output->GetPointData()->AddArray(gradients);
504
      }
505
      if(divergence)
506
      {
507
        output->GetPointData()->AddArray(divergence);
508
      }
509
      if(vorticity)
510
      {
511
        output->GetPointData()->AddArray(vorticity);
512
      }
513
      if(qCriterion)
514
      {
515
        output->GetPointData()->AddArray(qCriterion);
516
      }
517
    }
518
    else // this->FasterApproximation
519
    {
520 521 522
      // The cell computation is faster and works off of point data anyway.  The
      // faster approximation is to use the cell algorithm and then convert the
      // result to point data.
523
      vtkSmartPointer<vtkDataArray> cellGradients = nullptr;
524
      if(gradients)
525
      {
526
        cellGradients.TakeReference(vtkDataArray::CreateDataArray(arrayType));
527 528 529
        cellGradients->SetName(gradients->GetName());
        cellGradients->SetNumberOfComponents(3*array->GetNumberOfComponents());
        cellGradients->SetNumberOfTuples(input->GetNumberOfCells());
530
      }
531
      vtkSmartPointer<vtkDataArray> cellDivergence = nullptr;
532
      if(divergence)
533
      {
534
        cellDivergence.TakeReference(vtkDataArray::CreateDataArray(arrayType));
535 536
        cellDivergence->SetName(divergence->GetName());
        cellDivergence->SetNumberOfTuples(input->GetNumberOfCells());
537
      }
538
      vtkSmartPointer<vtkDataArray> cellVorticity = nullptr;
539
      if(vorticity)
540
      {
541
        cellVorticity.TakeReference(vtkDataArray::CreateDataArray(arrayType));
542 543 544
        cellVorticity->SetName(vorticity->GetName());
        cellVorticity->SetNumberOfComponents(3);
        cellVorticity->SetNumberOfTuples(input->GetNumberOfCells());
545
      }
546
      vtkSmartPointer<vtkDataArray> cellQCriterion = nullptr;
547
      if(qCriterion)
548
      {
549
        cellQCriterion.TakeReference(vtkDataArray::CreateDataArray(arrayType));
550 551
        cellQCriterion->SetName(qCriterion->GetName());
        cellQCriterion->SetNumberOfTuples(input->GetNumberOfCells());
552
      }
553

554 555 556
      switch (arrayType)
      { // ok to use template macro here since we made the output arrays ourselves
        vtkFloatingPointTemplateMacro(
557
          ComputeCellGradientsUG(
558
            input, array,
559
            (cellGradients == nullptr ? nullptr :
560
             static_cast<VTK_TT *>(cellGradients->GetVoidPointer(0))),
561
            numberOfInputComponents,
562
            (vorticity == nullptr ? nullptr :
563
             static_cast<VTK_TT *>(cellVorticity->GetVoidPointer(0))),
564
            (qCriterion == nullptr ? nullptr :
565
             static_cast<VTK_TT *>(cellQCriterion->GetVoidPointer(0))),
566
            (divergence == nullptr ? nullptr :
567
             static_cast<VTK_TT *>(cellDivergence->GetVoidPointer(0)))));
568
      }
569

570
      // We need to convert cell Array to points Array.
571 572
      vtkSmartPointer<vtkDataSet> dummy;
      dummy.TakeReference(input->NewInstance());
573
      dummy->CopyStructure(input);
574
      if(cellGradients)
575
      {
576
        dummy->GetCellData()->AddArray(cellGradients);
577
      }
578
      if(divergence)
579
      {
580
        dummy->GetCellData()->AddArray(cellDivergence);
581
      }
582
      if(vorticity)
583
      {
584
        dummy->GetCellData()->AddArray(cellVorticity);
585
      }
586
      if(qCriterion)
587
      {
588
        dummy->GetCellData()->AddArray(cellQCriterion);
589
      }
590

591
      vtkNew<vtkCellDataToPointData> cd2pd;
592
      cd2pd->SetInputData(dummy);
593
      cd2pd->PassCellDataOff();
594
      cd2pd->SetContributingCellOption(this->ContributingCellOption);
595 596 597
      cd2pd->Update();

      // Set the gradients array in the output and cleanup.
598
      if(gradients)
599
      {
600 601
        output->GetPointData()->AddArray(
          cd2pd->GetOutput()->GetPointData()->GetArray(gradients->GetName()));
602
      }
603
      if(qCriterion)
604
      {
605 606
        output->GetPointData()->AddArray(
          cd2pd->GetOutput()->GetPointData()->GetArray(qCriterion->GetName()));
607
      }
608
      if(divergence)
609
      {
610 611
        output->GetPointData()->AddArray(
          cd2pd->GetOutput()->GetPointData()->GetArray(divergence->GetName()));
612
      }
613
      if(vorticity)
614
      {
615 616
        output->GetPointData()->AddArray(
          cd2pd->GetOutput()->GetPointData()->GetArray(vorticity->GetName()));
617 618
      }
    }
619
  }
luz.paz's avatar
luz.paz committed
620
  else  // fieldAssociation == vtkDataObject::FIELD_ASSOCIATION_CELLS
621
  {
622
    // We need to convert cell Array to points Array.
623 624
    vtkSmartPointer<vtkDataSet> dummy;
    dummy.TakeReference(input->NewInstance());
625
    dummy->CopyStructure(input);
626
    dummy->GetCellData()->SetScalars(array);
627

628
    vtkNew<vtkCellDataToPointData> cd2pd;
629
    cd2pd->SetInputData(dummy);
630
    cd2pd->PassCellDataOff();
631
    cd2pd->SetContributingCellOption(this->ContributingCellOption);
632 633 634 635 636
    cd2pd->Update();
    vtkDataArray *pointScalars
      = cd2pd->GetOutput()->GetPointData()->GetScalars();
    pointScalars->Register(this);

637 638 639 640
    switch (arrayType)
    { // ok to use template macro here since we made the output arrays ourselves
      vtkFloatingPointTemplateMacro(ComputeCellGradientsUG(
                         input, pointScalars,
641
                         (gradients == nullptr ? nullptr :
642
                          static_cast<VTK_TT *>(gradients->GetVoidPointer(0))),
643
                         numberOfInputComponents,
644
                         (vorticity == nullptr ? nullptr :
645
                          static_cast<VTK_TT *>(vorticity->GetVoidPointer(0))),
646
                         (qCriterion == nullptr ? nullptr :
647
                          static_cast<VTK_TT *>(qCriterion->GetVoidPointer(0))),
648
                         (divergence == nullptr ? nullptr :
649
                          static_cast<VTK_TT *>(divergence->GetVoidPointer(0)))));
650
    }
651

652
    if(gradients)
653
    {
654
      output->GetCellData()->AddArray(gradients);
655
    }
656
    if(vorticity)
657
    {
658
      output->GetCellData()->AddArray(vorticity);
659
    }
660
    if(divergence)
661
    {
662
      output->GetCellData()->AddArray(divergence);
663
    }
664
    if(qCriterion)
665
    {
666
      output->GetCellData()->AddArray(qCriterion);
667
    }
668 669
    pointScalars->UnRegister(this);
  }
670 671 672 673 674

  return 1;
}

//-----------------------------------------------------------------------------
675
int vtkGradientFilter::ComputeRegularGridGradient(
676
  vtkDataArray* array, int fieldAssociation, bool computeVorticity,
677
  bool computeQCriterion, bool computeDivergence, vtkDataSet* output)
678
{
679
  int arrayType = this->GetOutputArrayType(array);
680
  int numberOfInputComponents = array->GetNumberOfComponents();
681
  vtkSmartPointer<vtkDataArray> gradients = nullptr;
682
  if(this->ComputeGradient)
683
  {
684
    gradients.TakeReference(vtkDataArray::CreateDataArray(arrayType));
685 686 687
    gradients->SetNumberOfComponents(3*numberOfInputComponents);
    gradients->SetNumberOfTuples(array->GetNumberOfTuples());
    if (this->ResultArrayName)
688
    {
689
      gradients->SetName(this->ResultArrayName);
690
    }
691
    else
692
    {
693
      gradients->SetName("Gradients");
694
    }
695
  }
696
  vtkSmartPointer<vtkDataArray> divergence = nullptr;
697
  if(computeDivergence)
698
  {
699
    divergence.TakeReference(vtkDataArray::CreateDataArray(arrayType));
700 701
    divergence->SetNumberOfTuples(array->GetNumberOfTuples());
    if (this->DivergenceArrayName)
702
    {
703
      divergence->SetName(this->DivergenceArrayName);
704
    }
705
    else
706
    {
707
      divergence->SetName("Divergence");
708
    }
709
  }
710 711
  vtkSmartPointer<vtkDataArray> vorticity;
  if(computeVorticity)
712
  {
713
    vorticity.TakeReference(vtkDataArray::CreateDataArray(arrayType));
714 715 716
    vorticity->SetNumberOfComponents(3);
    vorticity->SetNumberOfTuples(array->GetNumberOfTuples());
    if (this->VorticityArrayName)
717
    {
718
      vorticity->SetName(this->VorticityArrayName);
719
    }
720
    else
721
    {
722 723
      vorticity->SetName("Vorticity");
    }
724
  }
725
  vtkSmartPointer<vtkDataArray> qCriterion;
726
  if(computeQCriterion)
727
  {
728
    qCriterion.TakeReference(vtkDataArray::CreateDataArray(arrayType));
729
    qCriterion->SetNumberOfTuples(array->GetNumberOfTuples());
730
    if (this->QCriterionArrayName)
731
    {
732
      qCriterion->SetName(this->QCriterionArrayName);
733
    }
734
    else
735
    {
736
      qCriterion->SetName("Q-criterion");
737
    }
738
  }
739

740
  if(vtkStructuredGrid* structuredGrid = vtkStructuredGrid::SafeDownCast(output))
741
  {
742 743 744 745
    switch (arrayType)
    { // ok to use template macro here since we made the output arrays ourselves
      vtkFloatingPointTemplateMacro(ComputeGradientsSG(
                         structuredGrid, array,
746
                         (gradients == nullptr ? nullptr :
747
                          static_cast<VTK_TT *>(gradients->GetVoidPointer(0))),
748
                         numberOfInputComponents, fieldAssociation,
749
                         (vorticity == nullptr ? nullptr :
750
                          static_cast<VTK_TT *>(vorticity->GetVoidPointer(0))),
751
                         (qCriterion == nullptr ? nullptr :
752
                          static_cast<VTK_TT *>(qCriterion->GetVoidPointer(0))),
753
                         (divergence == nullptr ? nullptr :
754 755
                          static_cast<VTK_TT *>(divergence->GetVoidPointer(0)))));

756
    }
757
  }
758
  else if(vtkImageData* imageData = vtkImageData::SafeDownCast(output))
759
  {
760 761 762 763
    switch (arrayType)
    { // ok to use template macro here since we made the output arrays ourselves
      vtkFloatingPointTemplateMacro(ComputeGradientsSG(
                         imageData, array,
764
                         (gradients == nullptr ? nullptr :
765
                          static_cast<VTK_TT *>(gradients->GetVoidPointer(0))),
766
                         numberOfInputComponents, fieldAssociation,
767
                         (vorticity == nullptr ? nullptr :
768
                          static_cast<VTK_TT *>(vorticity->GetVoidPointer(0))),
769
                         (qCriterion == nullptr ? nullptr :
770
                          static_cast<VTK_TT *>(qCriterion->GetVoidPointer(0))),
771
                         (divergence == nullptr ? nullptr :
772
                          static_cast<VTK_TT *>(divergence->GetVoidPointer(0)))));
773
    }
774
  }
775
  else if(vtkRectilinearGrid* rectilinearGrid = vtkRectilinearGrid::SafeDownCast(output))
776
  {
777 778 779 780
    switch (arrayType)
    { // ok to use template macro here since we made the output arrays ourselves
      vtkFloatingPointTemplateMacro(ComputeGradientsSG(
                         rectilinearGrid, array,
781
                         (gradients == nullptr ? nullptr :
782
                          static_cast<VTK_TT *>(gradients->GetVoidPointer(0))),
783
                         numberOfInputComponents, fieldAssociation,
784
                         (vorticity == nullptr ? nullptr :
785
                          static_cast<VTK_TT *>(vorticity->GetVoidPointer(0))),
786
                         (qCriterion == nullptr ? nullptr :
787
                          static_cast<VTK_TT *>(qCriterion->GetVoidPointer(0))),
788
                         (divergence == nullptr ? nullptr :
789 790
                          static_cast<VTK_TT *>(divergence->GetVoidPointer(0)))));

791
    }
792
  }
793
  if(fieldAssociation == vtkDataObject::FIELD_ASSOCIATION_POINTS)
794
  {
795
    if(gradients)
796
    {
797
      output->GetPointData()->AddArray(gradients);
798
    }
799
    if(vorticity)
800
    {
801
      output->GetPointData()->AddArray(vorticity);
802
    }
803
    if(qCriterion)
804
    {
805
      output->GetPointData()->AddArray(qCriterion);
806
    }
807
    if(divergence)
808
    {
809
      output->GetPointData()->AddArray(divergence);
810
    }
811
  }
812
  else if(fieldAssociation == vtkDataObject::FIELD_ASSOCIATION_CELLS)
813
  {
814
    if(gradients)
815
    {
816
      output->GetCellData()->AddArray(gradients);
817
    }
818
    if(vorticity)
819
    {
820
      output->GetCellData()->AddArray(vorticity);
821
    }
822
    if(qCriterion)
823
    {
824
      output->GetCellData()->AddArray(qCriterion);
825
    }
826
    if(divergence)
827
    {
828
      output->GetCellData()->AddArray(divergence);
829
    }
830
  }
831
  else
832
  {
833
    vtkErrorMacro("Bad fieldAssociation value " << fieldAssociation << endl);
834
  }
835

836 837 838
  return 1;
}

839 840 841 842 843 844 845 846 847 848 849
//-----------------------------------------------------------------------------
int vtkGradientFilter::GetOutputArrayType(vtkDataArray* array)
{
  int retType = VTK_DOUBLE;
  switch (array->GetDataType())
  {
    vtkTemplateMacro(retType = GetOutputDataType(static_cast<VTK_TT>(0)));
  }
  return retType;
}

850
namespace {
851
//-----------------------------------------------------------------------------
852 853
  template<class data_type>
  void ComputePointGradientsUG(
854
    vtkDataSet *structure, vtkDataArray *array, data_type *gradients,
855
    int numberOfInputComponents, data_type* vorticity, data_type* qCriterion,
856
    data_type* divergence, int highestCellDimension, int contributingCellOption)
857
  {
858
    vtkNew<vtkIdList> currentPoint;
859
    currentPoint->SetNumberOfIds(1);
860
    vtkNew<vtkIdList> cellsOnPoint;
861

862
    vtkIdType numpts = structure->GetNumberOfPoints();
863

864
    int numberOfOutputComponents = 3*numberOfInputComponents;
865
    std::vector<data_type> g(numberOfOutputComponents);
866

867 868 869 870
    // if we are doing patches for contributing cell dimensions we want to keep track of
    // the maximum expected dimension so we can exit out of the check loop quicker
    const int maxCellDimension = structure->IsA("vtkPolyData") ? 2 : 3;

871
    for (vtkIdType point = 0; point < numpts; point++)
872
    {
873 874 875 876
      currentPoint->SetId(0, point);
      double pointcoords[3];
      structure->GetPoint(point, pointcoords);
      // Get all cells touching this point.
877
      structure->GetCellNeighbors(-1, currentPoint, cellsOnPoint);
878 879
      vtkIdType numCellNeighbors = cellsOnPoint->GetNumberOfIds();

880
      for(int i=0;i<numberOfOutputComponents;i++)
881
      {
882
        g[i] = 0;
883
      }
884

885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902
      if (contributingCellOption == vtkGradientFilter::Patch)
      {
        highestCellDimension = 0;
        for (vtkIdType neighbor = 0; neighbor < numCellNeighbors; neighbor++)
        {
          int cellDimension = structure->GetCell(cellsOnPoint->GetId(neighbor))->GetCellDimension();
          if (cellDimension > highestCellDimension)
          {
            highestCellDimension = cellDimension;
            if (highestCellDimension == maxCellDimension)
            {
              break;
            }
          }
        }
      }
      vtkIdType numValidCellNeighbors = 0;

903 904 905
      // Iterate on all cells and find all points connected to current point
      // by an edge.
      for (vtkIdType neighbor = 0; neighbor < numCellNeighbors; neighbor++)
906
      {
907
        vtkCell *cell = structure->GetCell(cellsOnPoint->GetId(neighbor));
908
        if (cell->GetCellDimension() >= highestCellDimension)
909
        {
910 911 912 913
          int subId;
          double parametricCoord[3];
          if(GetCellParametricData(point, pointcoords, cell,
                                   subId, parametricCoord))
914
          {
915 916
            numValidCellNeighbors++;
            for(int inputComponent=0;inputComponent<numberOfInputComponents;inputComponent++)
917
            {
918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935
              int numberOfCellPoints = cell->GetNumberOfPoints();
              std::vector<double> values(numberOfCellPoints);
              // Get values of Array at cell points.
              for (int i = 0; i < numberOfCellPoints; i++)
              {
                values[i] = array->GetComponent(cell->GetPointId(i), inputComponent);
              }

              double derivative[3];
              // Get derivative of cell at point.
              cell->Derivatives(subId, parametricCoord, &values[0], 1, derivative);

              g[inputComponent*3] += static_cast<data_type>(derivative[0]);
              g[inputComponent*3+1] += static_cast<data_type>(derivative[1]);
              g[inputComponent*3+2] += static_cast<data_type>(derivative[2]);
            } // iterating over Components
          } // if(GetCellParametricData())
        } // if(cell->GetCellDimension () >= highestCellDimension
936
      } // iterating over neighbors
937

938
      if (numValidCellNeighbors > 0)
939
      {
940
        for(int i=0;i<3*numberOfInputComponents;i++)
941
        {
942
          g[i] /= numValidCellNeighbors;
943
        }
944

945
        if(vorticity)
946
        {
947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962
          ComputeVorticityFromGradient(&g[0], vorticity+3*point);
        }
        if(qCriterion)
        {
          ComputeQCriterionFromGradient(&g[0], qCriterion+point);
        }
        if(divergence)
        {
          ComputeDivergenceFromGradient(&g[0], divergence+point);
        }
        if(gradients)
        {
          for(int i=0;i<numberOfOutputComponents;i++)
          {
            gradients[point*numberOfOutputComponents+i] = g[i];
          }
963
        }
964 965
      }
    }  // iterating over points in grid
966
  }
967

968
//-----------------------------------------------------------------------------
969
  int GetCellParametricData(vtkIdType pointId, double pointCoord[3],
970 971 972 973 974 975 976
                            vtkCell *cell, int &subId, double parametricCoord[3])
  {
    // Watch out for degenerate cells.  They make the derivative calculation
    // fail.
    vtkIdList *pointIds = cell->GetPointIds();
    int timesPointRegistered = 0;
    for (int i = 0; i < pointIds->GetNumberOfIds(); i++)
977
    {
978
      if (pointId == pointIds->GetId(i))
979
      {
980 981
        timesPointRegistered++;
      }
982
    }
983
    if (timesPointRegistered != 1)
984
    {
985 986
      // The cell should have the point exactly once.  Not good.
      return