Commit 556b9227 authored by Thomas Otahal's avatar Thomas Otahal

Created a filter for the vector magnitude worklet.

Modified the vector magnitude worklet to accept VecAll instead
of Vec3 and return Scalar. Modified the Magnitude() and Sqrt() functions
to return FloatDefault for all inputs except for Float64.

Perhaps we should modify other functions in Math.h and VectorAnalysis.h to
return float types for intergral arguments instead of integral types?
parent 589285eb
......@@ -740,6 +740,18 @@ vtkm::Float64 Pow(vtkm::Float64 x, vtkm::Float64 y) {
/// Compute the square root of \p x.
///
template<typename T>
static inline VTKM_EXEC_CONT
vtkm::FloatDefault Sqrt(T x) {
#ifdef VTKM_CUDA
if(std::is_same< vtkm::FloatDefault, vtkm::Float64 >::value)
return VTKM_CUDA_MATH_FUNCTION_64(sqrt)(static_cast<vtkm::FloatDefault>(x));
else
return VTKM_CUDA_MATH_FUNCTION_32(sqrt)(static_cast<vtkm::FloatDefault>(x));
#else
return std::sqrt(static_cast<vtkm::FloatDefault>(x));
#endif
}
static inline VTKM_EXEC_CONT
vtkm::Float32 Sqrt(vtkm::Float32 x) {
#ifdef VTKM_CUDA
......@@ -758,8 +770,8 @@ vtkm::Float64 Sqrt(vtkm::Float64 x) {
}
template<typename T, vtkm::IdComponent N>
static inline VTKM_EXEC_CONT
vtkm::Vec<T,N> Sqrt(const vtkm::Vec<T,N> &x) {
vtkm::Vec<T,N> result;
vtkm::Vec<vtkm::FloatDefault,N> Sqrt(const vtkm::Vec<T,N> &x) {
vtkm::Vec<vtkm::FloatDefault,N> result;
for (vtkm::IdComponent index = 0; index < N; index++)
{
result[index] = vtkm::Sqrt(x[index]);
......@@ -768,24 +780,52 @@ vtkm::Vec<T,N> Sqrt(const vtkm::Vec<T,N> &x) {
}
template<typename T>
static inline VTKM_EXEC_CONT
vtkm::Vec<T,4> Sqrt(const vtkm::Vec<T,4> &x) {
return vtkm::Vec<T,4>(vtkm::Sqrt(x[0]),
vtkm::Sqrt(x[1]),
vtkm::Sqrt(x[2]),
vtkm::Sqrt(x[3]));
vtkm::Vec<vtkm::FloatDefault,4> Sqrt(const vtkm::Vec<T,4> &x) {
return vtkm::Vec<vtkm::FloatDefault,4>(vtkm::Sqrt(x[0]),
vtkm::Sqrt(x[1]),
vtkm::Sqrt(x[2]),
vtkm::Sqrt(x[3]));
}
template<typename T>
static inline VTKM_EXEC_CONT
vtkm::Vec<T,3> Sqrt(const vtkm::Vec<T,3> &x) {
return vtkm::Vec<T,3>(vtkm::Sqrt(x[0]),
vtkm::Sqrt(x[1]),
vtkm::Sqrt(x[2]));
vtkm::Vec<vtkm::FloatDefault,3> Sqrt(const vtkm::Vec<T,3> &x) {
return vtkm::Vec<vtkm::FloatDefault,3>(vtkm::Sqrt(x[0]),
vtkm::Sqrt(x[1]),
vtkm::Sqrt(x[2]));
}
template<typename T>
static inline VTKM_EXEC_CONT
vtkm::Vec<T,2> Sqrt(const vtkm::Vec<T,2> &x) {
return vtkm::Vec<T,2>(vtkm::Sqrt(x[0]),
vtkm::Sqrt(x[1]));
vtkm::Vec<vtkm::FloatDefault,2> Sqrt(const vtkm::Vec<T,2> &x) {
return vtkm::Vec<vtkm::FloatDefault,2>(vtkm::Sqrt(x[0]),
vtkm::Sqrt(x[1]));
}
template<vtkm::IdComponent N>
static inline VTKM_EXEC_CONT
vtkm::Vec<vtkm::Float64,N> Sqrt(const vtkm::Vec<vtkm::Float64,N> &x) {
vtkm::Vec<vtkm::Float64,N> result;
for (vtkm::IdComponent index = 0; index < N; index++)
{
result[index] = vtkm::Sqrt(x[index]);
}
return result;
}
static inline VTKM_EXEC_CONT
vtkm::Vec<vtkm::Float64,4> Sqrt(const vtkm::Vec<vtkm::Float64,4> &x) {
return vtkm::Vec<vtkm::Float64,4>(vtkm::Sqrt(x[0]),
vtkm::Sqrt(x[1]),
vtkm::Sqrt(x[2]),
vtkm::Sqrt(x[3]));
}
static inline VTKM_EXEC_CONT
vtkm::Vec<vtkm::Float64,3> Sqrt(const vtkm::Vec<vtkm::Float64,3> &x) {
return vtkm::Vec<vtkm::Float64,3>(vtkm::Sqrt(x[0]),
vtkm::Sqrt(x[1]),
vtkm::Sqrt(x[2]));
}
static inline VTKM_EXEC_CONT
vtkm::Vec<vtkm::Float64,2> Sqrt(const vtkm::Vec<vtkm::Float64,2> &x) {
return vtkm::Vec<vtkm::Float64,2>(vtkm::Sqrt(x[0]),
vtkm::Sqrt(x[1]));
}
/// Compute the reciprocal square root of \p x. The result of this function is
......
......@@ -92,11 +92,18 @@ T MagnitudeTemplate(T x, vtkm::TypeTraitsScalarTag)
template<typename T>
VTKM_EXEC_CONT
typename vtkm::VecTraits<T>::ComponentType
MagnitudeTemplate(const T &x, vtkm::TypeTraitsVectorTag)
vtkm::FloatDefault MagnitudeTemplate(const T &x, vtkm::TypeTraitsVectorTag)
{
return vtkm::Sqrt(vtkm::MagnitudeSquared(x));
}
template<vtkm::IdComponent Size>
VTKM_EXEC_CONT
vtkm::Float64 MagnitudeTemplate(const vtkm::Vec<vtkm::Float64, Size> &x)
{
return vtkm::Sqrt(vtkm::MagnitudeSquared(x));
}
} // namespace detail
/// \brief Returns the magnitude of a vector.
......@@ -109,12 +116,17 @@ MagnitudeTemplate(const T &x, vtkm::TypeTraitsVectorTag)
///
template<typename T>
VTKM_EXEC_CONT
typename vtkm::VecTraits<T>::ComponentType
Magnitude(const T &x)
vtkm::FloatDefault Magnitude(const T &x)
{
return detail::MagnitudeTemplate(
x, typename vtkm::TypeTraits<T>::DimensionalityTag());
}
template<vtkm::IdComponent Size>
VTKM_EXEC_CONT
vtkm::Float64 Magnitude(const vtkm::Vec<vtkm::Float64, Size> &x)
{
return detail::MagnitudeTemplate(x);
}
// ----------------------------------------------------------------------------
namespace detail {
......
......@@ -40,6 +40,7 @@ set(headers
ResultDataSet.h
ResultField.h
Threshold.h
VectorMagnitude.h
VertexClustering.h
)
......@@ -58,6 +59,7 @@ set(header_template_sources
PointAverage.hxx
PointElevation.hxx
Threshold.hxx
VectorMagnitude.hxx
VertexClustering.hxx
)
......
//============================================================================
// 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.
//============================================================================
#ifndef vtk_m_filter_VectorMagnitude_h
#define vtk_m_filter_VectorMagnitude_h
#include <vtkm/filter/FilterField.h>
#include <vtkm/worklet/Magnitude.h>
namespace vtkm {
namespace filter {
class VectorMagnitude : public vtkm::filter::FilterField<VectorMagnitude>
{
public:
VTKM_CONT
VectorMagnitude();
template<typename T, typename StorageType, typename DerivedPolicy, typename DeviceAdapter>
VTKM_CONT
vtkm::filter::ResultField DoExecute(const vtkm::cont::DataSet &input,
const vtkm::cont::ArrayHandle<T, StorageType>& field,
const vtkm::filter::FieldMetadata& fieldMeta,
const vtkm::filter::PolicyBase<DerivedPolicy>& policy,
const DeviceAdapter& tag);
private:
vtkm::worklet::Magnitude Worklet;
};
}
} // namespace vtkm::filter
#include <vtkm/filter/VectorMagnitude.hxx>
#endif // vtk_m_filter_VectorMagnitude_h
//============================================================================
// 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/DispatcherMapField.h>
namespace vtkm {
namespace filter {
//-----------------------------------------------------------------------------
inline VTKM_CONT
VectorMagnitude::VectorMagnitude():
vtkm::filter::FilterField<VectorMagnitude>(),
Worklet()
{
this->SetOutputFieldName("magnitude");
}
//-----------------------------------------------------------------------------
template<typename T,
typename StorageType,
typename DerivedPolicy,
typename DeviceAdapter>
inline VTKM_CONT
vtkm::filter::ResultField
VectorMagnitude::DoExecute(const vtkm::cont::DataSet &inDataSet,
const vtkm::cont::ArrayHandle<T,StorageType> &field,
const vtkm::filter::FieldMetadata &fieldMetadata,
const vtkm::filter::PolicyBase<DerivedPolicy>&,
const DeviceAdapter&)
{
typedef typename vtkm::VecTraits<T>::ComponentType ComponentType;
vtkm::cont::ArrayHandle<ComponentType> outArray;
vtkm::worklet::DispatcherMapField<vtkm::worklet::Magnitude,
DeviceAdapter > dispatcher(this->Worklet);
dispatcher.Invoke(field, outArray);
return vtkm::filter::ResultField(inDataSet,
outArray,
this->GetOutputFieldName(),
fieldMetadata.GetAssociation(),
fieldMetadata.GetCellSetName());
}
}
} // namespace vtkm::filter
......@@ -30,6 +30,7 @@ set(unit_tests
UnitTestPointAverageFilter.cxx
UnitTestPointElevationFilter.cxx
UnitTestThresholdFilter.cxx
UnitTestVectorMagnitudeFilter.cxx
UnitTestVertexClusteringFilter.cxx
)
......
//============================================================================
// 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/filter/VectorMagnitude.h>
#include <vtkm/cont/testing/MakeTestDataSet.h>
#include <vtkm/cont/testing/Testing.h>
#include <vector>
namespace {
void TestVectorMagnitude()
{
std::cout << "Testing VectorMagnitude Filter" << 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(nVerts);
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);
vtkm::cont::DataSetFieldAdd::AddPointField(dataSet, "vec_pointvar", input);
vtkm::filter::ResultField result;
vtkm::filter::VectorMagnitude vm;
result = vm.Execute(dataSet, dataSet.GetField("vec_pointvar"));
VTKM_TEST_ASSERT( result.IsValid(), "result should be valid" );
VTKM_TEST_ASSERT(result.GetField().GetName() == "magnitude",
"Output field has wrong name.");
VTKM_TEST_ASSERT(result.GetField().GetAssociation() ==
vtkm::cont::Field::ASSOC_POINTS,
"Output field has wrong association");
vtkm::cont::ArrayHandle<vtkm::Float64> resultArrayHandle;
const bool valid = result.FieldAs(resultArrayHandle);
if(valid)
{
for (vtkm::Id i = 0; i < resultArrayHandle.GetNumberOfValues(); ++i)
{
VTKM_TEST_ASSERT(test_equal(std::sqrt(3*vars[i]*vars[i]),
resultArrayHandle.GetPortalConstControl().Get(i)),
"Wrong result for Magnitude worklet");
}
}
}
}
int UnitTestVectorMagnitudeFilter(int, char *[])
{
return vtkm::cont::testing::Testing::Run(TestVectorMagnitude);
}
......@@ -30,13 +30,13 @@ namespace worklet {
class Magnitude : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(FieldIn<Vec3>, FieldOut<>);
typedef void ControlSignature(FieldIn<VecAll>, FieldOut<Scalar>);
typedef void ExecutionSignature(_1, _2);
template<typename T>
template<typename T, typename T2>
VTKM_EXEC
void operator()(const vtkm::Vec<T,3> &inValue,
T &outValue) const
void operator()(const T &inValue,
T2 &outValue) const
{
outValue = vtkm::Magnitude(inValue);
}
......
......@@ -26,6 +26,7 @@ set(unit_tests
UnitTestExternalFaces.cxx
UnitTestFieldHistogram.cxx
UnitTestFieldStatistics.cxx
UnitTestMagnitude.cxx
UnitTestMarchingCubes.cxx
UnitTestPointElevation.cxx
UnitTestPointGradient.cxx
......
//============================================================================
// 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/Magnitude.h>
#include <vtkm/worklet/DispatcherMapField.h>
#include <vtkm/cont/testing/Testing.h>
namespace {
void TestMagnitude()
{
std::cout << "Testing Magnitude Worklet" << std::endl;
vtkm::worklet::Magnitude magnitudeWorklet;
typedef vtkm::cont::ArrayHandle<vtkm::Float64> ArrayReturnType;
typedef vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Int32, 4> > ArrayVectorType;
typedef ArrayVectorType::PortalControl PortalType;
ArrayVectorType pythagoreanTriples;
pythagoreanTriples.Allocate(5);
PortalType pt = pythagoreanTriples.GetPortalControl();
pt.Set(0, vtkm::make_Vec(3, 4, 5, 0));
pt.Set(1, vtkm::make_Vec(5, 12, 13, 0));
pt.Set(2, vtkm::make_Vec(8, 15, 17, 0));
pt.Set(3, vtkm::make_Vec(7, 24, 25, 0));
pt.Set(4, vtkm::make_Vec(9, 40, 41, 0));
vtkm::worklet::DispatcherMapField<vtkm::worklet::Magnitude>
dispatcher(magnitudeWorklet);
ArrayReturnType result;
dispatcher.Invoke(pythagoreanTriples,
result);
for (vtkm::Id i = 0; i < result.GetNumberOfValues(); ++i)
{
VTKM_TEST_ASSERT(
test_equal(std::sqrt(pt.Get(i)[0]*pt.Get(i)[0] +
pt.Get(i)[1]*pt.Get(i)[1] +
pt.Get(i)[2]*pt.Get(i)[2]),
result.GetPortalConstControl().Get(i)),
"Wrong result for Magnitude worklet");
}
}
}
int UnitTestMagnitude(int, char *[])
{
return vtkm::cont::testing::Testing::Run(TestMagnitude);
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment