vtkmGradient.cxx 9.3 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
//=============================================================================
//
//  Copyright (c) Kitware, Inc.
//  All rights reserved.
//  See LICENSE.txt 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 2012 Sandia Corporation.
//  Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
//  the U.S. Government retains certain rights in this software.
//
//=============================================================================
#include "vtkmGradient.h"

18
#include "vtkCellData.h"
19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41
#include "vtkDataSet.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkNew.h"
#include "vtkObjectFactory.h"
#include "vtkPointData.h"

#include "vtkmlib/ArrayConverters.h"
#include "vtkmlib/DataSetConverters.h"
#include "vtkmlib/PolyDataConverter.h"
#include "vtkmlib/Storage.h"

#include "vtkmCellSetExplicit.h"
#include "vtkmCellSetSingleType.h"
#include "vtkmFilterPolicy.h"

#include <vtkm/filter/Gradient.h>
#include <vtkm/filter/PointAverage.h>

vtkStandardNewMacro(vtkmGradient)

namespace
{
42
struct GradientTypes
43
    : vtkm::ListTagBase<
44 45
                        vtkm::Float32,
                        vtkm::Float64,
46 47 48 49 50 51 52 53 54
                        vtkm::Vec<vtkm::Float32,3>,
                        vtkm::Vec<vtkm::Float64,3>,
                        vtkm::Vec< vtkm::Vec<vtkm::Float32,3>, 3>,
                        vtkm::Vec< vtkm::Vec<vtkm::Float64,3>, 3>
                        >
{
};

//------------------------------------------------------------------------------
55 56
class vtkmGradientFilterPolicy
      : public vtkm::filter::PolicyBase<vtkmGradientFilterPolicy>
57 58
  {
  public:
59
    using FieldTypeList = GradientTypes;
60

61 62 63
    using StructuredCellSetList = tovtkm::CellListStructuredInVTK;
    using UnstructuredCellSetList = tovtkm::CellListUnstructuredInVTK;
    using AllCellSetList = tovtkm::CellListAllInVTK;
64
  };
Sujin Philip's avatar
Sujin Philip committed
65 66 67 68 69 70 71 72 73 74 75 76 77

vtkm::cont::DataSet CopyDataSetStructure(const vtkm::cont::DataSet& ds)
{
  vtkm::cont::DataSet cp;
  for (vtkm::IdComponent i = 0; i < ds.GetNumberOfCoordinateSystems(); ++i)
  {
    cp.AddCoordinateSystem(ds.GetCoordinateSystem(i));
  }
  for (vtkm::IdComponent i = 0; i < ds.GetNumberOfCellSets(); ++i)
  {
    cp.AddCellSet(ds.GetCellSet(i));
  }
  return cp;
78 79
}

Sujin Philip's avatar
Sujin Philip committed
80 81
} // anonymous namespace

82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99
//------------------------------------------------------------------------------
vtkmGradient::vtkmGradient()
{
}

//------------------------------------------------------------------------------
vtkmGradient::~vtkmGradient()
{
}

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

//------------------------------------------------------------------------------
int vtkmGradient::RequestData(vtkInformation* request,
100 101
                              vtkInformationVector** inputVector,
                              vtkInformationVector* outputVector)
102 103 104 105 106 107 108 109 110 111 112 113 114 115 116
{
  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()));

  output->ShallowCopy(input);

  // grab the input array to process to determine the field want to compute
  // the gradient for
  int association = this->GetInputArrayAssociation(0, inputVector);
  vtkDataArray* inputArray = this->GetInputArrayToProcess(0, inputVector);
Sujin Philip's avatar
Sujin Philip committed
117 118
  if (inputArray == nullptr || inputArray->GetName() == nullptr ||
      inputArray->GetName()[0] == '\0')
119
  {
Sujin Philip's avatar
Sujin Philip committed
120 121
    vtkErrorMacro("Invalid input array.");
    return 0;
122
  }
