From 134f24496b753c5fd2f1b5833029571db6cbc0fb Mon Sep 17 00:00:00 2001 From: Robert Maynard Date: Thu, 1 Dec 2016 14:20:47 -0500 Subject: [PATCH] Add a point gradient worklet. --- vtkm/exec/arg/ThreadIndicesTopologyMap.h | 6 +- vtkm/worklet/Gradient.h | 195 ++++++++++++++++-- vtkm/worklet/testing/CMakeLists.txt | 1 + vtkm/worklet/testing/UnitTestCellGradient.cxx | 3 - .../worklet/testing/UnitTestPointGradient.cxx | 190 +++++++++++++++++ 5 files changed, 373 insertions(+), 22 deletions(-) create mode 100644 vtkm/worklet/testing/UnitTestPointGradient.cxx diff --git a/vtkm/exec/arg/ThreadIndicesTopologyMap.h b/vtkm/exec/arg/ThreadIndicesTopologyMap.h index 0a5aaf681..43232e0c4 100644 --- a/vtkm/exec/arg/ThreadIndicesTopologyMap.h +++ b/vtkm/exec/arg/ThreadIndicesTopologyMap.h @@ -79,7 +79,7 @@ public: const OutToInArrayType& inToOut, const VisitArrayType& visit, const ConnectivityType& connectivity, - vtkm::Id globalThreadIndexOffset=0) + vtkm::Id globalThreadIndexOffset=0) : Superclass(threadIndex, inToOut.Get(threadIndex), visit.Get(threadIndex), globalThreadIndexOffset), CellShape(detail::CellShapeInitializer::GetDefault()) @@ -212,7 +212,7 @@ public: const OutToInArrayType& inToOut, const VisitArrayType& visit, const ConnectivityType& connectivity, - vtkm::Id globalThreadIndexOffset=0) + vtkm::Id globalThreadIndexOffset=0) { this->InputIndex = inToOut.Get(threadIndex); @@ -251,7 +251,7 @@ public: VTKM_SUPPRESS_EXEC_WARNINGS VTKM_EXEC ThreadIndicesTopologyMap(vtkm::Id threadIndex, - vtkm::Id inIndex, + vtkm::Id vtkmNotUsed(inIndex), vtkm::IdComponent visitIndex, const ConnectivityType& connectivity, vtkm::Id globalThreadIndexOffset=0) diff --git a/vtkm/worklet/Gradient.h b/vtkm/worklet/Gradient.h index 1eb02f2b5..c4895c7cc 100644 --- a/vtkm/worklet/Gradient.h +++ b/vtkm/worklet/Gradient.h @@ -21,14 +21,18 @@ #ifndef vtk_m_worklet_Gradient_h #define vtk_m_worklet_Gradient_h -#include #include +#include #include #include +#include #include -#include #include + + +#include + namespace vtkm { namespace worklet { @@ -90,6 +94,7 @@ struct CellGradient : vtkm::worklet::WorkletMapPointToCell { vtkm::Vec center = vtkm::exec::ParametricCoordinatesCenter(pointCount, shape, *this); + outputField = vtkm::exec::CellDerivative(field, wCoords, center, shape, *this); } @@ -105,23 +110,181 @@ struct CellGradient : vtkm::worklet::WorkletMapPointToCell std::false_type) const { //this is invalid - // std::cout << "calling invalid compute" << std::endl; - // using InValueType = typename FieldInVecType::ComponentType; - // using InDimensionTag = typename TypeTraits::DimensionalityTag; - - // using OutValueType = typename VecTraits::ComponentType; - // using OutDimensionTag = typename TypeTraits::DimensionalityTag; - - // std::cout << typeid(InDimensionTag).name() << '\n'; - // std::cout << typeid(FieldInVecType).name() << '\n'; - - // std::cout << typeid(OutDimensionTag).name() << '\n'; - // std::cout << typeid(FieldOutType).name() << '\n'; - - // std::cout << std::endl; } }; +struct PointGradient : public vtkm::worklet::WorkletMapCellToPoint +{ + typedef void ControlSignature(CellSetIn, + WholeCellSetIn, + WholeArrayIn pointCoordinates, + WholeArrayIn inputField, + FieldOutPoint outputField); + + typedef void ExecutionSignature(CellCount, CellIndices, WorkIndex, _2, _3, _4, _5); + typedef _1 InputDomain; + + template + VTKM_EXEC void operator()(const vtkm::IdComponent& numCells, + const FromIndexType& cellIds, + const vtkm::Id& pointId, + const CellSetInType& geometry, + const WholeCoordinatesIn& pointCoordinates, + const WholeFieldIn& inputField, + FieldOutType& outputField) const + { + //To confirm that we have the proper input and output types we need + //to verify that input type matches the output vtkm::Vec 'T' type. + //For example: + // input is float => output is vtkm::Vec + // input is vtkm::Vec => output is vtkm::Vec< vtkm::Vec< float > > + + //Grab the dimension tag for the input + using InValueType = typename WholeFieldIn::ValueType; + using InDimensionTag = typename TypeTraits::DimensionalityTag; + + //grab the dimension tag for the output component type + using OutValueType = typename FieldOutType::ComponentType; + using OutDimensionTag = typename TypeTraits::DimensionalityTag; + + //Verify that input and output dimension tags match + using Matches = typename std::is_same::type; + this->Compute(numCells, cellIds, pointId, geometry, pointCoordinates, + inputField, outputField, Matches()); + } + + template + VTKM_EXEC void Compute(const vtkm::IdComponent& numCells, + const FromIndexType& cellIds, + const vtkm::Id& pointId, + const CellSetInType& geometry, + const WholeCoordinatesIn& pointCoordinates, + const WholeFieldIn& inputField, + FieldOutType& outputField, + std::true_type) const + { + using ThreadIndices = vtkm::exec::arg::ThreadIndicesTopologyMap; + using ValueType = typename WholeFieldIn::ValueType; + using CellShapeTag = typename CellSetInType::CellShapeTag; + + vtkm::Vec gradient( ValueType(0.0) ); + for (vtkm::IdComponent i = 0; i < numCells; ++i) + { + const vtkm::Id cellId = cellIds[i]; + ThreadIndices cellIndices(cellId, cellId, 0, geometry); + + const CellShapeTag cellShape = cellIndices.GetCellShape(); + + // compute the parametric coordinates for the current point + const auto wCoords = this->GetValues(cellIndices, pointCoordinates); + const auto field = this->GetValues(cellIndices, inputField); + + const vtkm::IdComponent pointIndexForCell = + this->GetPointIndexForCell(cellIndices, pointId); + + this->ComputeGradient(cellShape, pointIndexForCell, wCoords, field, gradient); + } + + using OutValueType = typename FieldOutType::ComponentType; + outputField[0] = static_cast(gradient[0] / numCells); + outputField[1] = static_cast(gradient[1] / numCells); + outputField[2] = static_cast(gradient[2] / numCells); + } + + template + VTKM_EXEC void Compute(const vtkm::IdComponent&, + const FromIndexType&, + const vtkm::Id&, + const CellSetInType&, + const WholeCoordinatesIn&, + const WholeFieldIn&, + FieldOutType&, + std::false_type) const + { + //this is invalid, as the input and output types don't match. + //e.g input is float => output is vtkm::Vec< vtkm::Vec< float > > + } + +private: + template + inline VTKM_EXEC + void ComputeGradient(CellShapeTag cellShape, + const vtkm::IdComponent& pointIndexForCell, + const PointCoordVecType& wCoords, + const FieldInVecType& field, + vtkm::Vec& gradient) const + { + vtkm::Vec pCoords; + vtkm::exec::ParametricCoordinatesPoint( + wCoords.GetNumberOfComponents(), pointIndexForCell, pCoords, cellShape, *this); + + //we need to add this to a return value + gradient += vtkm::exec::CellDerivative(field, + wCoords, pCoords, + cellShape, *this); + } + + template + VTKM_EXEC + vtkm::IdComponent GetPointIndexForCell( + const vtkm::exec::arg::ThreadIndicesTopologyMap& indices, + vtkm::Id pointId) const + { + vtkm::IdComponent result = 0; + const auto& topo = indices.GetIndicesFrom(); + for (vtkm::IdComponent i = 0; i < topo.GetNumberOfComponents(); ++i) + { + if (topo[i] == pointId) + { + result = i; + } + } + return result; + } + + //This is fairly complex so that we can trigger code to extract + //VecRectilinearPointCoordinates when using structured connectivity, and + //uniform point coordinates. + //c++14 would make the return type simply auto + template + VTKM_EXEC + typename vtkm::exec::arg::Fetch< + vtkm::exec::arg::FetchTagArrayTopologyMapIn, + vtkm::exec::arg::AspectTagDefault, + vtkm::exec::arg::ThreadIndicesTopologyMap, + typename WholeFieldIn::PortalType>::ValueType + GetValues(const vtkm::exec::arg::ThreadIndicesTopologyMap& indices, + const WholeFieldIn& in) const + { + //the current problem is that when the topology is structured + //we are passing in an vtkm::Id when it wants a Id2 or an Id3 that + //represents the flat index of the topology + using ExecObjectType = typename WholeFieldIn::PortalType; + using Fetch = vtkm::exec::arg::Fetch< + vtkm::exec::arg::FetchTagArrayTopologyMapIn, + vtkm::exec::arg::AspectTagDefault, + vtkm::exec::arg::ThreadIndicesTopologyMap, + ExecObjectType>; + Fetch fetch; + return fetch.Load(indices,in.GetPortal()); + } + +}; + + } } // namespace vtkm::worklet diff --git a/vtkm/worklet/testing/CMakeLists.txt b/vtkm/worklet/testing/CMakeLists.txt index a19e316d6..776839558 100644 --- a/vtkm/worklet/testing/CMakeLists.txt +++ b/vtkm/worklet/testing/CMakeLists.txt @@ -28,6 +28,7 @@ set(unit_tests UnitTestFieldStatistics.cxx UnitTestMarchingCubes.cxx UnitTestPointElevation.cxx + UnitTestPointGradient.cxx UnitTestRemoveUnusedPoints.cxx UnitTestScatterCounting.cxx UnitTestSplatKernels.cxx diff --git a/vtkm/worklet/testing/UnitTestCellGradient.cxx b/vtkm/worklet/testing/UnitTestCellGradient.cxx index 01de8de0e..7403ba2da 100644 --- a/vtkm/worklet/testing/UnitTestCellGradient.cxx +++ b/vtkm/worklet/testing/UnitTestCellGradient.cxx @@ -120,9 +120,6 @@ void TestCellGradientUniform3DWithVectorField() vtkm::Vec< vtkm::Vec, 3> e = expected[i]; vtkm::Vec< vtkm::Vec, 3> r = result.GetPortalConstControl().Get(i); - std::cout << "0: " << e[0] << " " << r[0] << std::endl; - std::cout << "1: " << e[1] << " " << r[1] << std::endl; - std::cout << "2: " << e[2] << " " << r[2] << std::endl; VTKM_TEST_ASSERT( test_equal(e[0],r[0]), "Wrong result for vec field CellGradient worklet on 3D uniform data"); diff --git a/vtkm/worklet/testing/UnitTestPointGradient.cxx b/vtkm/worklet/testing/UnitTestPointGradient.cxx new file mode 100644 index 000000000..31b719bf5 --- /dev/null +++ b/vtkm/worklet/testing/UnitTestPointGradient.cxx @@ -0,0 +1,190 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.txt for details. +// This software is distributed WITHOUT ANY WARRANTY; without even +// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +// PURPOSE. See the above copyright notice for more information. +// +// Copyright 2014 Sandia Corporation. +// Copyright 2014 UT-Battelle, LLC. +// Copyright 2014 Los Alamos National Security. +// +// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, +// the U.S. Government retains certain rights in this software. +// +// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National +// Laboratory (LANL), the U.S. Government retains certain rights in +// this software. +//============================================================================ + +#include +#include + +#include +#include + +namespace { + +template +struct PointGrad +{ + PointGrad(const vtkm::cont::DataSet& data, + const std::string& fieldName, + vtkm::cont::ArrayHandle< vtkm::Vec >& result): + Data(data), + FieldName(fieldName), + Result(result) + { + } + + template + void operator()(const CellSetType& cellset ) const + { + vtkm::worklet::DispatcherMapTopology dispatcher; + dispatcher.Invoke(cellset, //topology to iterate on a per point basis + cellset, //whole cellset in + this->Data.GetCoordinateSystem(), + this->Data.GetField(this->FieldName), + this->Result); + } + + vtkm::cont::DataSet Data; + std::string FieldName; + vtkm::cont::ArrayHandle< vtkm::Vec > Result; +}; + +void TestPointGradientUniform2D() +{ + std::cout << "Testing PointGradient Worklet on 2D structured data" << std::endl; + + vtkm::cont::testing::MakeTestDataSet testDataSet; + vtkm::cont::DataSet dataSet = testDataSet.Make2DUniformDataSet0(); + + vtkm::cont::ArrayHandle< vtkm::Vec > result; + + PointGrad func(dataSet, "pointvar", result); + vtkm::cont::CastAndCall(dataSet.GetCellSet(), func); + + vtkm::Vec expected[2] = { {10,30,0}, {10,30,0} }; + for (int i = 0; i < 2; ++i) + { + VTKM_TEST_ASSERT( + test_equal(result.GetPortalConstControl().Get(i), expected[i]), + "Wrong result for PointGradient worklet on 2D uniform data"); + } +} + + +void TestPointGradientUniform3D() +{ + std::cout << "Testing PointGradient Worklet on 3D strucutred data" << std::endl; + + vtkm::cont::testing::MakeTestDataSet testDataSet; + vtkm::cont::DataSet dataSet = testDataSet.Make3DUniformDataSet0(); + + vtkm::cont::ArrayHandle< vtkm::Vec > result; + + PointGrad func(dataSet, "pointvar", result); + vtkm::cont::CastAndCall(dataSet.GetCellSet(), func); + + vtkm::Vec expected[4] = { {10.0,30,60.1}, + {10.0,30.1,60.1}, + {10.0,30.1,60.2}, + {10.1,30,60.2}, + }; + for (int i = 0; i < 4; ++i) + { + VTKM_TEST_ASSERT( + test_equal(result.GetPortalConstControl().Get(i), expected[i]), + "Wrong result for PointGradient worklet on 3D uniform data"); + } + +} + +void TestPointGradientUniform3DWithVectorField() +{ + std::cout << "Testing PointGradient Worklet with a vector field on 3D strucutred data" << std::endl; + vtkm::cont::testing::MakeTestDataSet testDataSet; + vtkm::cont::DataSet dataSet = testDataSet.Make3DUniformDataSet0(); + + //Verify that we can compute the gradient of a 3 component vector + const int nVerts = 18; + vtkm::Float64 vars[nVerts] = {10.1, 20.1, 30.1, 40.1, 50.2, + 60.2, 70.2, 80.2, 90.3, 100.3, + 110.3, 120.3, 130.4, 140.4, + 150.4, 160.4, 170.5, 180.5}; + std::vector< vtkm::Vec > vec(18); + for(std::size_t i=0; i < vec.size(); ++i) + { + vec[i] = vtkm::make_Vec(vars[i],vars[i],vars[i]); + } + vtkm::cont::ArrayHandle< vtkm::Vec > input = + vtkm::cont::make_ArrayHandle(vec); + //we need to add Vec3 array to the dataset + vtkm::cont::DataSetFieldAdd::AddPointField(dataSet, "vec_pointvar", input); + + vtkm::cont::ArrayHandle< vtkm::Vec< vtkm::Vec, 3> > result; + PointGrad< vtkm::Vec > func(dataSet, "vec_pointvar", result); + vtkm::cont::CastAndCall(dataSet.GetCellSet(), func); + + vtkm::Vec< vtkm::Vec, 3> expected[4] = { + { {10.0,10.0,10.0}, {30.0,30.0,30.0}, {60.1,60.1,60.1} }, + { {10.0,10.0,10.0}, {30.1,30.1,30.1}, {60.1,60.1,60.1} }, + { {10.0,10.0,10.0}, {30.1,30.1,30.1}, {60.2,60.2,60.2} }, + { {10.1,10.1,10.1}, {30.0,30.0,30.0}, {60.2,60.2,60.2} } + }; + for (int i = 0; i < 4; ++i) + { + vtkm::Vec< vtkm::Vec, 3> e = expected[i]; + vtkm::Vec< vtkm::Vec, 3> r = result.GetPortalConstControl().Get(i); + + VTKM_TEST_ASSERT( + test_equal(e[0],r[0]), + "Wrong result for vec field PointGradient worklet on 3D uniform data"); + VTKM_TEST_ASSERT( + test_equal(e[1],r[1]), + "Wrong result for vec field PointGradient worklet on 3D uniform data"); + VTKM_TEST_ASSERT( + test_equal(e[2],r[2]), + "Wrong result for vec field PointGradient worklet on 3D uniform data"); + } + +} + +void TestPointGradientExplicit() +{ + std::cout << "Testing PointGradient Worklet on Explicit data" << std::endl; + + vtkm::cont::testing::MakeTestDataSet testDataSet; + vtkm::cont::DataSet dataSet = testDataSet.Make3DExplicitDataSet0(); + + vtkm::cont::ArrayHandle< vtkm::Vec > result; + + PointGrad func(dataSet, "pointvar", result); + vtkm::cont::CastAndCall(dataSet.GetCellSet(), func); + + vtkm::Vec expected[2] = { {10.f,10.1f,0.0f}, {10.f,10.1f,0.0f} }; + for (int i = 0; i < 2; ++i) + { + VTKM_TEST_ASSERT( + test_equal(result.GetPortalConstControl().Get(i), expected[i]), + "Wrong result for PointGradient worklet on 3D explicit data"); + } +} + + +void TestPointGradient() +{ + TestPointGradientUniform2D(); + TestPointGradientUniform3D(); + TestPointGradientUniform3DWithVectorField(); + TestPointGradientExplicit(); +} + +} + +int UnitTestPointGradient(int, char *[]) +{ + return vtkm::cont::testing::Testing::Run(TestPointGradient); +}