VTK-m Accelerators now support unstructured grids with only voxel elements.

When we encounter a unstructured grid with only voxel elements we will
construct a new VTK-m cell set which uses hexahedron to represent the voxel.
This doesn't resolve the issue of having voxel inside a mixed unstructured
grid.

Accelerators/Vtkm/Testing/Cxx/TestVTKMGradient.cxx

@@ -258,7 +258,7 @@ namespace
     vtkDoubleArray* correctPointArray = vtkArrayDownCast<vtkDoubleArray>(
       vtkDataSet::SafeDownCast(
         correctPointGradients->GetOutput())->GetPointData()->GetArray(resultName));
-    std::cout << "testing point gradients" << std::endl;
+
     if(!IsGradientCorrect(gradPointArray, correctPointArray))
     {
       return EXIT_FAILURE;
@@ -274,7 +274,6 @@ namespace
     pointVorticity->SetComputeDivergence(1);
     pointVorticity->Update();
 
-    std::cout << "testing point vorticity" << std::endl;
     // point stuff
     vtkDoubleArray* vorticityPointArray = vtkArrayDownCast<vtkDoubleArray>(
       vtkDataSet::SafeDownCast(
@@ -284,7 +283,6 @@ namespace
       return EXIT_FAILURE;
     }
 
-    std::cout << "testing point divergence" << std::endl;
     vtkDoubleArray* divergencePointArray = vtkArrayDownCast<vtkDoubleArray>(
       vtkDataSet::SafeDownCast(
         pointVorticity->GetOutput())->GetPointData()->GetArray("Divergence"));
@@ -293,7 +291,6 @@ namespace
       return EXIT_FAILURE;
     }
 
-    std::cout << "testing point QCriterion" << std::endl;
     vtkDoubleArray* qCriterionPointArray = vtkArrayDownCast<vtkDoubleArray>(
       vtkDataSet::SafeDownCast(
         pointVorticity->GetOutput())->GetPointData()->GetArray("Q-criterion"));

Accelerators/Vtkm/Testing/Cxx/TestVTKMGradientAndVorticity.cxx

@@ -323,11 +323,11 @@ namespace
       vtkDataSet::SafeDownCast(
         correctCellGradients->GetOutput())->GetCellData()->GetArray(resultName));
 
