Shuenhoy/vtk-m
1
0
Fork 0
mirror of https://gitlab.kitware.com/vtk/vtk-m synced 2024-09-08 13:23:51 +00:00
Code Issues 2 Actions Packages Projects Releases Wiki Activity

Add a point gradient worklet.

Browse Source
This commit is contained in:
Robert Maynard 2016-12-01 14:20:47 -05:00
parent 06af2fbcbf
commit 134f24496b
5 changed files with 373 additions and 22 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
vtkm/exec/arg/ThreadIndicesTopologyMap.h
View File

@ -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<CellShapeTag>::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)

195
vtkm/worklet/Gradient.h
View File

@ -21,14 +21,18 @@
#ifndef vtk_m_worklet_Gradient_h
#define vtk_m_worklet_Gradient_h
#include <vtkm/exec/ExecutionWholeArray.h>
#include <vtkm/exec/CellDerivative.h>
#include <vtkm/exec/ExecutionWholeArray.h>
#include <vtkm/exec/ParametricCoordinates.h>
#include <vtkm/CellTraits.h>
#include <vtkm/VecFromPortal.h>
#include <vtkm/VecFromPortalPermute.h>
#include <vtkm/VecTraits.h>
#include <vtkm/worklet/WorkletMapTopology.h>
#include <vtkm/exec/CellDerivative.h>
namespace vtkm {
namespace worklet {
@ -90,6 +94,7 @@ struct CellGradient : vtkm::worklet::WorkletMapPointToCell
{
vtkm::Vec<vtkm::FloatDefault, 3> 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<InValueType>::DimensionalityTag;
// using OutValueType = typename VecTraits<typename FieldOutType::ComponentType>::ComponentType;
// using OutDimensionTag = typename TypeTraits<OutValueType>::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<Point,Cell>,
WholeArrayIn<Vec3> pointCoordinates,
WholeArrayIn<FieldCommon> inputField,
FieldOutPoint<GradientOutTypes> outputField);
typedef void ExecutionSignature(CellCount, CellIndices, WorkIndex, _2, _3, _4, _5);
typedef _1 InputDomain;
template <typename FromIndexType,
typename CellSetInType,
typename WholeCoordinatesIn,
typename WholeFieldIn,
typename FieldOutType>
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> 'T' type.
//For example:
// input is float => output is vtkm::Vec<float>
// input is vtkm::Vec<float> => output is vtkm::Vec< vtkm::Vec< float > >
//Grab the dimension tag for the input
using InValueType = typename WholeFieldIn::ValueType;
using InDimensionTag = typename TypeTraits<InValueType>::DimensionalityTag;
//grab the dimension tag for the output component type
using OutValueType = typename FieldOutType::ComponentType;
using OutDimensionTag = typename TypeTraits<OutValueType>::DimensionalityTag;
//Verify that input and output dimension tags match
using Matches = typename std::is_same<InDimensionTag, OutDimensionTag>::type;
this->Compute(numCells, cellIds, pointId, geometry, pointCoordinates,
inputField, outputField, Matches());
}
template <typename FromIndexType,
typename CellSetInType,
typename WholeCoordinatesIn,
typename WholeFieldIn,
typename FieldOutType>
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<CellSetInType>;
using ValueType = typename WholeFieldIn::ValueType;
using CellShapeTag = typename CellSetInType::CellShapeTag;
vtkm::Vec<ValueType, 3> 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<OutValueType>(gradient[0] / numCells);
outputField[1] = static_cast<OutValueType>(gradient[1] / numCells);
outputField[2] = static_cast<OutValueType>(gradient[2] / numCells);
}
template <typename FromIndexType,
typename CellSetInType,
typename WholeCoordinatesIn,
typename WholeFieldIn,
typename FieldOutType>
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 <typename CellShapeTag, typename PointCoordVecType,
typename FieldInVecType, typename OutValueType>
inline VTKM_EXEC
void ComputeGradient(CellShapeTag cellShape,
const vtkm::IdComponent& pointIndexForCell,
const PointCoordVecType& wCoords,
const FieldInVecType& field,
vtkm::Vec<OutValueType, 3>& gradient) const
{
vtkm::Vec<vtkm::FloatDefault, 3> 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<typename CellSetInType>
VTKM_EXEC
vtkm::IdComponent GetPointIndexForCell(
const vtkm::exec::arg::ThreadIndicesTopologyMap<CellSetInType>& 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<typename CellSetInType, typename WholeFieldIn>
VTKM_EXEC
typename vtkm::exec::arg::Fetch<
vtkm::exec::arg::FetchTagArrayTopologyMapIn,
vtkm::exec::arg::AspectTagDefault,
vtkm::exec::arg::ThreadIndicesTopologyMap<CellSetInType>,
typename WholeFieldIn::PortalType>::ValueType
GetValues(const vtkm::exec::arg::ThreadIndicesTopologyMap<CellSetInType>& 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<CellSetInType>,
ExecObjectType>;
Fetch fetch;
return fetch.Load(indices,in.GetPortal());
}
};
}
} // namespace vtkm::worklet