123

Sujin Philip's avatar
Sujin Philip committed
124
  try
125
  {
Sujin Philip's avatar
Sujin Philip committed
126 127 128 129 130 131 132
    // convert the input dataset to a vtkm::cont::DataSet. We explicitly drop
    // all arrays from the conversion as this algorithm doesn't change topology
    // and therefore doesn't need input fields converted through the VTK-m filter
    auto in = tovtkm::Convert(input, tovtkm::FieldsFlag::None);
    vtkm::cont::Field field = tovtkm::Convert(inputArray, association);
    in.AddField(field);

133 134
    const bool fieldIsPoint = field.GetAssociation() == vtkm::cont::Field::Association::POINTS;
    const bool fieldIsCell = field.GetAssociation() == vtkm::cont::Field::Association::CELL_SET;
Sujin Philip's avatar
Sujin Philip committed
135 136 137 138 139 140 141 142 143 144 145 146 147
    const bool fieldIsVec = (inputArray->GetNumberOfComponents() == 3);
    const bool fieldIsScalar = inputArray->GetDataType() == VTK_FLOAT ||
                              inputArray->GetDataType() == VTK_DOUBLE;
    const bool fieldValid = (fieldIsPoint || fieldIsCell) &&
                            fieldIsScalar &&
                            (field.GetName() != std::string());

    if (!fieldValid)
    {
      vtkWarningMacro(<< "Unsupported field type\n"
                      << "Falling back to vtkGradientFilter.");
      return this->Superclass::RequestData(request, inputVector, outputVector);
    }
148

Sujin Philip's avatar
Sujin Philip committed
149 150 151 152 153 154 155 156 157 158 159 160 161
    vtkmGradientFilterPolicy policy;
    auto passNoFields =
      vtkm::filter::FieldSelection(vtkm::filter::FieldSelection::MODE_NONE);
    vtkm::filter::Gradient filter;
    filter.SetFieldsToPass(passNoFields);
    filter.SetColumnMajorOrdering();

    if (fieldIsVec)
    { //this properties are only valid when processing a vec<3> field
      filter.SetComputeDivergence(this->ComputeDivergence != 0);
      filter.SetComputeVorticity(this->ComputeVorticity != 0);
      filter.SetComputeQCriterion(this->ComputeQCriterion != 0);
    }
162

Sujin Philip's avatar
Sujin Philip committed
163 164 165 166
    if (this->ResultArrayName)
    {
      filter.SetOutputFieldName(this->ResultArrayName);
    }
167

Sujin Philip's avatar
Sujin Philip committed
168 169 170 171
    if (this->DivergenceArrayName)
    {
      filter.SetDivergenceName(this->DivergenceArrayName);
    }
172

Sujin Philip's avatar
Sujin Philip committed
173 174 175 176
    if (this->VorticityArrayName)
    {
      filter.SetVorticityName(this->VorticityArrayName);
    }
177

Sujin Philip's avatar
Sujin Philip committed
178 179 180 181 182 183 184 185
    if (this->QCriterionArrayName)
    {
      filter.SetQCriterionName(this->QCriterionArrayName);
    }
    else
    {
      filter.SetQCriterionName("Q-criterion");
    }
186

Sujin Philip's avatar
Sujin Philip committed
187 188 189 190 191 192
    // Run the VTK-m Gradient Filter
    // -----------------------------
    vtkm::cont::DataSet result;
    if (fieldIsPoint)
    {
      filter.SetComputePointGradient(!this->FasterApproximation);
193
      filter.SetActiveField(field.GetName(), vtkm::cont::Field::Association::POINTS);
Sujin Philip's avatar
Sujin Philip committed
194 195 196 197 198 199
      result = filter.Execute(in, policy);

      //When we have faster approximation enabled the VTK-m gradient will output
      //a cell field not a point field. So at that point we will need to convert
      //back to a point field
      if (this->FasterApproximation)
200
      {
Sujin Philip's avatar
Sujin Philip committed
201 202 203 204 205 206 207 208 209
        vtkm::filter::PointAverage cellToPoint;
        cellToPoint.SetFieldsToPass(passNoFields);

        auto c2pIn = result;
        result = CopyDataSetStructure(result);

        if (this->ComputeGradient)
        {
          cellToPoint.SetActiveField(
210
            filter.GetOutputFieldName(), vtkm::cont::Field::Association::CELL_SET);
Sujin Philip's avatar
Sujin Philip committed
211 212 213 214 215 216
          auto ds = cellToPoint.Execute(c2pIn, policy);
          result.AddField(ds.GetField(0));
        }
        if (this->ComputeDivergence && fieldIsVec)
        {
          cellToPoint.SetActiveField(
217
            filter.GetDivergenceName(), vtkm::cont::Field::Association::CELL_SET);
Sujin Philip's avatar
Sujin Philip committed
218 219 220 221 222 223
          auto ds = cellToPoint.Execute(c2pIn, policy);
          result.AddField(ds.GetField(0));
        }
        if (this->ComputeVorticity && fieldIsVec)
        {
          cellToPoint.SetActiveField(filter.GetVorticityName(),
224
            vtkm::cont::Field::Association::CELL_SET);
Sujin Philip's avatar
Sujin Philip committed
225 226 227 228 229 230
          auto ds = cellToPoint.Execute(c2pIn, policy);
          result.AddField(ds.GetField(0));
        }
        if (this->ComputeQCriterion && fieldIsVec)
        {
          cellToPoint.SetActiveField(filter.GetQCriterionName(),
231
            vtkm::cont::Field::Association::CELL_SET);
Sujin Philip's avatar
Sujin Philip committed
232 233 234
          auto ds = cellToPoint.Execute(c2pIn, policy);
          result.AddField(ds.GetField(0));
        }
235
      }
Sujin Philip's avatar
Sujin Philip committed
236 237 238 239 240 241 242 243 244 245 246
    }
    else
    {
      //we need to convert the field to be a point field
      vtkm::filter::PointAverage cellToPoint;
      cellToPoint.SetFieldsToPass(passNoFields);
      cellToPoint.SetActiveField(field.GetName(), field.GetAssociation());
      cellToPoint.SetOutputFieldName(field.GetName());
      in = cellToPoint.Execute(in, policy);

      filter.SetComputePointGradient(false);
247
      filter.SetActiveField(field.GetName(), vtkm::cont::Field::Association::POINTS);
Sujin Philip's avatar
Sujin Philip committed
248 249
      result = filter.Execute(in, policy);
    }
250

251 252 253 254 255 256 257 258 259 260 261 262 263 264 265
    // Remove gradient field from result if it was not requested.
    auto requestedResult = result;
    if (!this->ComputeGradient)
    {
      requestedResult = CopyDataSetStructure(result);
      vtkm::Id numOfFields = static_cast<vtkm::Id>(result.GetNumberOfFields());
      for (vtkm::Id i=0; i < numOfFields; ++i)
      {
        if (result.GetField(i).GetName() != filter.GetOutputFieldName())
        {
          requestedResult.AddField(result.GetField(i));
        }
      }
    }

266
    // convert arrays back to VTK
Sujin Philip's avatar
Sujin Philip committed
267
    if (!fromvtkm::ConvertArrays(result, output))
268 269 270 271
    {
      vtkErrorMacro(<< "Unable to convert VTKm DataSet back to VTK");
      return 0;
    }
272
  }
Sujin Philip's avatar
Sujin Philip committed
273 274 275 276 277 278
  catch (const vtkm::cont::Error& e)
  {
    vtkWarningMacro(<< "VTK-m error: " << e.GetMessage()
                    << "Falling back to serial implementation.");
    return this->Superclass::RequestData(request, inputVector, outputVector);
  }
279 280 281

  return 1;
}