-    if(!grid->IsA("vtkUnstructuredGrid"))
+    if(!grid->IsA("vtkStructuredGrid"))
     {
-      // ignore cell gradients if this is an unstructured grid
-      // because the accuracy is so lousy
-      std::cout << "testing cell gradients" << std::endl;
+      // ignore cell gradients on structured grids as the version for
+      // VTK-m differs from VTK. Once VTK-m is able to do stencil based
+      // gradients for point and cells, we can enable this check.
       if(!IsGradientCorrect(gradCellArray, correctCellArray))
       {
         return EXIT_FAILURE;
@@ -341,7 +341,7 @@ namespace
     vtkDoubleArray* correctPointArray = vtkArrayDownCast<vtkDoubleArray>(
       vtkDataSet::SafeDownCast(
         correctPointGradients->GetOutput())->GetPointData()->GetArray(resultName));
-    std::cout << "testing point gradients" << std::endl;
+
     if(!IsGradientCorrect(gradPointArray, correctPointArray))
     {
       return EXIT_FAILURE;
@@ -366,7 +366,6 @@ namespace
     pointVorticity->SetComputeDivergence(1);
     pointVorticity->Update();
 
-    std::cout << "testing cell vorticity" << std::endl;
     // cell stuff
     vtkDoubleArray* vorticityCellArray = vtkArrayDownCast<vtkDoubleArray>(
       vtkDataSet::SafeDownCast(
@@ -376,7 +375,6 @@ namespace
       return EXIT_FAILURE;
     }
 
-    std::cout << "testing point vorticity" << std::endl;
     // point stuff
     vtkDoubleArray* vorticityPointArray = vtkArrayDownCast<vtkDoubleArray>(
       vtkDataSet::SafeDownCast(
@@ -386,7 +384,6 @@ namespace
       return EXIT_FAILURE;
     }
 
-    std::cout << "testing point divergence" << std::endl;
     vtkDoubleArray* divergencePointArray = vtkArrayDownCast<vtkDoubleArray>(
       vtkDataSet::SafeDownCast(
         pointVorticity->GetOutput())->GetPointData()->GetArray("Divergence"));
@@ -395,7 +392,6 @@ namespace
       return EXIT_FAILURE;
     }
 
-    std::cout << "testing point QCriterion" << std::endl;
     vtkDoubleArray* qCriterionPointArray = vtkArrayDownCast<vtkDoubleArray>(
       vtkDataSet::SafeDownCast(
         pointVorticity->GetOutput())->GetPointData()->GetArray("Q-criterion"));

Accelerators/Vtkm/vtkmCellSetSingleType.cxx

@@ -29,11 +29,11 @@
 
 
 namespace {
-class ComputeReverseMapping : public vtkm::worklet::WorkletMapField
+class ComputeSingleTypeReverseMapping : public vtkm::worklet::WorkletMapField
 {
   vtkm::Id NumberOfPointsPerCell;
 public:
-  ComputeReverseMapping(vtkm::Id numberOfPointsPerCell):
+  ComputeSingleTypeReverseMapping(vtkm::Id numberOfPointsPerCell):
     NumberOfPointsPerCell(numberOfPointsPerCell)
   {
 
@@ -140,7 +140,7 @@ typename vtkm::exec::ReverseConnectivityVTK<Device>
     pointIdKey.Allocate(rconnSize);
     this->RConn.Allocate(rconnSize);
 
-    vtkm::worklet::DispatcherMapField<ComputeReverseMapping, Device> dispatcher( ComputeReverseMapping(this->DetermineNumberOfPoints()));
+    vtkm::worklet::DispatcherMapField<ComputeSingleTypeReverseMapping, Device> dispatcher( ComputeSingleTypeReverseMapping(this->DetermineNumberOfPoints()));
     dispatcher.Invoke(vtkm::cont::make_ArrayHandleCounting(0, 1, numberOfCells),
                       this->Connectivity, pointIdKey, this->RConn);
     Algorithm::SortByKey(pointIdKey, this->RConn);

Accelerators/Vtkm/vtkmlib/CellSetConverters.cxx

@@ -25,6 +25,7 @@
 #include "vtkmFilterPolicy.h"
 
 #include <vtkm/cont/serial/DeviceAdapterSerial.h>
+#include <vtkm/cont/tbb/DeviceAdapterTBB.h>
 #include <vtkm/worklet/DispatcherMapTopology.h>
 
 #include "vtkIdTypeArray.h"
@@ -32,6 +33,30 @@
 
 namespace tovtkm {
 
+namespace {
+
+struct ReorderHex : public vtkm::exec::FunctorBase
+{
+  ReorderHex(): Data(nullptr) {}
+  ReorderHex(vtkCellArray* fc): Data(fc->GetPointer()) {}
+
+  void operator()(vtkm::Id index) const
+  {
+    const std::size_t offset = (index * 9) + 1;
+    vtkIdType t = this->Data[offset+3];
+    this->Data[offset+3] = this->Data[offset+2];
+    this->Data[offset+2] = t;
+
+    t = this->Data[offset+7];
+    this->Data[offset+7] = this->Data[offset+6];
+    this->Data[offset+6] = t;
+  }
+
+  vtkIdType* Data;
+};
+
+}
+
 // convert a cell array of a single type to a vtkm CellSetSingleType
 vtkm::cont::DynamicCellSet ConvertSingleType(vtkCellArray* cells, int cellType,
                                              vtkIdType numberOfPoints)
@@ -62,6 +87,30 @@ vtkm::cont::DynamicCellSet ConvertSingleType(vtkCellArray* cells, int cellType,
     c.Fill(numberOfPoints, handle);
     return vtkm::cont::DynamicCellSet(c);
   }
+  case VTK_VOXEL:
+  {
+    //This is encountered when you have an unstructured grid that was
+    //cleaned / thresholded from an image data. At that point VTK should
+    //have converted to hex's as downstream the point coordinates can be
+    //transformed and invalidate the voxel requirements
+
+    //We need to construct a new array that has the correct ordering
+    //divide the array by 4, gets the number of times we need to flip values
+    using Algorithms = typename vtkm::cont::DeviceAdapterAlgorithm<vtkm::cont::DeviceAdapterTagTBB> ;
+
+    //construct through vtkm so that the memory is properly
+    //de-allocated when the DynamicCellSet is destroyed
+    vtkm::cont::ArrayHandle<vtkm::Id, tovtkm::vtkCellArrayContainerTag> fixedCells;
+    fixedCells.Allocate(cells->GetSize());
+    fixedCells.GetStorage().VTKArray()->DeepCopy(cells);
+
+    ReorderHex reorder_hex(fixedCells.GetStorage().VTKArray());
+    Algorithms::Schedule(reorder_hex, cells->GetNumberOfCells());
+
+    CellSetType c(vtkm::CellShapeTagHexahedron(), "cells");
+    c.Fill(numberOfPoints, fixedCells);
+    return vtkm::cont::DynamicCellSet(c);
+  }
   case VTK_QUAD:
   {
     CellSetType c(vtkm::CellShapeTagQuad(), "cells");