1
vtkm/worklet/testing/CMakeLists.txt
View File

@ -28,6 +28,7 @@ set(unit_tests
UnitTestFieldStatistics.cxx
UnitTestMarchingCubes.cxx
UnitTestPointElevation.cxx
UnitTestPointGradient.cxx
UnitTestRemoveUnusedPoints.cxx
UnitTestScatterCounting.cxx
UnitTestSplatKernels.cxx

3
vtkm/worklet/testing/UnitTestCellGradient.cxx
View File

@ -120,9 +120,6 @@ void TestCellGradientUniform3DWithVectorField()
vtkm::Vec< vtkm::Vec<vtkm::Float64,3>, 3> e = expected[i];
vtkm::Vec< vtkm::Vec<vtkm::Float64,3>, 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");

190
vtkm/worklet/testing/UnitTestPointGradient.cxx Normal file
View File

@ -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 <vtkm/worklet/Gradient.h>
#include <vtkm/worklet/DispatcherMapTopology.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/cont/testing/MakeTestDataSet.h>
namespace {
template<typename T>
struct PointGrad
{
PointGrad(const vtkm::cont::DataSet& data,
const std::string& fieldName,
vtkm::cont::ArrayHandle< vtkm::Vec<T,3> >& result):
Data(data),
FieldName(fieldName),
Result(result)
{
}
template<typename CellSetType>
void operator()(const CellSetType& cellset ) const
{
vtkm::worklet::DispatcherMapTopology<vtkm::worklet::PointGradient> 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<T,3> > 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<vtkm::Float32,3> > result;
PointGrad<vtkm::Float32> func(dataSet, "pointvar", result);
vtkm::cont::CastAndCall(dataSet.GetCellSet(), func);
vtkm::Vec<vtkm::Float32,3> 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<vtkm::Float64,3> > result;
PointGrad<vtkm::Float64> func(dataSet, "pointvar", result);
vtkm::cont::CastAndCall(dataSet.GetCellSet(), func);
vtkm::Vec<vtkm::Float64,3> 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<vtkm::Float64,3> > 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<vtkm::Float64,3> > 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<vtkm::Float64,3>, 3> > result;
PointGrad< vtkm::Vec<vtkm::Float64,3> > func(dataSet, "vec_pointvar", result);
vtkm::cont::CastAndCall(dataSet.GetCellSet(), func);
vtkm::Vec< vtkm::Vec<vtkm::Float64,3>, 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<vtkm::Float64,3>, 3> e = expected[i];
vtkm::Vec< vtkm::Vec<vtkm::Float64,3>, 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<vtkm::Float32,3> > result;
PointGrad<vtkm::Float32> func(dataSet, "pointvar", result);
vtkm::cont::CastAndCall(dataSet.GetCellSet(), func);
vtkm::Vec<vtkm::Float32,3> 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);
}