From d3512879d1c476fe3d4a4be8c6f05ef18e881e91 Mon Sep 17 00:00:00 2001 From: Robert Maynard Date: Thu, 19 Mar 2020 08:32:43 -0400 Subject: [PATCH 1/5] Remove unneeded template parameters from Gradient. --- vtkm/worklet/Gradient.h | 17 ++++++----------- vtkm/worklet/contour/MarchingCells.h | 13 ++++--------- vtkm/worklet/gradient/CellGradient.h | 4 ---- vtkm/worklet/gradient/PointGradient.h | 4 ---- vtkm/worklet/gradient/StructuredPointGradient.h | 16 ++++++---------- 5 files changed, 16 insertions(+), 38 deletions(-) diff --git a/vtkm/worklet/Gradient.h b/vtkm/worklet/Gradient.h index f786ec99e..d4293cf4d 100644 --- a/vtkm/worklet/Gradient.h +++ b/vtkm/worklet/Gradient.h @@ -50,7 +50,7 @@ struct DeducedPointGrad template void operator()(const CellSetType& cellset) const { - vtkm::worklet::DispatcherMapTopology> dispatcher; + vtkm::worklet::DispatcherMapTopology dispatcher; dispatcher.Invoke(cellset, //topology to iterate on a per point basis cellset, //whole cellset in *this->Points, @@ -60,7 +60,7 @@ struct DeducedPointGrad void operator()(const vtkm::cont::CellSetStructured<3>& cellset) const { - vtkm::worklet::DispatcherPointNeighborhood> dispatcher; + vtkm::worklet::DispatcherPointNeighborhood dispatcher; dispatcher.Invoke(cellset, //topology to iterate on a per point basis *this->Points, *this->Field, @@ -71,7 +71,7 @@ struct DeducedPointGrad void operator()(const vtkm::cont::CellSetPermutation, PermIterType>& cellset) const { - vtkm::worklet::DispatcherPointNeighborhood> dispatcher; + vtkm::worklet::DispatcherPointNeighborhood dispatcher; dispatcher.Invoke(cellset, //topology to iterate on a per point basis *this->Points, *this->Field, @@ -80,7 +80,7 @@ struct DeducedPointGrad void operator()(const vtkm::cont::CellSetStructured<2>& cellset) const { - vtkm::worklet::DispatcherPointNeighborhood> dispatcher; + vtkm::worklet::DispatcherPointNeighborhood dispatcher; dispatcher.Invoke(cellset, //topology to iterate on a per point basis *this->Points, *this->Field, @@ -91,7 +91,7 @@ struct DeducedPointGrad void operator()(const vtkm::cont::CellSetPermutation, PermIterType>& cellset) const { - vtkm::worklet::DispatcherPointNeighborhood> dispatcher; + vtkm::worklet::DispatcherPointNeighborhood dispatcher; dispatcher.Invoke(cellset, //topology to iterate on a per point basis *this->Points, *this->Field, @@ -241,12 +241,7 @@ public: const vtkm::cont::ArrayHandle& field, GradientOutputFields& extraOutput) { - using DispatcherType = - vtkm::worklet::DispatcherMapTopology>; - - vtkm::worklet::gradient::CellGradient worklet; - DispatcherType dispatcher(worklet); - + vtkm::worklet::DispatcherMapTopology dispatcher; dispatcher.Invoke(cells, coords, field, extraOutput); return extraOutput.Gradient; } diff --git a/vtkm/worklet/contour/MarchingCells.h b/vtkm/worklet/contour/MarchingCells.h index 7d3d04b4d..7bd058b49 100644 --- a/vtkm/worklet/contour/MarchingCells.h +++ b/vtkm/worklet/contour/MarchingCells.h @@ -429,8 +429,7 @@ public: const WholeFieldIn& inputField, NormalType& normal) const { - using T = typename WholeFieldIn::ValueType; - vtkm::worklet::gradient::PointGradient gradient; + vtkm::worklet::gradient::PointGradient gradient; gradient(numCells, cellIds, pointId, geometry, pointCoordinates, inputField, normal); } @@ -446,8 +445,6 @@ public: const WholeFieldIn& inputField, NormalType& normal) const { - using T = typename WholeFieldIn::ValueType; - //Optimization for structured cellsets so we can call StructuredPointGradient //and have way faster gradients vtkm::exec::ConnectivityStructured pointGeom(geometry); @@ -459,7 +456,7 @@ public: vtkm::exec::FieldNeighborhood points(pointPortal, boundary); vtkm::exec::FieldNeighborhood field(fieldPortal, boundary); - vtkm::worklet::gradient::StructuredPointGradient gradient; + vtkm::worklet::gradient::StructuredPointGradient gradient; gradient(boundary, points, field, normal); } }; @@ -506,8 +503,7 @@ public: const WholeWeightsIn& weights, NormalType& normal) const { - using T = typename WholeFieldIn::ValueType; - vtkm::worklet::gradient::PointGradient gradient; + vtkm::worklet::gradient::PointGradient gradient; NormalType grad1; gradient(numCells, cellIds, pointId, geometry, pointCoordinates, inputField, grad1); @@ -531,7 +527,6 @@ public: const WholeWeightsIn& weights, NormalType& normal) const { - using T = typename WholeFieldIn::ValueType; //Optimization for structured cellsets so we can call StructuredPointGradient //and have way faster gradients vtkm::exec::ConnectivityStructured pointGeom(geometry); @@ -543,7 +538,7 @@ public: vtkm::exec::FieldNeighborhood points(pointPortal, boundary); vtkm::exec::FieldNeighborhood field(fieldPortal, boundary); - vtkm::worklet::gradient::StructuredPointGradient gradient; + vtkm::worklet::gradient::StructuredPointGradient gradient; NormalType grad1; gradient(boundary, points, field, grad1); diff --git a/vtkm/worklet/gradient/CellGradient.h b/vtkm/worklet/gradient/CellGradient.h index 3668e58f9..dfa3f29d8 100644 --- a/vtkm/worklet/gradient/CellGradient.h +++ b/vtkm/worklet/gradient/CellGradient.h @@ -24,10 +24,6 @@ namespace worklet namespace gradient { -template -using CellGradientInType = vtkm::List; - -template struct CellGradient : vtkm::worklet::WorkletVisitCellsWithPoints { using ControlSignature = void(CellSetIn, diff --git a/vtkm/worklet/gradient/PointGradient.h b/vtkm/worklet/gradient/PointGradient.h index 2f30d89aa..e77322bae 100644 --- a/vtkm/worklet/gradient/PointGradient.h +++ b/vtkm/worklet/gradient/PointGradient.h @@ -25,10 +25,6 @@ namespace worklet namespace gradient { -template -using PointGradientInType = vtkm::List; - -template struct PointGradient : public vtkm::worklet::WorkletVisitPointsWithCells { using ControlSignature = void(CellSetIn, diff --git a/vtkm/worklet/gradient/StructuredPointGradient.h b/vtkm/worklet/gradient/StructuredPointGradient.h index 79964b937..44a921af4 100644 --- a/vtkm/worklet/gradient/StructuredPointGradient.h +++ b/vtkm/worklet/gradient/StructuredPointGradient.h @@ -22,10 +22,6 @@ namespace worklet namespace gradient { -template -using StructuredPointGradientInType = vtkm::List; - -template struct StructuredPointGradient : public vtkm::worklet::WorkletPointNeighborhood { @@ -51,9 +47,9 @@ struct StructuredPointGradient : public vtkm::worklet::WorkletPointNeighborhood vtkm::Vec xi, eta, zeta; this->Jacobian(inputPoints, boundary, xi, eta, zeta); //store the metrics in xi,eta,zeta - T dxi = inputField.Get(1, 0, 0) - inputField.Get(-1, 0, 0); - T deta = inputField.Get(0, 1, 0) - inputField.Get(0, -1, 0); - T dzeta = inputField.Get(0, 0, 1) - inputField.Get(0, 0, -1); + auto dxi = inputField.Get(1, 0, 0) - inputField.Get(-1, 0, 0); + auto deta = inputField.Get(0, 1, 0) - inputField.Get(0, -1, 0); + auto dzeta = inputField.Get(0, 0, 1) - inputField.Get(0, 0, -1); dxi = (boundary.IsRadiusInXBoundary(1) ? dxi * 0.5f : dxi); deta = (boundary.IsRadiusInYBoundary(1) ? deta * 0.5f : deta); @@ -86,9 +82,9 @@ struct StructuredPointGradient : public vtkm::worklet::WorkletPointNeighborhood r[1] = (boundary.IsRadiusInYBoundary(1) ? r[1] * 0.5f : r[1]); r[2] = (boundary.IsRadiusInZBoundary(1) ? r[2] * 0.5f : r[2]); - const T dx = inputField.Get(1, 0, 0) - inputField.Get(-1, 0, 0); - const T dy = inputField.Get(0, 1, 0) - inputField.Get(0, -1, 0); - const T dz = inputField.Get(0, 0, 1) - inputField.Get(0, 0, -1); + auto dx = inputField.Get(1, 0, 0) - inputField.Get(-1, 0, 0); + auto dy = inputField.Get(0, 1, 0) - inputField.Get(0, -1, 0); + auto dz = inputField.Get(0, 0, 1) - inputField.Get(0, 0, -1); outputGradient[0] = static_cast(dx * r[0]); outputGradient[1] = static_cast(dy * r[1]); From 1995e520477572e67a9cc086cc0d4ff72f306837 Mon Sep 17 00:00:00 2001 From: Robert Maynard Date: Thu, 19 Mar 2020 16:22:35 -0400 Subject: [PATCH 2/5] Correct minor issues from refactoring MarchingCells worklets --- vtkm/worklet/contour/CommonState.h | 2 +- vtkm/worklet/contour/MarchingCells.h | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/vtkm/worklet/contour/CommonState.h b/vtkm/worklet/contour/CommonState.h index 25a628ccd..237194633 100644 --- a/vtkm/worklet/contour/CommonState.h +++ b/vtkm/worklet/contour/CommonState.h @@ -29,7 +29,7 @@ struct CommonState } bool MergeDuplicatePoints = true; - bool GenerateNormals = true; + bool GenerateNormals = false; vtkm::cont::ArrayHandle InterpolationWeights; vtkm::cont::ArrayHandle InterpolationEdgeIds; vtkm::cont::ArrayHandle CellIdMap; diff --git a/vtkm/worklet/contour/MarchingCells.h b/vtkm/worklet/contour/MarchingCells.h index 7bd058b49..7e47ea58f 100644 --- a/vtkm/worklet/contour/MarchingCells.h +++ b/vtkm/worklet/contour/MarchingCells.h @@ -611,8 +611,8 @@ vtkm::cont::CellSetSingleType<> execute( const CellSetType& cells, const CoordinateSystem& coordinateSystem, const vtkm::cont::ArrayHandle& inputField, - vtkm::cont::ArrayHandle, StorageTagVertices> vertices, - vtkm::cont::ArrayHandle, StorageTagNormals> normals, + vtkm::cont::ArrayHandle, StorageTagVertices>& vertices, + vtkm::cont::ArrayHandle, StorageTagNormals>& normals, vtkm::worklet::contour::CommonState& sharedState) { using vtkm::worklet::marching_cells::ClassifyCell; From c2bab6c6bbf94851e092ef94682f11a4836cc4c5 Mon Sep 17 00:00:00 2001 From: Robert Maynard Date: Wed, 11 Mar 2020 17:27:54 -0400 Subject: [PATCH 3/5] Introduce FlyingEdges for 3d structured datasets --- vtkm/worklet/Contour.h | 92 ++-- vtkm/worklet/contour/CMakeLists.txt | 9 +- vtkm/worklet/contour/FlyingEdges.h | 221 +++++++++ vtkm/worklet/contour/FlyingEdgesHelpers.h | 198 ++++++++ vtkm/worklet/contour/FlyingEdgesPass1.h | 128 +++++ vtkm/worklet/contour/FlyingEdgesPass2.h | 193 ++++++++ vtkm/worklet/contour/FlyingEdgesPass4.h | 447 ++++++++++++++++++ vtkm/worklet/contour/FlyingEdgesPass5.h | 106 +++++ vtkm/worklet/contour/FlyingEdgesTables.h | 413 ++++++++++++++++ vtkm/worklet/contour/MarchingCells.h | 24 +- .../gradient/StructuredPointGradient.h | 7 + 11 files changed, 1778 insertions(+), 60 deletions(-) create mode 100644 vtkm/worklet/contour/FlyingEdges.h create mode 100644 vtkm/worklet/contour/FlyingEdgesHelpers.h create mode 100644 vtkm/worklet/contour/FlyingEdgesPass1.h create mode 100644 vtkm/worklet/contour/FlyingEdgesPass2.h create mode 100644 vtkm/worklet/contour/FlyingEdgesPass4.h create mode 100644 vtkm/worklet/contour/FlyingEdgesPass5.h create mode 100644 vtkm/worklet/contour/FlyingEdgesTables.h diff --git a/vtkm/worklet/Contour.h b/vtkm/worklet/Contour.h index 6a7c748f3..eb41754f9 100644 --- a/vtkm/worklet/Contour.h +++ b/vtkm/worklet/Contour.h @@ -13,9 +13,11 @@ #include #include +#include #include #include +#include #include @@ -24,6 +26,41 @@ namespace vtkm namespace worklet { + +namespace contour +{ +struct DeduceCoordType +{ + template + void operator()(const CoordinateType& coords, + const CellSetType& cells, + vtkm::cont::CellSetSingleType<>& result, + Args&&... args) const + { + result = marching_cells::execute(cells, coords, std::forward(args)...); + } + + template + void operator()(const vtkm::cont::ArrayHandleUniformPointCoordinates& coords, + const vtkm::cont::CellSetStructured<3>& cells, + vtkm::cont::CellSetSingleType<>& result, + Args&&... args) const + { + result = flying_edges::execute(cells, coords, std::forward(args)...); + } +}; + +struct DeduceCellType +{ + template + void operator()(const CellSetType& cells, CoordinateType&& coordinateSystem, Args&&... args) const + { + vtkm::cont::CastAndCall( + coordinateSystem, contour::DeduceCoordType{}, cells, std::forward(args)...); + } +}; +} + /// \brief Compute the isosurface of a given 3D data set, supports all /// linear cell types class Contour @@ -67,15 +104,15 @@ public: vtkm::cont::CellSetSingleType<> outputCells; vtkm::cont::CastAndCall(cells, - DeduceCellType{}, - this, + contour::DeduceCellType{}, + coordinateSystem, outputCells, isovalues, numIsoValues, - coordinateSystem, input, vertices, - normals); + normals, + this->SharedState); return outputCells; } @@ -100,15 +137,15 @@ public: vtkm::cont::CellSetSingleType<> outputCells; vtkm::cont::CastAndCall(cells, - DeduceCellType{}, - this, + contour::DeduceCellType{}, + coordinateSystem, outputCells, isovalues, numIsoValues, - coordinateSystem, input, vertices, - normals); + normals, + this->SharedState); return outputCells; } @@ -149,45 +186,6 @@ public: void ReleaseCellMapArrays() { this->SharedState.CellIdMap.ReleaseResources(); } private: - struct DeduceCellType - { - template - void operator()(const CellSetType& cells, - ContourAlg* contour, - vtkm::cont::CellSetSingleType<>& result, - Args&&... args) const - { - result = contour->DoRun(cells, std::forward(args)...); - } - }; - - //---------------------------------------------------------------------------- - template - vtkm::cont::CellSetSingleType<> DoRun( - const CellSetType& cells, - const ValueType* isovalues, - const vtkm::Id numIsoValues, - const CoordinateSystem& coordinateSystem, - const vtkm::cont::ArrayHandle& inputField, - vtkm::cont::ArrayHandle, StorageTagVertices> vertices, - vtkm::cont::ArrayHandle, StorageTagNormals> normals) - { - return worklet::marching_cells::execute(isovalues, - numIsoValues, - cells, - coordinateSystem, - inputField, - vertices, - normals, - this->SharedState); - } - vtkm::worklet::contour::CommonState SharedState; }; } diff --git a/vtkm/worklet/contour/CMakeLists.txt b/vtkm/worklet/contour/CMakeLists.txt index 375bb7dea..e386d78f9 100644 --- a/vtkm/worklet/contour/CMakeLists.txt +++ b/vtkm/worklet/contour/CMakeLists.txt @@ -11,8 +11,15 @@ set(headers CommonState.h FieldPropagation.h - MarchingCells.h + FlyingEdges.h + FlyingEdgesHelpers.h + FlyingEdgesPass1.h + FlyingEdgesPass2.h + FlyingEdgesPass4.h + FlyingEdgesPass5.h + FlyingEdgesTables.h MarchingCellTables.h + MarchingCells.h ) #----------------------------------------------------------------------------- diff --git a/vtkm/worklet/contour/FlyingEdges.h b/vtkm/worklet/contour/FlyingEdges.h new file mode 100644 index 000000000..80d4502b9 --- /dev/null +++ b/vtkm/worklet/contour/FlyingEdges.h @@ -0,0 +1,221 @@ + +//============================================================================ +// 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. +//============================================================================ + +#ifndef vtk_m_worklet_contour_flyingedges_h +#define vtk_m_worklet_contour_flyingedges_h + +#include +#include +#include +#include +#include + +#include +#include + +namespace vtkm +{ +namespace worklet +{ +namespace flying_edges +{ + +namespace detail +{ +inline vtkm::cont::CellSetStructured<3> make_metaDataMesh3D(SumXAxis, const vtkm::Id3& pdims) +{ + vtkm::cont::CellSetStructured<3> metaDataMesh; + metaDataMesh.SetPointDimensions(vtkm::Id3{ pdims[1], pdims[2], 1 }); + return metaDataMesh; +} +inline vtkm::cont::CellSetStructured<2> make_metaDataMesh2D(SumXAxis, const vtkm::Id3& pdims) +{ + vtkm::cont::CellSetStructured<2> metaDataMesh; + metaDataMesh.SetPointDimensions(vtkm::Id2{ pdims[1], pdims[2] }); + return metaDataMesh; +} + +inline vtkm::cont::CellSetStructured<3> make_metaDataMesh3D(SumYAxis, const vtkm::Id3& pdims) +{ + vtkm::cont::CellSetStructured<3> metaDataMesh; + metaDataMesh.SetPointDimensions(vtkm::Id3{ pdims[0], pdims[2], 1 }); + return metaDataMesh; +} +inline vtkm::cont::CellSetStructured<2> make_metaDataMesh2D(SumYAxis, const vtkm::Id3& pdims) +{ + vtkm::cont::CellSetStructured<2> metaDataMesh; + metaDataMesh.SetPointDimensions(vtkm::Id2{ pdims[0], pdims[2] }); + return metaDataMesh; +} + +template +vtkm::Id extend_by(vtkm::cont::ArrayHandle& handle, vtkm::Id size) +{ + vtkm::Id oldLen = handle.GetNumberOfValues(); + if (oldLen == 0) + { + handle.Allocate(size); + } + else + { + vtkm::cont::ArrayHandle tempHandle; + tempHandle.Allocate(oldLen + size); + vtkm::cont::Algorithm::CopySubRange(handle, 0, oldLen, tempHandle); + handle = tempHandle; + } + return oldLen; +} +} + + +//---------------------------------------------------------------------------- +template +vtkm::cont::CellSetSingleType<> execute( + const vtkm::cont::CellSetStructured<3>& cells, + const vtkm::cont::ArrayHandleUniformPointCoordinates& coordinateSystem, + const ValueType* isovalues, + const vtkm::Id numIsoValues, + const vtkm::cont::ArrayHandle& inputField, + vtkm::cont::ArrayHandle, StorageTagVertices>& points, + vtkm::cont::ArrayHandle, StorageTagNormals>& normals, + vtkm::worklet::contour::CommonState& sharedState) +{ + //Tasks: + //2. Refactor how we map fields. + // We need the ability unload everything in SharedState after + // we have mapped all fields + //3. Support switching AxisToSum by running this whole thing in a TryExecute + // Passes 5 can ignore this + + using AxisToSum = SumXAxis; + + vtkm::cont::Invoker invoke; + + auto pdims = cells.GetPointDimensions(); + + vtkm::cont::ArrayHandle edgeCases; + edgeCases.Allocate(coordinateSystem.GetNumberOfValues()); + + vtkm::cont::CellSetStructured<3> metaDataMesh3D = detail::make_metaDataMesh3D(AxisToSum{}, pdims); + vtkm::cont::CellSetStructured<2> metaDataMesh2D = detail::make_metaDataMesh2D(AxisToSum{}, pdims); + + vtkm::cont::ArrayHandle metaDataLinearSums; //per point of metaDataMesh + vtkm::cont::ArrayHandle metaDataMin; //per point of metaDataMesh + vtkm::cont::ArrayHandle metaDataMax; //per point of metaDataMesh + vtkm::cont::ArrayHandle metaDataNumTris; //per cell of metaDataMesh + + auto metaDataSums = vtkm::cont::make_ArrayHandleGroupVec<3>(metaDataLinearSums); + + // Since sharedState can be re-used between invocations of contour, + // we need to make sure we reset the size of the Interpolation + // arrays so we don't execute Pass5 over an array that is too large + sharedState.InterpolationEdgeIds.Shrink(0); + sharedState.InterpolationWeights.Shrink(0); + sharedState.CellIdMap.Shrink(0); + + vtkm::cont::ArrayHandle triangle_topology; + for (vtkm::Id i = 0; i < numIsoValues; ++i) + { + ValueType isoval = isovalues[i]; + + //---------------------------------------------------------------------------- + // PASS 1: Process all of the voxel edges that compose each row. Determine the + // edges case classification, count the number of edge intersections, and + // figure out where intersections along the row begins and ends + // (i.e., gather information for computational trimming). + // + // We mark everything as below as it is faster than having the worklet to it + vtkm::cont::Algorithm::Fill(edgeCases, static_cast(FlyingEdges3D::Below)); + ComputePass1 worklet1(isoval, pdims); + invoke(worklet1, metaDataMesh3D, metaDataSums, metaDataMin, metaDataMax, edgeCases, inputField); + + //---------------------------------------------------------------------------- + // PASS 2: Process a single row of voxels/cells. Count the number of other + // axis intersections by topological reasoning from previous edge cases. + // Determine the number of primitives (i.e., triangles) generated from this + // row. Use computational trimming to reduce work. + ComputePass2 worklet2(pdims); + invoke( + worklet2, metaDataMesh2D, metaDataSums, metaDataMin, metaDataMax, metaDataNumTris, edgeCases); + + //---------------------------------------------------------------------------- + // PASS 3: Compute the number of points and triangles that each edge + // row needs to generate by using exclusive scans. + vtkm::cont::Algorithm::ScanExtended(metaDataNumTris, metaDataNumTris); + auto sumTris = + vtkm::cont::ArrayGetValue(metaDataNumTris.GetNumberOfValues() - 1, metaDataNumTris); + if (sumTris > 0) + { + detail::extend_by(triangle_topology, 3 * sumTris); + detail::extend_by(sharedState.CellIdMap, sumTris); + + //By using the current length of points array as the starting value we + //handle multiple contours, by propagating the correct write offset + auto newPointSize = vtkm::cont::Algorithm::ScanExclusive( + metaDataLinearSums, metaDataLinearSums, vtkm::Sum{}, points.GetNumberOfValues()); + detail::extend_by(sharedState.InterpolationEdgeIds, newPointSize); + detail::extend_by(sharedState.InterpolationWeights, newPointSize); + + //---------------------------------------------------------------------------- + // PASS 4: Process voxel rows and generate topology, and interpolation state + ComputePass4 worklet4(isoval, pdims); + invoke(worklet4, + metaDataMesh2D, + metaDataSums, + metaDataMin, + metaDataMax, + metaDataNumTris, + edgeCases, + inputField, + triangle_topology, + sharedState.InterpolationEdgeIds, + sharedState.InterpolationWeights, + sharedState.CellIdMap); + } + + //---------------------------------------------------------------------------- + // PASS 5: Convert the edge interpolation information to point and normals + vtkm::Vec3f origin, spacing; + { //extract out the origin and spacing as these are needed for Pass5 to properly + //interpolate the new points + auto portal = coordinateSystem.ReadPortal(); + origin = portal.GetOrigin(); + spacing = portal.GetSpacing(); + } + if (sharedState.GenerateNormals) + { + normals.Allocate(sharedState.InterpolationEdgeIds.GetNumberOfValues()); + } + + ComputePass5 worklet5(pdims, origin, spacing, sharedState.GenerateNormals); + invoke(worklet5, + sharedState.InterpolationEdgeIds, + sharedState.InterpolationWeights, + points, + inputField, + normals); + } + + vtkm::cont::CellSetSingleType<> outputCells; + outputCells.Fill(points.GetNumberOfValues(), vtkm::CELL_SHAPE_TRIANGLE, 3, triangle_topology); + return outputCells; +} + +} //namespace flying_edges +} +} + +#endif diff --git a/vtkm/worklet/contour/FlyingEdgesHelpers.h b/vtkm/worklet/contour/FlyingEdgesHelpers.h new file mode 100644 index 000000000..a3c164fe6 --- /dev/null +++ b/vtkm/worklet/contour/FlyingEdgesHelpers.h @@ -0,0 +1,198 @@ + +//============================================================================ +// 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. +//============================================================================ + + +#ifndef vtk_m_worklet_contour_flyingedges_helpers_h +#define vtk_m_worklet_contour_flyingedges_helpers_h + +#include + +namespace vtkm +{ +namespace worklet +{ +namespace flying_edges +{ + +struct FlyingEdges3D +{ +public: + // Edge case table values. + enum EdgeClass + { + Below = 0, // below isovalue + Above = 1, // above isovalue + LeftAbove = 1, // left vertex is above isovalue + RightAbove = 2, // right vertex is above isovalue + BothAbove = 3 // entire edge is above isovalue + }; + enum CellClass + { + Interior = 0, + MinBoundary = 1, + MaxBoundary = 2 + }; +}; + +struct SumXAxis +{ + static constexpr vtkm::Id xindex = 0; + static constexpr vtkm::Id yindex = 1; + static constexpr vtkm::Id zindex = 2; +}; +struct SumYAxis +{ + static constexpr vtkm::Id xindex = 1; + static constexpr vtkm::Id yindex = 0; + static constexpr vtkm::Id zindex = 2; +}; + +VTKM_EXEC inline vtkm::Id compute_num_pts(SumXAxis, vtkm::Id nx, vtkm::Id vtkmNotUsed(ny)) +{ + return nx; +} +VTKM_EXEC inline vtkm::Id compute_num_pts(SumYAxis, vtkm::Id vtkmNotUsed(nx), vtkm::Id ny) +{ + return ny; +} + +VTKM_EXEC inline vtkm::Id3 compute_ijk(SumXAxis, const vtkm::Id3& executionSpaceIJK) +{ + return vtkm::Id3{ 0, executionSpaceIJK[0], executionSpaceIJK[1] }; +} +VTKM_EXEC inline vtkm::Id3 compute_ijk(SumYAxis, const vtkm::Id3& executionSpaceIJK) +{ + return vtkm::Id3{ executionSpaceIJK[0], 0, executionSpaceIJK[1] }; +} + + +VTKM_EXEC inline vtkm::Id3 compute_cdims(SumXAxis, + const vtkm::Id3& executionSpacePDims, + vtkm::Id numOfXPoints) +{ + return vtkm::Id3{ numOfXPoints - 1, executionSpacePDims[0] - 1, executionSpacePDims[1] - 1 }; +} +VTKM_EXEC inline vtkm::Id3 compute_cdims(SumYAxis, + const vtkm::Id3& executionSpacePDims, + vtkm::Id numOfYPoints) +{ + return vtkm::Id3{ executionSpacePDims[0] - 1, numOfYPoints - 1, executionSpacePDims[1] - 1 }; +} +VTKM_EXEC inline vtkm::Id3 compute_pdims(SumXAxis, + const vtkm::Id3& executionSpacePDims, + vtkm::Id numOfXPoints) +{ + return vtkm::Id3{ numOfXPoints, executionSpacePDims[0], executionSpacePDims[1] }; +} +VTKM_EXEC inline vtkm::Id3 compute_pdims(SumYAxis, + const vtkm::Id3& executionSpacePDims, + vtkm::Id numOfYPoints) +{ + return vtkm::Id3{ executionSpacePDims[0], numOfYPoints, executionSpacePDims[1] }; +} + +VTKM_EXEC inline vtkm::Id compute_start(SumXAxis, const vtkm::Id3& ijk, const vtkm::Id3& dims) +{ + return (dims[0] * ijk[1]) + ((dims[0] * dims[1]) * ijk[2]); +} +VTKM_EXEC inline vtkm::Id compute_start(SumYAxis, const vtkm::Id3& ijk, const vtkm::Id3& dims) +{ + return ijk[0] + ((dims[0] * dims[1]) * ijk[2]); +} + +VTKM_EXEC inline vtkm::Id4 compute_neighbor_starts(SumXAxis, + const vtkm::Id3& ijk, + const vtkm::Id3& pdims) +{ + //Optimized form of + // return vtkm::Id4 { compute_start(sx, ijk, pdims), + // compute_start(sx, ijk + vtkm::Id3{ 0, 1, 0 }, pdims), + // compute_start(sx, ijk + vtkm::Id3{ 0, 0, 1 }, pdims), + // compute_start(sx, ijk + vtkm::Id3{ 0, 1, 1 }, pdims) }; + const auto sliceSize = (pdims[0] * pdims[1]); + const auto rowPos = (pdims[0] * ijk[1]); + return vtkm::Id4{ rowPos + (sliceSize * ijk[2]), + rowPos + pdims[0] + (sliceSize * ijk[2]), + rowPos + (sliceSize * (ijk[2] + 1)), + rowPos + pdims[0] + (sliceSize * (ijk[2] + 1)) }; +} +VTKM_EXEC inline vtkm::Id4 compute_neighbor_starts(SumYAxis, + const vtkm::Id3& ijk, + const vtkm::Id3& pdims) +{ + //Optimized form of + // return vtkm::Id4{ compute_start(sy, ijk, pdims), + // compute_start(sy, ijk + vtkm::Id3{ 1, 0, 0 }, pdims), + // compute_start(sy, ijk + vtkm::Id3{ 0, 0, 1 }, pdims), + // compute_start(sy, ijk + vtkm::Id3{ 1, 0, 1 }, pdims) }; + const auto sliceSize = (pdims[0] * pdims[1]); + return vtkm::Id4{ ijk[0] + (sliceSize * ijk[2]), + ijk[0] + 1 + (sliceSize * ijk[2]), + ijk[0] + (sliceSize * (ijk[2] + 1)), + ijk[0] + 1 + (sliceSize * (ijk[2] + 1)) }; +} + + + +VTKM_EXEC inline vtkm::Id compute_inc(SumXAxis, const vtkm::Id3&) +{ + return 1; +} +VTKM_EXEC inline vtkm::Id compute_inc(SumYAxis, const vtkm::Id3& dims) +{ + return dims[0]; +} + +//---------------------------------------------------------------------------- +template +VTKM_EXEC inline vtkm::UInt8 getEdgeCase(const WholeEdgeField& edges, + const vtkm::Id4& startPos, + vtkm::Id inc) +{ + vtkm::UInt8 e0 = edges.Get(startPos[0] + inc); + vtkm::UInt8 e1 = edges.Get(startPos[1] + inc); + vtkm::UInt8 e2 = edges.Get(startPos[2] + inc); + vtkm::UInt8 e3 = edges.Get(startPos[3] + inc); + return static_cast(e0 | (e1 << 2) | (e2 << 4) | (e3 << 6)); +} + +//---------------------------------------------------------------------------- +template +VTKM_EXEC inline void adjustTrimBounds(vtkm::Id rightMax, + const WholeEdgeField& edges, + const vtkm::Id4& startPos, + vtkm::Id inc, + vtkm::Id& left, + vtkm::Id& right) +{ + + vtkm::UInt8 e0 = edges.Get(startPos[0] + (left * inc)); + vtkm::UInt8 e1 = edges.Get(startPos[1] + (left * inc)); + vtkm::UInt8 e2 = edges.Get(startPos[2] + (left * inc)); + vtkm::UInt8 e3 = edges.Get(startPos[3] + (left * inc)); + if ((e0 & 0x1) != (e1 & 0x1) || (e1 & 0x1) != (e2 & 0x1) || (e2 & 0x1) != (e3 & 0x1)) + { + left = 0; + } + + e0 = edges.Get(startPos[0] + (right * inc)); + e1 = edges.Get(startPos[1] + (right * inc)); + e2 = edges.Get(startPos[2] + (right * inc)); + e3 = edges.Get(startPos[3] + (right * inc)); + if ((e0 & 0x2) != (e1 & 0x2) || (e1 & 0x2) != (e2 & 0x2) || (e2 & 0x2) != (e3 & 0x2)) + { + right = rightMax; + } +} +} +} +} +#endif diff --git a/vtkm/worklet/contour/FlyingEdgesPass1.h b/vtkm/worklet/contour/FlyingEdgesPass1.h new file mode 100644 index 000000000..66c12525a --- /dev/null +++ b/vtkm/worklet/contour/FlyingEdgesPass1.h @@ -0,0 +1,128 @@ + +//============================================================================ +// 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. +//============================================================================ + + +#ifndef vtk_m_worklet_contour_flyingedges_pass1_h +#define vtk_m_worklet_contour_flyingedges_pass1_h + +#include + +namespace vtkm +{ +namespace worklet +{ +namespace flying_edges +{ + +/* +* Understanding Pass1 in general +* +* PASS 1: Process all of the voxel edges that compose each row. Determine the +* edges case classification, count the number of edge intersections, and +* figure out where intersections along the row begins and ends +* (i.e., gather information for computational trimming). +* +* So in general the algorithm selects a primary axis to stride ( X or Y). +* It does this by forming a plane along the other two axes and marching +* over the sum/primary axis. +* +* So for SumXAxis, this means that we form a YZ plane and march the +* X axis along each point. As we march we are looking at the X axis edge +* that is formed by the current and next point. +* +* So for SumYAxis, this means that we form a XZ plane and march the +* Y axis along each point. As we march we are looking at the Y axis edge +* that is formed by the current and next point. +* +*/ +template +struct ComputePass1 : public vtkm::worklet::WorkletPointNeighborhood +{ + + vtkm::Id NumberOfPoints = 0; + T IsoValue; + + ComputePass1() {} + ComputePass1(T value, const vtkm::Id3& pdims) + : NumberOfPoints(compute_num_pts(AxisToSum{}, pdims[0], pdims[1])) + , IsoValue(value) + { + } + + using ControlSignature = void(CellSetIn, + FieldOut axis_sum, + FieldOut axis_min, + FieldOut axis_max, + WholeArrayInOut edgeData, + WholeArrayIn data); + using ExecutionSignature = void(Boundary, _2, _3, _4, _5, _6); + using InputDomain = _1; + + template + VTKM_EXEC void operator()(const vtkm::exec::BoundaryState& boundary, + vtkm::Id3& axis_sum, + vtkm::Id& axis_min, + vtkm::Id& axis_max, + WholeEdgeField& edges, + const WholeDataField& field) const + { + + const vtkm::Id3 ijk = compute_ijk(AxisToSum{}, boundary.IJK); + const vtkm::Id3 dims = compute_pdims(AxisToSum{}, boundary.PointDimensions, NumberOfPoints); + const vtkm::Id startPos = compute_start(AxisToSum{}, ijk, dims); + const vtkm::Id offset = compute_inc(AxisToSum{}, dims); + + const T value = this->IsoValue; + axis_min = this->NumberOfPoints; + axis_max = 0; + T s1 = field.Get(startPos); + T s0 = s1; + axis_sum = { 0, 0, 0 }; + for (vtkm::Id i = 0; i < NumberOfPoints - 1; ++i) + { + s0 = s1; + s1 = field.Get(startPos + (offset * (i + 1))); + + // We don't explicit write the Below case as that ruins performance. + // It is better to initially fill everything as Below and only + // write the exceptions + vtkm::UInt8 edgeCase = FlyingEdges3D::Below; + if (s0 >= value) + { + edgeCase = FlyingEdges3D::LeftAbove; + } + if (s1 >= value) + { + edgeCase |= FlyingEdges3D::RightAbove; + } + if (edgeCase != FlyingEdges3D::Below) + { + edges.Set(startPos + (offset * i), edgeCase); + } + + if (edgeCase == FlyingEdges3D::LeftAbove || edgeCase == FlyingEdges3D::RightAbove) + { + axis_sum[AxisToSum::xindex] += 1; // increment number of intersections along axis + axis_max = i + 1; + if (axis_min == this->NumberOfPoints) + { + axis_min = i; + } + } + } + } +}; +} +} +} + + +#endif diff --git a/vtkm/worklet/contour/FlyingEdgesPass2.h b/vtkm/worklet/contour/FlyingEdgesPass2.h new file mode 100644 index 000000000..52ca518d2 --- /dev/null +++ b/vtkm/worklet/contour/FlyingEdgesPass2.h @@ -0,0 +1,193 @@ + +//============================================================================ +// 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. +//============================================================================ + + +#ifndef vtk_m_worklet_contour_flyingedges_pass2_h +#define vtk_m_worklet_contour_flyingedges_pass2_h + +#include +#include + +namespace vtkm +{ +namespace worklet +{ +namespace flying_edges +{ + +template +struct ComputePass2 : public vtkm::worklet::WorkletVisitCellsWithPoints +{ + vtkm::Id3 PointDims; + + ComputePass2() {} + explicit ComputePass2(const vtkm::Id3& pdims) + : PointDims(pdims) + { + } + + using ControlSignature = void(CellSetIn, + WholeArrayInOut axis_sums, + FieldInPoint axis_mins, + FieldInPoint axis_maxs, + FieldOutCell cell_tri_count, + WholeArrayIn edgeData); + using ExecutionSignature = void(ThreadIndices, _2, _3, _4, _5, _6); + using InputDomain = _1; + + template + VTKM_EXEC void operator()(const ThreadIndices& threadIndices, + const WholeSumField& axis_sums, + const FieldInPointId& axis_mins, + const FieldInPointId& axis_maxs, + vtkm::Int32& cell_tri_count, + const WholeEdgeField& edges) const + { + // Pass 2. Traverse all cells in the meta data plane. This allows us to + // easily grab the four edge cases bounding this voxel-row + + // find adjusted trim values. + vtkm::Id left = vtkm::Min(axis_mins[0], axis_mins[1]); + left = vtkm::Min(left, axis_mins[2]); + left = vtkm::Min(left, axis_mins[3]); + + vtkm::Id right = vtkm::Max(axis_maxs[0], axis_maxs[1]); + right = vtkm::Max(right, axis_maxs[2]); + right = vtkm::Max(right, axis_maxs[3]); + + const vtkm::Id3 ijk = compute_ijk(AxisToSum{}, threadIndices.GetInputIndex3D()); + const vtkm::Id3 pdims = this->PointDims; + + const vtkm::Id4 startPos = compute_neighbor_starts(AxisToSum{}, ijk, pdims); + const vtkm::Id axis_inc = compute_inc(AxisToSum{}, pdims); + + vtkm::Vec onBoundary(false, false, false); //updated in for-loop + onBoundary[AxisToSum::yindex] = (ijk[AxisToSum::yindex] >= (pdims[AxisToSum::yindex] - 2)); + onBoundary[AxisToSum::zindex] = (ijk[AxisToSum::zindex] >= (pdims[AxisToSum::zindex] - 2)); + + cell_tri_count = 0; + vtkm::Id3 sums = axis_sums.Get(threadIndices.GetIndicesIncident()[0]); + vtkm::Id3 adj_row_sum(0, 0, 0); + vtkm::Id3 adj_col_sum(0, 0, 0); + if (onBoundary[AxisToSum::yindex]) + { + adj_row_sum = axis_sums.Get(threadIndices.GetIndicesIncident()[1]); + } + if (onBoundary[AxisToSum::zindex]) + { + adj_col_sum = axis_sums.Get(threadIndices.GetIndicesIncident()[3]); + } + + if (left == pdims[AxisToSum::xindex] && right == 0) + { + //verify that we have nothing to generate and early terminate. + bool mins_same = (axis_mins[0] == axis_mins[1] && axis_mins[0] == axis_mins[2] && + axis_mins[0] == axis_mins[3]); + bool maxs_same = (axis_maxs[0] == axis_maxs[1] && axis_maxs[0] == axis_maxs[2] && + axis_maxs[0] == axis_maxs[3]); + if (mins_same && maxs_same) + { + return; + } + else + { + left = 0; + right = pdims[AxisToSum::xindex] - 1; + } + } + + // The trim edges may need adjustment if the contour travels between rows + // of edges (without intersecting these edges). This means checking + // whether the trim faces at (left,rightR) made up of the edges intersect + // the contour. Basically just an intersection operation. + adjustTrimBounds(pdims[AxisToSum::xindex] - 1, edges, startPos, axis_inc, left, right); + + for (vtkm::Id i = left; i < right; ++i) // run along the trimmed voxels + { + vtkm::UInt8 edgeCase = getEdgeCase(edges, startPos, (axis_inc * i)); + vtkm::UInt8 numTris = data::GetNumberOfPrimitives(edgeCase); + if (numTris > 0) + { + cell_tri_count += numTris; + + // Count the number of y- and z-points to be generated. Pass# 1 counted + // the number of x-intersections along the x-edges. Now we count all + // intersections on the y- and z-voxel axes. + auto* edgeUses = data::GetEdgeUses(edgeCase); + + onBoundary[AxisToSum::xindex] = (i >= (pdims[AxisToSum::xindex] - 2)); + + // row axes edge always counted + sums[AxisToSum::yindex] += edgeUses[4]; + // col axes edge always counted + sums[AxisToSum::zindex] += edgeUses[8]; + + // handle boundary + this->CountBoundaryEdgeUses(onBoundary, edgeUses, sums, adj_row_sum, adj_col_sum); + } + } + + axis_sums.Set(threadIndices.GetIndicesIncident()[0], sums); + if (onBoundary[AxisToSum::yindex]) + { + axis_sums.Set(threadIndices.GetIndicesIncident()[1], adj_row_sum); + } + if (onBoundary[AxisToSum::zindex]) + { + axis_sums.Set(threadIndices.GetIndicesIncident()[3], adj_col_sum); + } + } + + //---------------------------------------------------------------------------- + // Count intersections along voxel axes. When traversing the volume across + // edges, the voxel axes on the boundary may be undefined near boundaries + // (because there are no fully-formed cells). Thus the voxel axes on the + // boundary are treated specially. + // + // Only on these boundaries do we write to the metaData of our neighbor + // as it is safe as those + VTKM_EXEC inline void CountBoundaryEdgeUses(vtkm::Vec onBoundary, + vtkm::UInt8 const* const edgeUses, + vtkm::Id3& sums, + vtkm::Id3& adj_row_sum, + vtkm::Id3& adj_col_sum) const + { + if (onBoundary[AxisToSum::xindex]) //+x boundary + { + sums[AxisToSum::yindex] += edgeUses[5]; + sums[AxisToSum::zindex] += edgeUses[9]; + if (onBoundary[AxisToSum::yindex]) //+x +y + { + adj_row_sum[AxisToSum::zindex] += edgeUses[11]; + } + if (onBoundary[AxisToSum::zindex]) //+x +z + { + adj_col_sum[AxisToSum::yindex] += edgeUses[7]; + } + } + if (onBoundary[AxisToSum::yindex]) //+y boundary + { + adj_row_sum[AxisToSum::zindex] += edgeUses[10]; + } + if (onBoundary[AxisToSum::zindex]) //+z boundary + { + adj_col_sum[AxisToSum::yindex] += edgeUses[6]; + } + } +}; +} +} +} + +#endif diff --git a/vtkm/worklet/contour/FlyingEdgesPass4.h b/vtkm/worklet/contour/FlyingEdgesPass4.h new file mode 100644 index 000000000..b5ce44f20 --- /dev/null +++ b/vtkm/worklet/contour/FlyingEdgesPass4.h @@ -0,0 +1,447 @@ + +//============================================================================ +// 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. +//============================================================================ + + +#ifndef vtk_m_worklet_contour_flyingedges_pass4_h +#define vtk_m_worklet_contour_flyingedges_pass4_h + + +#include +#include + +namespace vtkm +{ +namespace worklet +{ +namespace flying_edges +{ + +VTKM_EXEC inline vtkm::Id3 compute_incs3d(const vtkm::Id3& dims) +{ + return vtkm::Id3{ 1, dims[0], (dims[0] * dims[1]) }; +} + +template +VTKM_EXEC inline void interpolate_weight(T value, + const vtkm::Id2& iEdge, + vtkm::Id writeIndex, + const WholeField& field, + EdgeField& interpolatedEdgeIds, + WeightField& weights) +{ + interpolatedEdgeIds.Set(writeIndex, iEdge); + T s0 = field.Get(iEdge[0]); + T s1 = field.Get(iEdge[1]); + + auto t = (value - s0) / (s1 - s0); + weights.Set(writeIndex, static_cast(t)); +} + +VTKM_EXEC inline bool case_includes_axes(vtkm::UInt8 const* const edgeUses) +{ + return (edgeUses[0] != 0 || edgeUses[4] != 0 || edgeUses[8] != 0); +} + +template +VTKM_EXEC inline void generate_tris(vtkm::Id inputCellId, + vtkm::UInt8 edgeCase, + vtkm::UInt8 numTris, + vtkm::Id* edgeIds, + vtkm::Int32& triId, + const WholeConnField& conn, + const WholeCellIdField& cellIds) +{ + auto* edges = data::GetTriEdgeCases(edgeCase); + vtkm::Id edgeIndex = 1; + vtkm::Id index = static_cast(triId) * 3; + for (vtkm::UInt8 i = 0; i < numTris; ++i) + { + cellIds.Set(triId + i, inputCellId); + + //We use edgeIndex, edgeIndex+2, edgeIndex+1 to keep + //the same winding for the triangles that marching celss + //produced. By keeping the winding the same we make sure + //that 'fast' normals are consistent with the marching + //cells version + conn.Set(index, edgeIds[edges[edgeIndex]]); + conn.Set(index + 1, edgeIds[edges[edgeIndex + 2]]); + conn.Set(index + 2, edgeIds[edges[edgeIndex + 1]]); + index += 3; + edgeIndex += 3; + } + triId += numTris; +} + +// Helper function to set up the point ids on voxel edges. +//---------------------------------------------------------------------------- +template +VTKM_EXEC inline void init_voxelIds(AxisToSum, + vtkm::UInt8 edgeCase, + const FieldInPointId3& axis_sums, + vtkm::Id* edgeIds) +{ + auto* edgeUses = data::GetEdgeUses(edgeCase); + edgeIds[0] = axis_sums[0][AxisToSum::xindex]; // x-edges + edgeIds[1] = axis_sums[1][AxisToSum::xindex]; + edgeIds[2] = axis_sums[3][AxisToSum::xindex]; + edgeIds[3] = axis_sums[2][AxisToSum::xindex]; + edgeIds[4] = axis_sums[0][AxisToSum::yindex]; // y-edges + edgeIds[5] = edgeIds[4] + edgeUses[4]; + edgeIds[6] = axis_sums[3][AxisToSum::yindex]; + edgeIds[7] = edgeIds[6] + edgeUses[6]; + edgeIds[8] = axis_sums[0][AxisToSum::zindex]; // z-edges + edgeIds[9] = edgeIds[8] + edgeUses[8]; + edgeIds[10] = axis_sums[1][AxisToSum::zindex]; + edgeIds[11] = edgeIds[10] + edgeUses[10]; +} + +// Helper function to advance the point ids along voxel rows. +//---------------------------------------------------------------------------- +VTKM_EXEC inline void advance_voxelIds(vtkm::UInt8 const* const edgeUses, vtkm::Id* edgeIds) +{ + edgeIds[0] += edgeUses[0]; // x-edges + edgeIds[1] += edgeUses[1]; + edgeIds[2] += edgeUses[2]; + edgeIds[3] += edgeUses[3]; + edgeIds[4] += edgeUses[4]; // y-edges + edgeIds[5] = edgeIds[4] + edgeUses[5]; + edgeIds[6] += edgeUses[6]; + edgeIds[7] = edgeIds[6] + edgeUses[7]; + edgeIds[8] += edgeUses[8]; // z-edges + edgeIds[9] = edgeIds[8] + edgeUses[9]; + edgeIds[10] += edgeUses[10]; + edgeIds[11] = edgeIds[10] + edgeUses[11]; +} + + + +template +struct ComputePass4 : public vtkm::worklet::WorkletVisitCellsWithPoints +{ + + vtkm::Id3 PointDims; + T IsoValue; + + ComputePass4() {} + ComputePass4(T value, const vtkm::Id3& pdims) + : PointDims(pdims) + , IsoValue(value) + { + } + + using ControlSignature = void(CellSetIn, + FieldInPoint axis_sums, + FieldInPoint axis_mins, + FieldInPoint axis_maxs, + WholeArrayIn cell_tri_count, + WholeArrayIn edgeData, + WholeArrayIn data, + WholeArrayOut connectivity, + WholeArrayOut edgeIds, + WholeArrayOut weights, + WholeArrayOut inputCellIds); + using ExecutionSignature = + void(ThreadIndices, _2, _3, _4, _5, _6, _7, _8, _9, _10, _11, WorkIndex); + + template + VTKM_EXEC void operator()(const ThreadIndices& threadIndices, + const FieldInPointId3& axis_sums, + const FieldInPointId& axis_mins, + const FieldInPointId& axis_maxs, + const WholeTriField& cellTriCount, + const WholeEdgeField& edges, + const WholeDataField& field, + const WholeConnField& conn, + const WholeEdgeIdField& interpolatedEdgeIds, + const WholeWeightField& weights, + const WholeCellIdField& inputCellIds, + vtkm::Id oidx) const + { + //This works as cellTriCount was computed with ScanExtended + //and therefore has one more entry than the number of cells + vtkm::Int32 cell_tri_offset = cellTriCount.Get(oidx); + vtkm::Int32 next_tri_offset = cellTriCount.Get(oidx + 1); + if (cell_tri_offset == next_tri_offset) + { //we produce nothing + return; + } + + // find adjusted trim values. + vtkm::Id left = vtkm::Min(axis_mins[0], axis_mins[1]); + left = vtkm::Min(left, axis_mins[2]); + left = vtkm::Min(left, axis_mins[3]); + + vtkm::Id right = vtkm::Max(axis_maxs[0], axis_maxs[1]); + right = vtkm::Max(right, axis_maxs[2]); + right = vtkm::Max(right, axis_maxs[3]); + + vtkm::Id3 ijk = compute_ijk(AxisToSum{}, threadIndices.GetInputIndex3D()); + const vtkm::Id3 pdims = this->PointDims; + const vtkm::Id4 startPos = compute_neighbor_starts(AxisToSum{}, ijk, pdims); + const vtkm::Id axis_inc = compute_inc(AxisToSum{}, pdims); + + if (left == pdims[AxisToSum::xindex] && right == 0) + { + //verify that we have nothing to generate and early terminate. + bool mins_same = (axis_mins[0] == axis_mins[1] && axis_mins[0] == axis_mins[2] && + axis_mins[0] == axis_mins[3]); + bool maxs_same = (axis_maxs[0] == axis_maxs[1] && axis_maxs[0] == axis_maxs[2] && + axis_maxs[0] == axis_maxs[3]); + if (mins_same && maxs_same) + { + return; + } + else + { + left = 0; + right = pdims[AxisToSum::xindex] - 1; + } + } + + // The trim edges may need adjustment if the contour travels between rows + // of edges (without intersecting these edges). This means checking + // whether the trim faces at (left,right) made up of the edges intersect + // the contour. + adjustTrimBounds(pdims[AxisToSum::xindex] - 1, edges, startPos, axis_inc, left, right); + if (left == right) + { + return; + } + + const vtkm::UInt8 yLoc = + (ijk[AxisToSum::yindex] < 1 + ? FlyingEdges3D::MinBoundary + : (ijk[AxisToSum::yindex] >= (pdims[AxisToSum::yindex] - 2) ? FlyingEdges3D::MaxBoundary + : FlyingEdges3D::Interior)); + const vtkm::UInt8 zLoc = + (ijk[AxisToSum::zindex] < 1 + ? FlyingEdges3D::MinBoundary + : (ijk[AxisToSum::zindex] >= (pdims[AxisToSum::zindex] - 2) ? FlyingEdges3D::MaxBoundary + : FlyingEdges3D::Interior)); + const vtkm::UInt8 yzLoc = static_cast((yLoc << 2) | (zLoc << 4)); + + + + const vtkm::Id3 increments = compute_incs3d(pdims); + vtkm::Id edgeIds[12]; + auto edgeCase = getEdgeCase(edges, startPos, (axis_inc * left)); + init_voxelIds(AxisToSum{}, edgeCase, axis_sums, edgeIds); + for (vtkm::Id i = left; i < right; ++i) // run along the trimmed voxels + { + ijk[AxisToSum::xindex] = i; + edgeCase = getEdgeCase(edges, startPos, (axis_inc * i)); + vtkm::UInt8 numTris = data::GetNumberOfPrimitives(edgeCase); + if (numTris > 0) + { + //compute what the current cellId is + vtkm::Id cellId = compute_start(AxisToSum{}, ijk, pdims - vtkm::Id3{ 1, 1, 1 }); + + // Start by generating triangles for this case + generate_tris(cellId, edgeCase, numTris, edgeIds, cell_tri_offset, conn, inputCellIds); + + // Now generate edgeIds and weights along voxel axes if needed. Remember to take + // boundary into account. + vtkm::UInt8 loc = static_cast( + yzLoc | (i < 1 ? FlyingEdges3D::MinBoundary + : (i >= (pdims[AxisToSum::xindex] - 2) ? FlyingEdges3D::MaxBoundary + : FlyingEdges3D::Interior))); + auto* edgeUses = data::GetEdgeUses(edgeCase); + if (loc != FlyingEdges3D::Interior || case_includes_axes(edgeUses)) + { + this->GenerateWeights(loc, + field, + interpolatedEdgeIds, + weights, + startPos, + increments, + (axis_inc * i), + edgeUses, + edgeIds); + } + advance_voxelIds(edgeUses, edgeIds); + } + } + } + + //---------------------------------------------------------------------------- + template + VTKM_EXEC inline void GenerateWeights(vtkm::UInt8 loc, + const WholeDataField& field, + const WholeIEdgeField& interpolatedEdgeIds, + const WholeWeightField& weights, + const vtkm::Id4& startPos, + const vtkm::Id3& incs, + vtkm::Id offset, + vtkm::UInt8 const* const edgeUses, + vtkm::Id* edgeIds) const + { + vtkm::Id2 pos(startPos[0] + offset, 0); + //EdgesUses 0,4,8 work for Y axis + if (edgeUses[0]) + { // edgesUses[0] == i axes edge + pos[1] = startPos[0] + offset + incs[AxisToSum::xindex]; + interpolate_weight(this->IsoValue, pos, edgeIds[0], field, interpolatedEdgeIds, weights); + } + if (edgeUses[4]) + { // edgesUses[4] == j axes edge + pos[1] = startPos[1] + offset; + interpolate_weight(this->IsoValue, pos, edgeIds[4], field, interpolatedEdgeIds, weights); + } + if (edgeUses[8]) + { // edgesUses[8] == k axes edge + pos[1] = startPos[2] + offset; + interpolate_weight(this->IsoValue, pos, edgeIds[8], field, interpolatedEdgeIds, weights); + } + // On the boundary cells special work has to be done to cover the partial + // cell axes. These are boundary situations where the voxel axes is not + // fully formed. These situations occur on the +x,+y,+z volume + // boundaries. (The other cases fall through the default: case which is + // expected.) + // + // Note that loc is one of 27 regions in the volume, with (0,1,2) + // indicating (interior, min, max) along coordinate axes. + switch (loc) + { + case 2: + case 6: + case 18: + case 22: //+x + this->InterpolateEdge( + pos[0], incs, 5, edgeUses, edgeIds, field, interpolatedEdgeIds, weights); + this->InterpolateEdge( + pos[0], incs, 9, edgeUses, edgeIds, field, interpolatedEdgeIds, weights); + break; + case 8: + case 9: + case 24: + case 25: //+y + this->InterpolateEdge( + pos[0], incs, 1, edgeUses, edgeIds, field, interpolatedEdgeIds, weights); + this->InterpolateEdge( + pos[0], incs, 10, edgeUses, edgeIds, field, interpolatedEdgeIds, weights); + break; + case 32: + case 33: + case 36: + case 37: //+z + this->InterpolateEdge( + pos[0], incs, 2, edgeUses, edgeIds, field, interpolatedEdgeIds, weights); + this->InterpolateEdge( + pos[0], incs, 6, edgeUses, edgeIds, field, interpolatedEdgeIds, weights); + break; + case 10: + case 26: //+x +y + this->InterpolateEdge( + pos[0], incs, 1, edgeUses, edgeIds, field, interpolatedEdgeIds, weights); + this->InterpolateEdge( + pos[0], incs, 5, edgeUses, edgeIds, field, interpolatedEdgeIds, weights); + this->InterpolateEdge( + pos[0], incs, 9, edgeUses, edgeIds, field, interpolatedEdgeIds, weights); + this->InterpolateEdge( + pos[0], incs, 10, edgeUses, edgeIds, field, interpolatedEdgeIds, weights); + this->InterpolateEdge( + pos[0], incs, 11, edgeUses, edgeIds, field, interpolatedEdgeIds, weights); + break; + case 34: + case 38: //+x +z + this->InterpolateEdge( + pos[0], incs, 2, edgeUses, edgeIds, field, interpolatedEdgeIds, weights); + this->InterpolateEdge( + pos[0], incs, 5, edgeUses, edgeIds, field, interpolatedEdgeIds, weights); + this->InterpolateEdge( + pos[0], incs, 9, edgeUses, edgeIds, field, interpolatedEdgeIds, weights); + this->InterpolateEdge( + pos[0], incs, 6, edgeUses, edgeIds, field, interpolatedEdgeIds, weights); + this->InterpolateEdge( + pos[0], incs, 7, edgeUses, edgeIds, field, interpolatedEdgeIds, weights); + break; + case 40: + case 41: //+y +z + this->InterpolateEdge( + pos[0], incs, 1, edgeUses, edgeIds, field, interpolatedEdgeIds, weights); + this->InterpolateEdge( + pos[0], incs, 2, edgeUses, edgeIds, field, interpolatedEdgeIds, weights); + this->InterpolateEdge( + pos[0], incs, 3, edgeUses, edgeIds, field, interpolatedEdgeIds, weights); + this->InterpolateEdge( + pos[0], incs, 6, edgeUses, edgeIds, field, interpolatedEdgeIds, weights); + this->InterpolateEdge( + pos[0], incs, 10, edgeUses, edgeIds, field, interpolatedEdgeIds, weights); + break; + case 42: //+x +y +z happens no more than once per volume + this->InterpolateEdge( + pos[0], incs, 1, edgeUses, edgeIds, field, interpolatedEdgeIds, weights); + this->InterpolateEdge( + pos[0], incs, 2, edgeUses, edgeIds, field, interpolatedEdgeIds, weights); + this->InterpolateEdge( + pos[0], incs, 3, edgeUses, edgeIds, field, interpolatedEdgeIds, weights); + this->InterpolateEdge( + pos[0], incs, 5, edgeUses, edgeIds, field, interpolatedEdgeIds, weights); + this->InterpolateEdge( + pos[0], incs, 9, edgeUses, edgeIds, field, interpolatedEdgeIds, weights); + this->InterpolateEdge( + pos[0], incs, 10, edgeUses, edgeIds, field, interpolatedEdgeIds, weights); + this->InterpolateEdge( + pos[0], incs, 11, edgeUses, edgeIds, field, interpolatedEdgeIds, weights); + this->InterpolateEdge( + pos[0], incs, 6, edgeUses, edgeIds, field, interpolatedEdgeIds, weights); + this->InterpolateEdge( + pos[0], incs, 7, edgeUses, edgeIds, field, interpolatedEdgeIds, weights); + break; + default: // interior, or -x,-y,-z boundaries + return; + } + } + + // Indicate whether voxel axes need processing for this case. + //---------------------------------------------------------------------------- + template + VTKM_EXEC inline void InterpolateEdge(vtkm::Id currentIdx, + const vtkm::Id3& incs, + vtkm::Id edgeNum, + vtkm::UInt8 const* const edgeUses, + vtkm::Id* edgeIds, + const WholeField& field, + const WholeIEdgeField& interpolatedEdgeIds, + const WholeWeightField& weights) const + { + // if this edge is not used then get out + if (!edgeUses[edgeNum]) + { + return; + } + const vtkm::Id writeIndex = edgeIds[edgeNum]; + vtkm::Id2 iEdge; + + // build the edge information + vtkm::Vec verts = data::GetVertMap(edgeNum); + + vtkm::Id3 offsets = data::GetVertOffsets(AxisToSum{}, verts[0]); + iEdge[0] = currentIdx + vtkm::Dot(offsets, incs); + offsets = data::GetVertOffsets(AxisToSum{}, verts[1]); + iEdge[1] = currentIdx + vtkm::Dot(offsets, incs); + + interpolate_weight(this->IsoValue, iEdge, writeIndex, field, interpolatedEdgeIds, weights); + } +}; +} +} +} +#endif diff --git a/vtkm/worklet/contour/FlyingEdgesPass5.h b/vtkm/worklet/contour/FlyingEdgesPass5.h new file mode 100644 index 000000000..7873ad601 --- /dev/null +++ b/vtkm/worklet/contour/FlyingEdgesPass5.h @@ -0,0 +1,106 @@ + +//============================================================================ +// 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. +//============================================================================ + + +#ifndef vtk_m_worklet_contour_flyingedges_pass5_h +#define vtk_m_worklet_contour_flyingedges_pass5_h + + +#include +#include + +#include +#include + +namespace vtkm +{ +namespace worklet +{ +namespace flying_edges +{ + +template +struct ComputePass5 : public vtkm::worklet::WorkletMapField +{ + + vtkm::internal::ArrayPortalUniformPointCoordinates Coordinates; + bool GenerateNormals; + + ComputePass5() {} + ComputePass5(const vtkm::Id3& pdims, + const vtkm::Vec3f& origin, + const vtkm::Vec3f& spacing, + bool generateNormals) + : Coordinates(pdims, origin, spacing) + , GenerateNormals(generateNormals) + { + } + + using ControlSignature = void(FieldIn interpEdgeIds, + FieldIn interpWeight, + FieldOut points, + WholeArrayIn field, + WholeArrayOut normals); + using ExecutionSignature = void(_1, _2, _3, _4, _5, WorkIndex); + + template + VTKM_EXEC void operator()(const vtkm::Id2& interpEdgeIds, + vtkm::FloatDefault weight, + vtkm::Vec& outPoint, + const WholeInputField& field, + WholeNormalField& normals, + vtkm::Id oidx) const + { + { + vtkm::Vec3f point1 = this->Coordinates.Get(interpEdgeIds[0]); + vtkm::Vec3f point2 = this->Coordinates.Get(interpEdgeIds[1]); + outPoint = vtkm::Lerp(point1, point2, weight); + } + + if (this->GenerateNormals) + { + vtkm::Vec g0, g1; + const vtkm::Id3& dims = this->Coordinates.GetDimensions(); + vtkm::Id3 ijk{ interpEdgeIds[0] % dims[0], + (interpEdgeIds[0] / dims[0]) % dims[1], + interpEdgeIds[0] / (dims[0] * dims[1]) }; + + vtkm::worklet::gradient::StructuredPointGradient gradient; + vtkm::exec::BoundaryState boundary(ijk, dims); + vtkm::exec::FieldNeighborhood + coord_neighborhood(this->Coordinates, boundary); + + vtkm::exec::FieldNeighborhood field_neighborhood(field, boundary); + + + //compute the gradient at point 1 + gradient(boundary, coord_neighborhood, field_neighborhood, g0); + + //compute the gradient at point 2. This optimization can be optimized + boundary.IJK = vtkm::Id3{ interpEdgeIds[1] % dims[0], + (interpEdgeIds[1] / dims[0]) % dims[1], + interpEdgeIds[1] / (dims[0] * dims[1]) }; + gradient(boundary, coord_neighborhood, field_neighborhood, g1); + + vtkm::Vec3f n = vtkm::Lerp(g0, g1, weight); + const auto mag2 = vtkm::MagnitudeSquared(n); + if (mag2 > 0.) + { + n = n * vtkm::RSqrt(mag2); + } + normals.Set(oidx, n); + } + } +}; +} +} +} +#endif diff --git a/vtkm/worklet/contour/FlyingEdgesTables.h b/vtkm/worklet/contour/FlyingEdgesTables.h new file mode 100644 index 000000000..7a34dd473 --- /dev/null +++ b/vtkm/worklet/contour/FlyingEdgesTables.h @@ -0,0 +1,413 @@ + +//============================================================================ +// 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. +//============================================================================ + +#ifndef vtk_m_worklet_contour_flyingedges_tables_h +#define vtk_m_worklet_contour_flyingedges_tables_h + +#include + +namespace vtkm +{ +namespace worklet +{ +namespace flying_edges +{ +namespace data +{ + +VTKM_EXEC inline vtkm::UInt8 GetNumberOfPrimitives(vtkm::UInt8 edgecase) +{ + VTKM_STATIC_CONSTEXPR_ARRAY vtkm::UInt8 numTris[256] = { + 0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 2, 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 3, + 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 3, 2, 3, 3, 2, 3, 4, 4, 3, 3, 4, 4, 3, 4, 5, 5, 2, + 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 3, 2, 3, 3, 4, 3, 2, 4, 3, 3, 4, 4, 5, 4, 3, 5, 2, + 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 4, 3, 4, 4, 3, 4, 3, 5, 2, 4, 5, 5, 4, 5, 4, 2, 1, + 1, 2, 2, 3, 2, 3, 3, 4, 2, 3, 3, 4, 3, 4, 4, 3, 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 4, 5, 5, 4, + 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 2, 3, 4, 5, 3, 2, 3, 4, 4, 3, 4, 5, 5, 4, 4, 5, 3, 2, 5, 2, 4, 1, + 2, 3, 3, 4, 3, 4, 4, 5, 3, 4, 4, 5, 2, 3, 3, 2, 3, 4, 4, 5, 4, 3, 5, 4, 4, 5, 5, 2, 3, 2, 4, 1, + 3, 4, 4, 5, 4, 5, 5, 2, 4, 5, 3, 4, 3, 4, 2, 1, 2, 3, 3, 2, 3, 2, 4, 1, 3, 4, 2, 1, 2, 1, 1, 0 + }; + + return numTris[edgecase]; +} + +VTKM_EXEC inline vtkm::UInt8 const* GetEdgeUses(vtkm::UInt8 edgecase) +{ + VTKM_STATIC_CONSTEXPR_ARRAY vtkm::UInt8 edgeUses[128][12] = { + // This is [128][12] as idx 0 == idx 254, idx 1 == 253... + // + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, { 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0 }, + { 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0 }, { 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0 }, + { 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0 }, { 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0 }, + { 1, 1, 0, 0, 1, 1, 0, 0, 0, 1, 1, 0 }, { 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 1, 0 }, + { 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1 }, { 1, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 1 }, + { 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1 }, { 0, 1, 0, 0, 1, 0, 0, 0, 1, 1, 0, 1 }, + { 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1 }, { 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 1 }, + { 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 1 }, { 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1 }, + { 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0 }, { 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0 }, + { 1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 0 }, { 0, 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 0 }, + { 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0 }, { 1, 1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0 }, + { 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0 }, { 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0 }, + { 0, 1, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1 }, { 1, 1, 1, 0, 1, 1, 1, 0, 0, 0, 0, 1 }, + { 1, 1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 1 }, { 0, 1, 1, 0, 1, 0, 1, 0, 0, 1, 0, 1 }, + { 0, 0, 1, 0, 1, 1, 1, 0, 1, 0, 1, 1 }, { 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 1 }, + { 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1 }, { 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 1 }, + { 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0 }, { 1, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0, 0 }, + { 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0 }, { 0, 0, 1, 0, 1, 1, 0, 1, 1, 0, 0, 0 }, + { 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0 }, { 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0 }, + { 1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0 }, { 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0 }, + { 0, 1, 1, 0, 0, 1, 0, 1, 0, 1, 0, 1 }, { 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1 }, + { 1, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1 }, { 0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 1 }, + { 0, 0, 1, 0, 1, 1, 0, 1, 0, 1, 1, 1 }, { 1, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 1 }, + { 1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 1 }, { 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 1, 1 }, + { 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0 }, { 1, 0, 0, 0, 1, 0, 1, 1, 0, 1, 0, 0 }, + { 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0 }, { 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0 }, + { 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0 }, { 1, 1, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0 }, + { 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 1, 0 }, { 0, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0 }, + { 0, 1, 0, 0, 0, 1, 1, 1, 1, 1, 0, 1 }, { 1, 1, 0, 0, 1, 1, 1, 1, 0, 1, 0, 1 }, + { 1, 1, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1 }, { 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1 }, + { 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1 }, { 1, 0, 0, 0, 0, 1, 1, 1, 0, 1, 1, 1 }, + { 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1 }, { 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1 }, + { 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0 }, { 1, 0, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0 }, + { 1, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0 }, { 0, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0 }, + { 0, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0 }, { 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0 }, + { 1, 1, 0, 1, 1, 1, 1, 0, 0, 1, 0, 0 }, { 0, 1, 0, 1, 0, 1, 1, 0, 1, 1, 0, 0 }, + { 0, 1, 0, 1, 0, 1, 1, 0, 0, 0, 1, 1 }, { 1, 1, 0, 1, 1, 1, 1, 0, 1, 0, 1, 1 }, + { 1, 1, 0, 1, 0, 0, 1, 0, 0, 1, 1, 1 }, { 0, 1, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1 }, + { 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1 }, { 1, 0, 0, 1, 0, 1, 1, 0, 1, 0, 0, 1 }, + { 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1 }, { 0, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1 }, + { 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0 }, { 1, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 0 }, + { 1, 0, 1, 1, 0, 1, 0, 0, 1, 1, 1, 0 }, { 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0 }, + { 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 0, 0 }, { 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 0, 0 }, { 0, 1, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0 }, + { 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1 }, { 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1 }, + { 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1 }, { 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1 }, + { 0, 0, 1, 1, 1, 1, 0, 0, 1, 0, 0, 1 }, { 1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 1 }, + { 1, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 1 }, { 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 1 }, + { 0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 1, 0 }, { 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0 }, + { 1, 0, 1, 1, 0, 1, 1, 1, 0, 0, 1, 0 }, { 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0 }, + { 0, 1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0 }, { 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 0, 0 }, + { 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0 }, { 0, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0, 0 }, + { 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1 }, { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 }, + { 1, 1, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1 }, { 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1 }, + { 0, 0, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1 }, { 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1 }, + { 1, 0, 1, 1, 1, 0, 1, 1, 0, 0, 0, 1 }, { 0, 0, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1 }, + { 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1, 0 }, { 1, 0, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0 }, + { 1, 0, 0, 1, 0, 1, 0, 1, 1, 0, 1, 0 }, { 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0 }, + { 0, 1, 0, 1, 1, 0, 0, 1, 1, 1, 0, 0 }, { 1, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0 }, + { 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 0, 0 }, { 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0 }, + { 0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1 }, { 1, 1, 0, 1, 1, 1, 0, 1, 0, 1, 1, 1 }, + { 1, 1, 0, 1, 0, 0, 0, 1, 1, 0, 1, 1 }, { 0, 1, 0, 1, 1, 0, 0, 1, 0, 0, 1, 1 }, + { 0, 0, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1 }, { 1, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1 }, + { 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1 }, { 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1 }, + }; + + return edgecase < 128 ? edgeUses[edgecase] : edgeUses[127 - (edgecase - 128)]; +} + +VTKM_EXEC inline vtkm::UInt8 const* GetTriEdgeCases(vtkm::UInt8 edgecase) +{ + VTKM_STATIC_CONSTEXPR_ARRAY vtkm::UInt8 edgeCases[256][16] = { + // I expect we have some form on symmetry in this table + // that we can exploit to make it smaller + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 1, 0, 4, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 1, 0, 9, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 2, 5, 4, 8, 9, 5, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 1, 4, 1, 10, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 2, 0, 1, 10, 8, 0, 10, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 2, 5, 0, 9, 1, 10, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 3, 5, 1, 10, 5, 10, 9, 9, 10, 8, 0, 0, 0, 0, 0, 0 }, + { 1, 5, 11, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 2, 0, 4, 8, 5, 11, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 2, 9, 11, 1, 0, 9, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 3, 1, 4, 8, 1, 8, 11, 11, 8, 9, 0, 0, 0, 0, 0, 0 }, + { 2, 4, 5, 11, 10, 4, 11, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 3, 0, 5, 11, 0, 11, 8, 8, 11, 10, 0, 0, 0, 0, 0, 0 }, + { 3, 4, 0, 9, 4, 9, 10, 10, 9, 11, 0, 0, 0, 0, 0, 0 }, + { 2, 9, 11, 8, 11, 10, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 1, 2, 8, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 2, 2, 0, 4, 6, 2, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 2, 0, 9, 5, 8, 6, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 3, 2, 9, 5, 2, 5, 6, 6, 5, 4, 0, 0, 0, 0, 0, 0 }, + { 2, 8, 6, 2, 4, 1, 10, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 3, 10, 6, 2, 10, 2, 1, 1, 2, 0, 0, 0, 0, 0, 0, 0 }, + { 3, 9, 5, 0, 8, 6, 2, 1, 10, 4, 0, 0, 0, 0, 0, 0 }, + { 4, 2, 10, 6, 9, 10, 2, 9, 1, 10, 9, 5, 1, 0, 0, 0 }, + { 2, 5, 11, 1, 8, 6, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 3, 4, 6, 2, 4, 2, 0, 5, 11, 1, 0, 0, 0, 0, 0, 0 }, + { 3, 9, 11, 1, 9, 1, 0, 8, 6, 2, 0, 0, 0, 0, 0, 0 }, + { 4, 1, 9, 11, 1, 6, 9, 1, 4, 6, 6, 2, 9, 0, 0, 0 }, + { 3, 4, 5, 11, 4, 11, 10, 6, 2, 8, 0, 0, 0, 0, 0, 0 }, + { 4, 5, 11, 10, 5, 10, 2, 5, 2, 0, 6, 2, 10, 0, 0, 0 }, + { 4, 2, 8, 6, 9, 10, 0, 9, 11, 10, 10, 4, 0, 0, 0, 0 }, + { 3, 2, 10, 6, 2, 9, 10, 9, 11, 10, 0, 0, 0, 0, 0, 0 }, + { 1, 9, 2, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 2, 9, 2, 7, 0, 4, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 2, 0, 2, 7, 5, 0, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 3, 8, 2, 7, 8, 7, 4, 4, 7, 5, 0, 0, 0, 0, 0, 0 }, + { 2, 9, 2, 7, 1, 10, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 3, 0, 1, 10, 0, 10, 8, 2, 7, 9, 0, 0, 0, 0, 0, 0 }, + { 3, 0, 2, 7, 0, 7, 5, 1, 10, 4, 0, 0, 0, 0, 0, 0 }, + { 4, 1, 7, 5, 1, 8, 7, 1, 10, 8, 2, 7, 8, 0, 0, 0 }, + { 2, 5, 11, 1, 9, 2, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 3, 4, 8, 0, 5, 11, 1, 2, 7, 9, 0, 0, 0, 0, 0, 0 }, + { 3, 7, 11, 1, 7, 1, 2, 2, 1, 0, 0, 0, 0, 0, 0, 0 }, + { 4, 1, 7, 11, 4, 7, 1, 4, 2, 7, 4, 8, 2, 0, 0, 0 }, + { 3, 11, 10, 4, 11, 4, 5, 9, 2, 7, 0, 0, 0, 0, 0, 0 }, + { 4, 2, 7, 9, 0, 5, 8, 8, 5, 11, 8, 11, 10, 0, 0, 0 }, + { 4, 7, 0, 2, 7, 10, 0, 7, 11, 10, 10, 4, 0, 0, 0, 0 }, + { 3, 7, 8, 2, 7, 11, 8, 11, 10, 8, 0, 0, 0, 0, 0, 0 }, + { 2, 9, 8, 6, 7, 9, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 3, 9, 0, 4, 9, 4, 7, 7, 4, 6, 0, 0, 0, 0, 0, 0 }, + { 3, 0, 8, 6, 0, 6, 5, 5, 6, 7, 0, 0, 0, 0, 0, 0 }, + { 2, 5, 4, 7, 4, 6, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 3, 6, 7, 9, 6, 9, 8, 4, 1, 10, 0, 0, 0, 0, 0, 0 }, + { 4, 9, 6, 7, 9, 1, 6, 9, 0, 1, 1, 10, 6, 0, 0, 0 }, + { 4, 1, 10, 4, 0, 8, 5, 5, 8, 6, 5, 6, 7, 0, 0, 0 }, + { 3, 10, 5, 1, 10, 6, 5, 6, 7, 5, 0, 0, 0, 0, 0, 0 }, + { 3, 9, 8, 6, 9, 6, 7, 11, 1, 5, 0, 0, 0, 0, 0, 0 }, + { 4, 11, 1, 5, 9, 0, 7, 7, 0, 4, 7, 4, 6, 0, 0, 0 }, + { 4, 8, 1, 0, 8, 7, 1, 8, 6, 7, 11, 1, 7, 0, 0, 0 }, + { 3, 1, 7, 11, 1, 4, 7, 4, 6, 7, 0, 0, 0, 0, 0, 0 }, + { 4, 9, 8, 7, 8, 6, 7, 11, 4, 5, 11, 10, 4, 0, 0, 0 }, + { 5, 7, 0, 6, 7, 9, 0, 6, 0, 10, 5, 11, 0, 10, 0, 11 }, + { 5, 10, 0, 11, 10, 4, 0, 11, 0, 7, 8, 6, 0, 7, 0, 6 }, + { 2, 10, 7, 11, 6, 7, 10, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 1, 6, 10, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 2, 4, 8, 0, 10, 3, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 2, 0, 9, 5, 10, 3, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 3, 8, 9, 5, 8, 5, 4, 10, 3, 6, 0, 0, 0, 0, 0, 0 }, + { 2, 6, 4, 1, 3, 6, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 3, 6, 8, 0, 6, 0, 3, 3, 0, 1, 0, 0, 0, 0, 0, 0 }, + { 3, 1, 3, 6, 1, 6, 4, 0, 9, 5, 0, 0, 0, 0, 0, 0 }, + { 4, 5, 1, 3, 5, 3, 8, 5, 8, 9, 8, 3, 6, 0, 0, 0 }, + { 2, 11, 1, 5, 3, 6, 10, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 3, 5, 11, 1, 4, 8, 0, 3, 6, 10, 0, 0, 0, 0, 0, 0 }, + { 3, 1, 0, 9, 1, 9, 11, 3, 6, 10, 0, 0, 0, 0, 0, 0 }, + { 4, 3, 6, 10, 1, 4, 11, 11, 4, 8, 11, 8, 9, 0, 0, 0 }, + { 3, 11, 3, 6, 11, 6, 5, 5, 6, 4, 0, 0, 0, 0, 0, 0 }, + { 4, 11, 3, 6, 5, 11, 6, 5, 6, 8, 5, 8, 0, 0, 0, 0 }, + { 4, 0, 6, 4, 0, 11, 6, 0, 9, 11, 3, 6, 11, 0, 0, 0 }, + { 3, 6, 11, 3, 6, 8, 11, 8, 9, 11, 0, 0, 0, 0, 0, 0 }, + { 2, 3, 2, 8, 10, 3, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 3, 4, 10, 3, 4, 3, 0, 0, 3, 2, 0, 0, 0, 0, 0, 0 }, + { 3, 8, 10, 3, 8, 3, 2, 9, 5, 0, 0, 0, 0, 0, 0, 0 }, + { 4, 9, 3, 2, 9, 4, 3, 9, 5, 4, 10, 3, 4, 0, 0, 0 }, + { 3, 8, 4, 1, 8, 1, 2, 2, 1, 3, 0, 0, 0, 0, 0, 0 }, + { 2, 0, 1, 2, 2, 1, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 4, 5, 0, 9, 1, 2, 4, 1, 3, 2, 2, 8, 4, 0, 0, 0 }, + { 3, 5, 2, 9, 5, 1, 2, 1, 3, 2, 0, 0, 0, 0, 0, 0 }, + { 3, 3, 2, 8, 3, 8, 10, 1, 5, 11, 0, 0, 0, 0, 0, 0 }, + { 4, 5, 11, 1, 4, 10, 0, 0, 10, 3, 0, 3, 2, 0, 0, 0 }, + { 4, 2, 8, 10, 2, 10, 3, 0, 9, 1, 1, 9, 11, 0, 0, 0 }, + { 5, 11, 4, 9, 11, 1, 4, 9, 4, 2, 10, 3, 4, 2, 4, 3 }, + { 4, 8, 4, 5, 8, 5, 3, 8, 3, 2, 3, 5, 11, 0, 0, 0 }, + { 3, 11, 0, 5, 11, 3, 0, 3, 2, 0, 0, 0, 0, 0, 0, 0 }, + { 5, 2, 4, 3, 2, 8, 4, 3, 4, 11, 0, 9, 4, 11, 4, 9 }, + { 2, 11, 2, 9, 3, 2, 11, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 2, 2, 7, 9, 6, 10, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 3, 0, 4, 8, 2, 7, 9, 10, 3, 6, 0, 0, 0, 0, 0, 0 }, + { 3, 7, 5, 0, 7, 0, 2, 6, 10, 3, 0, 0, 0, 0, 0, 0 }, + { 4, 10, 3, 6, 8, 2, 4, 4, 2, 7, 4, 7, 5, 0, 0, 0 }, + { 3, 6, 4, 1, 6, 1, 3, 7, 9, 2, 0, 0, 0, 0, 0, 0 }, + { 4, 9, 2, 7, 0, 3, 8, 0, 1, 3, 3, 6, 8, 0, 0, 0 }, + { 4, 4, 1, 3, 4, 3, 6, 5, 0, 7, 7, 0, 2, 0, 0, 0 }, + { 5, 3, 8, 1, 3, 6, 8, 1, 8, 5, 2, 7, 8, 5, 8, 7 }, + { 3, 9, 2, 7, 11, 1, 5, 6, 10, 3, 0, 0, 0, 0, 0, 0 }, + { 4, 3, 6, 10, 5, 11, 1, 0, 4, 8, 2, 7, 9, 0, 0, 0 }, + { 4, 6, 10, 3, 7, 11, 2, 2, 11, 1, 2, 1, 0, 0, 0, 0 }, + { 5, 4, 8, 2, 4, 2, 7, 4, 7, 1, 11, 1, 7, 10, 3, 6 }, + { 4, 9, 2, 7, 11, 3, 5, 5, 3, 6, 5, 6, 4, 0, 0, 0 }, + { 5, 5, 11, 3, 5, 3, 6, 5, 6, 0, 8, 0, 6, 9, 2, 7 }, + { 5, 2, 11, 0, 2, 7, 11, 0, 11, 4, 3, 6, 11, 4, 11, 6 }, + { 4, 6, 11, 3, 6, 8, 11, 7, 11, 2, 2, 11, 8, 0, 0, 0 }, + { 3, 3, 7, 9, 3, 9, 10, 10, 9, 8, 0, 0, 0, 0, 0, 0 }, + { 4, 4, 10, 3, 0, 4, 3, 0, 3, 7, 0, 7, 9, 0, 0, 0 }, + { 4, 0, 8, 10, 0, 10, 7, 0, 7, 5, 7, 10, 3, 0, 0, 0 }, + { 3, 3, 4, 10, 3, 7, 4, 7, 5, 4, 0, 0, 0, 0, 0, 0 }, + { 4, 7, 9, 8, 7, 8, 1, 7, 1, 3, 4, 1, 8, 0, 0, 0 }, + { 3, 9, 3, 7, 9, 0, 3, 0, 1, 3, 0, 0, 0, 0, 0, 0 }, + { 5, 5, 8, 7, 5, 0, 8, 7, 8, 3, 4, 1, 8, 3, 8, 1 }, + { 2, 5, 3, 7, 1, 3, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 4, 5, 11, 1, 9, 10, 7, 9, 8, 10, 10, 3, 7, 0, 0, 0 }, + { 5, 0, 4, 10, 0, 10, 3, 0, 3, 9, 7, 9, 3, 5, 11, 1 }, + { 5, 10, 7, 8, 10, 3, 7, 8, 7, 0, 11, 1, 7, 0, 7, 1 }, + { 4, 3, 4, 10, 3, 7, 4, 1, 4, 11, 11, 4, 7, 0, 0, 0 }, + { 5, 5, 3, 4, 5, 11, 3, 4, 3, 8, 7, 9, 3, 8, 3, 9 }, + { 4, 11, 0, 5, 11, 3, 0, 9, 0, 7, 7, 0, 3, 0, 0, 0 }, + { 2, 0, 8, 4, 7, 11, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 1, 11, 3, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 1, 11, 7, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 2, 0, 4, 8, 7, 3, 11, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 2, 9, 5, 0, 7, 3, 11, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 3, 5, 4, 8, 5, 8, 9, 7, 3, 11, 0, 0, 0, 0, 0, 0 }, + { 2, 1, 10, 4, 11, 7, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 3, 10, 8, 0, 10, 0, 1, 11, 7, 3, 0, 0, 0, 0, 0, 0 }, + { 3, 0, 9, 5, 1, 10, 4, 7, 3, 11, 0, 0, 0, 0, 0, 0 }, + { 4, 7, 3, 11, 5, 1, 9, 9, 1, 10, 9, 10, 8, 0, 0, 0 }, + { 2, 5, 7, 3, 1, 5, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 3, 5, 7, 3, 5, 3, 1, 4, 8, 0, 0, 0, 0, 0, 0, 0 }, + { 3, 9, 7, 3, 9, 3, 0, 0, 3, 1, 0, 0, 0, 0, 0, 0 }, + { 4, 7, 8, 9, 7, 1, 8, 7, 3, 1, 4, 8, 1, 0, 0, 0 }, + { 3, 3, 10, 4, 3, 4, 7, 7, 4, 5, 0, 0, 0, 0, 0, 0 }, + { 4, 0, 10, 8, 0, 7, 10, 0, 5, 7, 7, 3, 10, 0, 0, 0 }, + { 4, 4, 3, 10, 0, 3, 4, 0, 7, 3, 0, 9, 7, 0, 0, 0 }, + { 3, 3, 9, 7, 3, 10, 9, 10, 8, 9, 0, 0, 0, 0, 0, 0 }, + { 2, 7, 3, 11, 2, 8, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 3, 2, 0, 4, 2, 4, 6, 3, 11, 7, 0, 0, 0, 0, 0, 0 }, + { 3, 5, 0, 9, 7, 3, 11, 8, 6, 2, 0, 0, 0, 0, 0, 0 }, + { 4, 11, 7, 3, 5, 6, 9, 5, 4, 6, 6, 2, 9, 0, 0, 0 }, + { 3, 4, 1, 10, 6, 2, 8, 11, 7, 3, 0, 0, 0, 0, 0, 0 }, + { 4, 7, 3, 11, 2, 1, 6, 2, 0, 1, 1, 10, 6, 0, 0, 0 }, + { 4, 0, 9, 5, 2, 8, 6, 1, 10, 4, 7, 3, 11, 0, 0, 0 }, + { 5, 9, 5, 1, 9, 1, 10, 9, 10, 2, 6, 2, 10, 7, 3, 11 }, + { 3, 3, 1, 5, 3, 5, 7, 2, 8, 6, 0, 0, 0, 0, 0, 0 }, + { 4, 5, 7, 1, 7, 3, 1, 4, 2, 0, 4, 6, 2, 0, 0, 0 }, + { 4, 8, 6, 2, 9, 7, 0, 0, 7, 3, 0, 3, 1, 0, 0, 0 }, + { 5, 6, 9, 4, 6, 2, 9, 4, 9, 1, 7, 3, 9, 1, 9, 3 }, + { 4, 8, 6, 2, 4, 7, 10, 4, 5, 7, 7, 3, 10, 0, 0, 0 }, + { 5, 7, 10, 5, 7, 3, 10, 5, 10, 0, 6, 2, 10, 0, 10, 2 }, + { 5, 0, 9, 7, 0, 7, 3, 0, 3, 4, 10, 4, 3, 8, 6, 2 }, + { 4, 3, 9, 7, 3, 10, 9, 2, 9, 6, 6, 9, 10, 0, 0, 0 }, + { 2, 11, 9, 2, 3, 11, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 3, 2, 3, 11, 2, 11, 9, 0, 4, 8, 0, 0, 0, 0, 0, 0 }, + { 3, 11, 5, 0, 11, 0, 3, 3, 0, 2, 0, 0, 0, 0, 0, 0 }, + { 4, 8, 5, 4, 8, 3, 5, 8, 2, 3, 3, 11, 5, 0, 0, 0 }, + { 3, 11, 9, 2, 11, 2, 3, 10, 4, 1, 0, 0, 0, 0, 0, 0 }, + { 4, 0, 1, 8, 1, 10, 8, 2, 11, 9, 2, 3, 11, 0, 0, 0 }, + { 4, 4, 1, 10, 0, 3, 5, 0, 2, 3, 3, 11, 5, 0, 0, 0 }, + { 5, 3, 5, 2, 3, 11, 5, 2, 5, 8, 1, 10, 5, 8, 5, 10 }, + { 3, 5, 9, 2, 5, 2, 1, 1, 2, 3, 0, 0, 0, 0, 0, 0 }, + { 4, 4, 8, 0, 5, 9, 1, 1, 9, 2, 1, 2, 3, 0, 0, 0 }, + { 2, 0, 2, 1, 2, 3, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 3, 8, 1, 4, 8, 2, 1, 2, 3, 1, 0, 0, 0, 0, 0, 0 }, + { 4, 9, 2, 3, 9, 3, 4, 9, 4, 5, 10, 4, 3, 0, 0, 0 }, + { 5, 8, 5, 10, 8, 0, 5, 10, 5, 3, 9, 2, 5, 3, 5, 2 }, + { 3, 4, 3, 10, 4, 0, 3, 0, 2, 3, 0, 0, 0, 0, 0, 0 }, + { 2, 3, 8, 2, 10, 8, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 3, 6, 3, 11, 6, 11, 8, 8, 11, 9, 0, 0, 0, 0, 0, 0 }, + { 4, 0, 4, 6, 0, 6, 11, 0, 11, 9, 3, 11, 6, 0, 0, 0 }, + { 4, 11, 6, 3, 5, 6, 11, 5, 8, 6, 5, 0, 8, 0, 0, 0 }, + { 3, 11, 6, 3, 11, 5, 6, 5, 4, 6, 0, 0, 0, 0, 0, 0 }, + { 4, 1, 10, 4, 11, 8, 3, 11, 9, 8, 8, 6, 3, 0, 0, 0 }, + { 5, 1, 6, 0, 1, 10, 6, 0, 6, 9, 3, 11, 6, 9, 6, 11 }, + { 5, 5, 0, 8, 5, 8, 6, 5, 6, 11, 3, 11, 6, 1, 10, 4 }, + { 4, 10, 5, 1, 10, 6, 5, 11, 5, 3, 3, 5, 6, 0, 0, 0 }, + { 4, 5, 3, 1, 5, 8, 3, 5, 9, 8, 8, 6, 3, 0, 0, 0 }, + { 5, 1, 9, 3, 1, 5, 9, 3, 9, 6, 0, 4, 9, 6, 9, 4 }, + { 3, 6, 0, 8, 6, 3, 0, 3, 1, 0, 0, 0, 0, 0, 0, 0 }, + { 2, 6, 1, 4, 3, 1, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 5, 8, 3, 9, 8, 6, 3, 9, 3, 5, 10, 4, 3, 5, 3, 4 }, + { 2, 0, 5, 9, 10, 6, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 4, 6, 0, 8, 6, 3, 0, 4, 0, 10, 10, 0, 3, 0, 0, 0 }, + { 1, 6, 3, 10, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 2, 10, 11, 7, 6, 10, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 3, 10, 11, 7, 10, 7, 6, 8, 0, 4, 0, 0, 0, 0, 0, 0 }, + { 3, 7, 6, 10, 7, 10, 11, 5, 0, 9, 0, 0, 0, 0, 0, 0 }, + { 4, 11, 7, 6, 11, 6, 10, 9, 5, 8, 8, 5, 4, 0, 0, 0 }, + { 3, 1, 11, 7, 1, 7, 4, 4, 7, 6, 0, 0, 0, 0, 0, 0 }, + { 4, 8, 0, 1, 8, 1, 7, 8, 7, 6, 11, 7, 1, 0, 0, 0 }, + { 4, 9, 5, 0, 7, 4, 11, 7, 6, 4, 4, 1, 11, 0, 0, 0 }, + { 5, 9, 1, 8, 9, 5, 1, 8, 1, 6, 11, 7, 1, 6, 1, 7 }, + { 3, 10, 1, 5, 10, 5, 6, 6, 5, 7, 0, 0, 0, 0, 0, 0 }, + { 4, 0, 4, 8, 5, 6, 1, 5, 7, 6, 6, 10, 1, 0, 0, 0 }, + { 4, 9, 7, 6, 9, 6, 1, 9, 1, 0, 1, 6, 10, 0, 0, 0 }, + { 5, 6, 1, 7, 6, 10, 1, 7, 1, 9, 4, 8, 1, 9, 1, 8 }, + { 2, 5, 7, 4, 4, 7, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 3, 0, 6, 8, 0, 5, 6, 5, 7, 6, 0, 0, 0, 0, 0, 0 }, + { 3, 9, 4, 0, 9, 7, 4, 7, 6, 4, 0, 0, 0, 0, 0, 0 }, + { 2, 9, 6, 8, 7, 6, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 3, 7, 2, 8, 7, 8, 11, 11, 8, 10, 0, 0, 0, 0, 0, 0 }, + { 4, 7, 2, 0, 7, 0, 10, 7, 10, 11, 10, 0, 4, 0, 0, 0 }, + { 4, 0, 9, 5, 8, 11, 2, 8, 10, 11, 11, 7, 2, 0, 0, 0 }, + { 5, 11, 2, 10, 11, 7, 2, 10, 2, 4, 9, 5, 2, 4, 2, 5 }, + { 4, 1, 11, 7, 4, 1, 7, 4, 7, 2, 4, 2, 8, 0, 0, 0 }, + { 3, 7, 1, 11, 7, 2, 1, 2, 0, 1, 0, 0, 0, 0, 0, 0 }, + { 5, 4, 1, 11, 4, 11, 7, 4, 7, 8, 2, 8, 7, 0, 9, 5 }, + { 4, 7, 1, 11, 7, 2, 1, 5, 1, 9, 9, 1, 2, 0, 0, 0 }, + { 4, 1, 5, 7, 1, 7, 8, 1, 8, 10, 2, 8, 7, 0, 0, 0 }, + { 5, 0, 10, 2, 0, 4, 10, 2, 10, 7, 1, 5, 10, 7, 10, 5 }, + { 5, 0, 7, 1, 0, 9, 7, 1, 7, 10, 2, 8, 7, 10, 7, 8 }, + { 2, 9, 7, 2, 1, 4, 10, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 3, 8, 7, 2, 8, 4, 7, 4, 5, 7, 0, 0, 0, 0, 0, 0 }, + { 2, 0, 7, 2, 5, 7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 4, 8, 7, 2, 8, 4, 7, 9, 7, 0, 0, 7, 4, 0, 0, 0 }, + { 1, 9, 7, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 3, 2, 6, 10, 2, 10, 9, 9, 10, 11, 0, 0, 0, 0, 0, 0 }, + { 4, 0, 4, 8, 2, 6, 9, 9, 6, 10, 9, 10, 11, 0, 0, 0 }, + { 4, 5, 10, 11, 5, 2, 10, 5, 0, 2, 6, 10, 2, 0, 0, 0 }, + { 5, 4, 2, 5, 4, 8, 2, 5, 2, 11, 6, 10, 2, 11, 2, 10 }, + { 4, 1, 11, 9, 1, 9, 6, 1, 6, 4, 6, 9, 2, 0, 0, 0 }, + { 5, 9, 6, 11, 9, 2, 6, 11, 6, 1, 8, 0, 6, 1, 6, 0 }, + { 5, 4, 11, 6, 4, 1, 11, 6, 11, 2, 5, 0, 11, 2, 11, 0 }, + { 2, 5, 1, 11, 8, 2, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 4, 2, 6, 10, 9, 2, 10, 9, 10, 1, 9, 1, 5, 0, 0, 0 }, + { 5, 9, 2, 6, 9, 6, 10, 9, 10, 5, 1, 5, 10, 0, 4, 8 }, + { 3, 10, 2, 6, 10, 1, 2, 1, 0, 2, 0, 0, 0, 0, 0, 0 }, + { 4, 10, 2, 6, 10, 1, 2, 8, 2, 4, 4, 2, 1, 0, 0, 0 }, + { 3, 2, 5, 9, 2, 6, 5, 6, 4, 5, 0, 0, 0, 0, 0, 0 }, + { 4, 2, 5, 9, 2, 6, 5, 0, 5, 8, 8, 5, 6, 0, 0, 0 }, + { 2, 2, 4, 0, 6, 4, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 1, 2, 6, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 2, 9, 8, 11, 11, 8, 10, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 3, 4, 9, 0, 4, 10, 9, 10, 11, 9, 0, 0, 0, 0, 0, 0 }, + { 3, 0, 11, 5, 0, 8, 11, 8, 10, 11, 0, 0, 0, 0, 0, 0 }, + { 2, 4, 11, 5, 10, 11, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 3, 1, 8, 4, 1, 11, 8, 11, 9, 8, 0, 0, 0, 0, 0, 0 }, + { 2, 9, 1, 11, 0, 1, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 4, 1, 8, 4, 1, 11, 8, 0, 8, 5, 5, 8, 11, 0, 0, 0 }, + { 1, 5, 1, 11, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 3, 5, 10, 1, 5, 9, 10, 9, 8, 10, 0, 0, 0, 0, 0, 0 }, + { 4, 4, 9, 0, 4, 10, 9, 5, 9, 1, 1, 9, 10, 0, 0, 0 }, + { 2, 0, 10, 1, 8, 10, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 1, 4, 10, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 2, 5, 8, 4, 9, 8, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 1, 0, 5, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 1, 0, 8, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, + }; + return edgeCases[edgecase]; +} + +VTKM_EXEC inline vtkm::Vec GetVertMap(vtkm::Id index) +{ + VTKM_STATIC_CONSTEXPR_ARRAY vtkm::Vec vertMap[12] = { + { 0, 1 }, { 2, 3 }, { 4, 5 }, { 6, 7 }, { 0, 2 }, { 1, 3 }, + { 4, 6 }, { 5, 7 }, { 0, 4 }, { 1, 5 }, { 2, 6 }, { 3, 7 }, + }; + return vertMap[index]; +} + +//x-axis +VTKM_EXEC inline vtkm::Id3 GetVertOffsets(SumXAxis, vtkm::UInt8 index) +{ + VTKM_STATIC_CONSTEXPR_ARRAY vtkm::Id3 offsetMap[8] = { + { 0, 0, 0 }, { 1, 0, 0 }, { 0, 1, 0 }, { 1, 1, 0 }, + { 0, 0, 1 }, { 1, 0, 1 }, { 0, 1, 1 }, { 1, 1, 1 }, + }; + return offsetMap[index]; +} +//y-axis +VTKM_EXEC inline vtkm::Id3 GetVertOffsets(SumYAxis, vtkm::UInt8 index) +{ + VTKM_STATIC_CONSTEXPR_ARRAY vtkm::Id3 offsetMap[8] = { + { 0, 0, 0 }, { 0, 1, 0 }, { 1, 0, 0 }, { 1, 1, 0 }, + { 0, 0, 1 }, { 0, 1, 1 }, { 1, 0, 1 }, { 1, 1, 1 }, + }; + return offsetMap[index]; +} + +} // namespace data +} +} +} +#endif diff --git a/vtkm/worklet/contour/MarchingCells.h b/vtkm/worklet/contour/MarchingCells.h index 7e47ea58f..0c6fbd060 100644 --- a/vtkm/worklet/contour/MarchingCells.h +++ b/vtkm/worklet/contour/MarchingCells.h @@ -597,19 +597,19 @@ struct GenerateNormals }; //---------------------------------------------------------------------------- -template vtkm::cont::CellSetSingleType<> execute( - const ValueType* isovalues, - const vtkm::Id numIsoValues, const CellSetType& cells, const CoordinateSystem& coordinateSystem, + const ValueType* isovalues, + const vtkm::Id numIsoValues, const vtkm::cont::ArrayHandle& inputField, vtkm::cont::ArrayHandle, StorageTagVertices>& vertices, vtkm::cont::ArrayHandle, StorageTagNormals>& normals, @@ -719,14 +719,14 @@ vtkm::cont::CellSetSingleType<> execute( //now that the vertices have been generated we can generate the normals if (sharedState.GenerateNormals) { - vtkm::cont::CastAndCall(coordinateSystem, - GenerateNormals{}, - invoker, - normals, - inputField, - cells, - sharedState.InterpolationEdgeIds, - sharedState.InterpolationWeights); + GenerateNormals genNorms; + genNorms(coordinateSystem, + invoker, + normals, + inputField, + cells, + sharedState.InterpolationEdgeIds, + sharedState.InterpolationWeights); } return outputCells; diff --git a/vtkm/worklet/gradient/StructuredPointGradient.h b/vtkm/worklet/gradient/StructuredPointGradient.h index 44a921af4..f0b2c9ea4 100644 --- a/vtkm/worklet/gradient/StructuredPointGradient.h +++ b/vtkm/worklet/gradient/StructuredPointGradient.h @@ -86,9 +86,16 @@ struct StructuredPointGradient : public vtkm::worklet::WorkletPointNeighborhood auto dy = inputField.Get(0, 1, 0) - inputField.Get(0, -1, 0); auto dz = inputField.Get(0, 0, 1) - inputField.Get(0, 0, -1); +#if (defined(VTKM_CUDA) && defined(VTKM_GCC)) +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wconversion" +#endif outputGradient[0] = static_cast(dx * r[0]); outputGradient[1] = static_cast(dy * r[1]); outputGradient[2] = static_cast(dz * r[2]); +#if (defined(VTKM_CUDA) && defined(VTKM_GCC)) +#pragma GCC diagnostic pop +#endif } //we need to pass the coordinates into this function, and instead From 1b256281f4c8310e33db27143743c373f8257a55 Mon Sep 17 00:00:00 2001 From: Robert Maynard Date: Mon, 16 Mar 2020 16:29:59 -0400 Subject: [PATCH 4/5] Convert contour worklet to expect isovalues to a std::vector This reduces an ugly component of the Contour API --- vtkm/filter/Contour.hxx | 11 +++-------- vtkm/worklet/Contour.h | 8 ++------ vtkm/worklet/contour/FlyingEdges.h | 5 ++--- vtkm/worklet/contour/MarchingCells.h | 12 +++++------- vtkm/worklet/testing/UnitTestContour.cxx | 14 ++++++-------- 5 files changed, 18 insertions(+), 32 deletions(-) diff --git a/vtkm/filter/Contour.hxx b/vtkm/filter/Contour.hxx index ae5495803..1a3532730 100644 --- a/vtkm/filter/Contour.hxx +++ b/vtkm/filter/Contour.hxx @@ -97,8 +97,7 @@ inline VTKM_CONT vtkm::cont::DataSet Contour::DoExecute( : !this->ComputeFastNormalsForUnstructured; if (this->GenerateNormals && generateHighQualityNormals) { - outputCells = this->Worklet.Run(&ivalues[0], - static_cast(ivalues.size()), + outputCells = this->Worklet.Run(ivalues, vtkm::filter::ApplyPolicyCellSet(cells, policy), coords.GetData(), field, @@ -107,12 +106,8 @@ inline VTKM_CONT vtkm::cont::DataSet Contour::DoExecute( } else { - outputCells = this->Worklet.Run(&ivalues[0], - static_cast(ivalues.size()), - vtkm::filter::ApplyPolicyCellSet(cells, policy), - coords.GetData(), - field, - vertices); + outputCells = this->Worklet.Run( + ivalues, vtkm::filter::ApplyPolicyCellSet(cells, policy), coords.GetData(), field, vertices); } if (this->GenerateNormals) diff --git a/vtkm/worklet/Contour.h b/vtkm/worklet/Contour.h index eb41754f9..0e13588f1 100644 --- a/vtkm/worklet/Contour.h +++ b/vtkm/worklet/Contour.h @@ -92,8 +92,7 @@ public: typename CoordinateType, typename StorageTagVertices> vtkm::cont::CellSetSingleType<> Run( - const ValueType* const isovalues, - const vtkm::Id numIsoValues, + const std::vector& isovalues, const CellSetType& cells, const CoordinateSystem& coordinateSystem, const vtkm::cont::ArrayHandle& input, @@ -108,7 +107,6 @@ public: coordinateSystem, outputCells, isovalues, - numIsoValues, input, vertices, normals, @@ -125,8 +123,7 @@ public: typename StorageTagVertices, typename StorageTagNormals> vtkm::cont::CellSetSingleType<> Run( - const ValueType* const isovalues, - const vtkm::Id numIsoValues, + const std::vector& isovalues, const CellSetType& cells, const CoordinateSystem& coordinateSystem, const vtkm::cont::ArrayHandle& input, @@ -141,7 +138,6 @@ public: coordinateSystem, outputCells, isovalues, - numIsoValues, input, vertices, normals, diff --git a/vtkm/worklet/contour/FlyingEdges.h b/vtkm/worklet/contour/FlyingEdges.h index 80d4502b9..835b2d137 100644 --- a/vtkm/worklet/contour/FlyingEdges.h +++ b/vtkm/worklet/contour/FlyingEdges.h @@ -86,8 +86,7 @@ template execute( const vtkm::cont::CellSetStructured<3>& cells, const vtkm::cont::ArrayHandleUniformPointCoordinates& coordinateSystem, - const ValueType* isovalues, - const vtkm::Id numIsoValues, + const std::vector& isovalues, const vtkm::cont::ArrayHandle& inputField, vtkm::cont::ArrayHandle, StorageTagVertices>& points, vtkm::cont::ArrayHandle, StorageTagNormals>& normals, @@ -127,7 +126,7 @@ vtkm::cont::CellSetSingleType<> execute( sharedState.CellIdMap.Shrink(0); vtkm::cont::ArrayHandle triangle_topology; - for (vtkm::Id i = 0; i < numIsoValues; ++i) + for (std::size_t i = 0; i < isovalues.size(); ++i) { ValueType isoval = isovalues[i]; diff --git a/vtkm/worklet/contour/MarchingCells.h b/vtkm/worklet/contour/MarchingCells.h index 0c6fbd060..c91164479 100644 --- a/vtkm/worklet/contour/MarchingCells.h +++ b/vtkm/worklet/contour/MarchingCells.h @@ -608,8 +608,7 @@ template execute( const CellSetType& cells, const CoordinateSystem& coordinateSystem, - const ValueType* isovalues, - const vtkm::Id numIsoValues, + const std::vector& isovalues, const vtkm::cont::ArrayHandle& inputField, vtkm::cont::ArrayHandle, StorageTagVertices>& vertices, vtkm::cont::ArrayHandle, StorageTagNormals>& normals, @@ -626,8 +625,7 @@ vtkm::cont::CellSetSingleType<> execute( // Setup the invoker vtkm::cont::Invoker invoker; - vtkm::cont::ArrayHandle isoValuesHandle = - vtkm::cont::make_ArrayHandle(isovalues, numIsoValues); + vtkm::cont::ArrayHandle isoValuesHandle = vtkm::cont::make_ArrayHandle(isovalues); // Call the ClassifyCell functor to compute the Marching Cubes case numbers // for each cell, and the number of vertices to be generated @@ -664,7 +662,7 @@ vtkm::cont::CellSetSingleType<> execute( triTable); } - if (numIsoValues <= 1 || !sharedState.MergeDuplicatePoints) + if (isovalues.size() <= 1 || !sharedState.MergeDuplicatePoints) { //release memory early that we are not going to need again contourIds.ReleaseResources(); } @@ -676,7 +674,7 @@ vtkm::cont::CellSetSingleType<> execute( // are updated. That is because MergeDuplicates will internally update // the InterpolationWeights and InterpolationOriginCellIds arrays to be the correct for the // output. But for InterpolationEdgeIds we need to do it manually once done - if (numIsoValues == 1) + if (isovalues.size() == 1) { marching_cells::MergeDuplicates(invoker, sharedState.InterpolationEdgeIds, //keys @@ -685,7 +683,7 @@ vtkm::cont::CellSetSingleType<> execute( originalCellIdsForPoints, //values connectivity); // computed using lower bounds } - else if (numIsoValues > 1) + else { marching_cells::MergeDuplicates( invoker, diff --git a/vtkm/worklet/testing/UnitTestContour.cxx b/vtkm/worklet/testing/UnitTestContour.cxx index 7dd49f6cc..d00c03e95 100644 --- a/vtkm/worklet/testing/UnitTestContour.cxx +++ b/vtkm/worklet/testing/UnitTestContour.cxx @@ -199,13 +199,12 @@ void TestContourUniformGrid() vtkm::worklet::Contour isosurfaceFilter; isosurfaceFilter.SetMergeDuplicatePoints(false); - vtkm::Float32 contourValue = 0.5f; + std::vector contourValue{ 0.5f }; vtkm::cont::ArrayHandle verticesArray; vtkm::cont::ArrayHandle normalsArray; vtkm::cont::ArrayHandle scalarsArray; - auto result = isosurfaceFilter.Run(&contourValue, - 1, + auto result = isosurfaceFilter.Run(contourValue, cellSet, dataSet.GetCoordinateSystem(), pointFieldArray, @@ -248,7 +247,7 @@ void TestContourExplicit() DataSetGenerator dataSetGenerator; vtkm::IdComponent Dimension = 10; - vtkm::Float32 contourValue = vtkm::Float32(.45); + std::vector contourValue{ 0.45f }; vtkm::cont::DataSet dataSet = dataSetGenerator.Make3DRadiantDataSet(Dimension); @@ -265,7 +264,7 @@ void TestContourExplicit() Contour.SetMergeDuplicatePoints(false); auto result = Contour.Run( - &contourValue, 1, cellSet, dataSet.GetCoordinateSystem(), contourArray, vertices, normals); + contourValue, cellSet, dataSet.GetCoordinateSystem(), contourArray, vertices, normals); DataHandle scalars; @@ -322,7 +321,7 @@ void TestContourClipped() vtkm::cont::ArrayHandle cellFieldArray; clipped.GetField("cellvar").GetData().CopyTo(cellFieldArray); - vtkm::Float32 contourValue = 0.5f; + std::vector contourValue{ 0.5f }; vtkm::cont::ArrayHandle verticesArray; vtkm::cont::ArrayHandle normalsArray; vtkm::cont::ArrayHandle scalarsArray; @@ -330,8 +329,7 @@ void TestContourClipped() vtkm::worklet::Contour isosurfaceFilter; isosurfaceFilter.SetMergeDuplicatePoints(false); - auto result = isosurfaceFilter.Run(&contourValue, - 1, + auto result = isosurfaceFilter.Run(contourValue, cellSet, clipped.GetCoordinateSystem(), pointFieldArray, From 91c56d680cc9b554de95ed351ebf41336ec286e5 Mon Sep 17 00:00:00 2001 From: Robert Maynard Date: Sun, 22 Mar 2020 13:25:26 -0400 Subject: [PATCH 5/5] Update Contour tests to be aware of flying edge output --- vtkm/filter/testing/UnitTestCleanGrid.cxx | 6 ++++++ vtkm/filter/testing/UnitTestContourFilter.cxx | 8 +++++--- vtkm/worklet/testing/UnitTestContour.cxx | 2 +- 3 files changed, 12 insertions(+), 4 deletions(-) diff --git a/vtkm/filter/testing/UnitTestCleanGrid.cxx b/vtkm/filter/testing/UnitTestCleanGrid.cxx index 17663bc0b..e52a9118d 100644 --- a/vtkm/filter/testing/UnitTestCleanGrid.cxx +++ b/vtkm/filter/testing/UnitTestCleanGrid.cxx @@ -72,6 +72,12 @@ void TestPointMerging() vtkm::cont::testing::MakeTestDataSet makeDataSet; vtkm::cont::DataSet baseData = makeDataSet.Make3DUniformDataSet3(vtkm::Id3(4, 4, 4)); + //Convert the baseData implicit points to explicit points, since the contour + //filter for uniform data always does point merging + vtkm::cont::ArrayHandle newcoords; + vtkm::cont::Algorithm::Copy(baseData.GetCoordinateSystem().GetData(), newcoords); + baseData.GetCoordinateSystem().SetData(newcoords); + vtkm::filter::Contour marchingCubes; marchingCubes.SetIsoValue(0.05); marchingCubes.SetMergeDuplicatePoints(false); diff --git a/vtkm/filter/testing/UnitTestContourFilter.cxx b/vtkm/filter/testing/UnitTestContourFilter.cxx index 1ee4f88a1..91daddf53 100644 --- a/vtkm/filter/testing/UnitTestContourFilter.cxx +++ b/vtkm/filter/testing/UnitTestContourFilter.cxx @@ -227,15 +227,17 @@ void TestContourUniformGrid() VTKM_TEST_ASSERT(cells.GetNumberOfCells() == 160, ""); } - //Now try with vertex merging disabled + //Now try with vertex merging disabled. Since this + //we use FlyingEdges we now which does point merging for free + //so we should see the number of points not change mc.SetMergeDuplicatePoints(false); mc.SetFieldsToPass(vtkm::filter::FieldSelection::MODE_ALL); result = mc.Execute(dataSet); { vtkm::cont::CoordinateSystem coords = result.GetCoordinateSystem(); - VTKM_TEST_ASSERT(coords.GetNumberOfPoints() == 480, - "Should have less coordinates than the unmerged version"); + VTKM_TEST_ASSERT(coords.GetNumberOfPoints() == 72, + "Shouldn't have less coordinates than the unmerged version"); //verify that the number of cells is correct (160) vtkm::cont::DynamicCellSet dcells = result.GetCellSet(); diff --git a/vtkm/worklet/testing/UnitTestContour.cxx b/vtkm/worklet/testing/UnitTestContour.cxx index d00c03e95..2855d5430 100644 --- a/vtkm/worklet/testing/UnitTestContour.cxx +++ b/vtkm/worklet/testing/UnitTestContour.cxx @@ -233,7 +233,7 @@ void TestContourUniformGrid() VTKM_TEST_ASSERT(result.GetNumberOfCells() == 160); - VTKM_TEST_ASSERT(verticesArray.GetNumberOfValues() == 480); + VTKM_TEST_ASSERT(verticesArray.GetNumberOfValues() == 72); } void TestContourExplicit()