mirror of
https://gitlab.kitware.com/vtk/vtk-m
synced 2024-09-16 17:22:55 +00:00
Allow Vector Gradient computation to output fancy fields and not gradient
Basically we can run the gradient worklet and compute the Divergence, Vorticity, and QCriterion without ever having to store the gradient tensor in global memory
This commit is contained in:
Robert Maynard
parent
85a62423bd
commit
8d04796b1f
@ -68,6 +68,13 @@ public:
|
||||
void SetComputeQCriterion(bool enable) { ComputeQCriterion = enable; }
|
||||
bool GetComputeQCriterion() const { return ComputeQCriterion; }
|
||||
|
||||
/// Add gradient field to the output data. The name of the array
|
||||
/// will be Gradients and will be a cell field unless \c ComputePointGradient
|
||||
/// is enabled. It is useful to turn this off when you are only interested
|
||||
/// in the results of Divergence, Vorticity, or QCriterion. The default is on.
|
||||
void SetComputeGradient(bool enable) { StoreGradient = enable; }
|
||||
bool GetComputeGradient() const { return StoreGradient; }
|
||||
|
||||
void SetDivergenceName(const std::string& name) { this->DivergenceName = name; }
|
||||
const std::string& GetDivergenceName() const { return this->DivergenceName; }
|
||||
|
||||
@ -89,7 +96,9 @@ private:
|
||||
bool ComputeDivergence;
|
||||
bool ComputeVorticity;
|
||||
bool ComputeQCriterion;
|
||||
bool StoreGradient;
|
||||
|
||||
std::string GradientsName;
|
||||
std::string DivergenceName;
|
||||
std::string VorticityName;
|
||||
std::string QCriterionName;
|
||||
|
||||
@ -23,7 +23,6 @@
|
||||
|
||||
namespace
|
||||
{
|
||||
|
||||
//-----------------------------------------------------------------------------
|
||||
template <typename HandleType>
|
||||
inline void add_field(vtkm::filter::ResultField& result,
|
||||
@ -43,51 +42,6 @@ inline void add_field(vtkm::filter::ResultField& result,
|
||||
}
|
||||
}
|
||||
|
||||
//-----------------------------------------------------------------------------
|
||||
template <typename T, typename S, typename DeviceAdapter>
|
||||
inline void add_extra_vec_fields(
|
||||
const vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Vec<T, 3>, 3>, S>& inField,
|
||||
const vtkm::filter::Gradient* const filter,
|
||||
vtkm::filter::ResultField& result,
|
||||
const DeviceAdapter&)
|
||||
{
|
||||
if (filter->GetComputeDivergence())
|
||||
{
|
||||
vtkm::cont::ArrayHandle<T> divergence;
|
||||
vtkm::worklet::DispatcherMapField<vtkm::worklet::Divergence, DeviceAdapter> dispatcher;
|
||||
dispatcher.Invoke(inField, divergence);
|
||||
|
||||
add_field(result, divergence, filter->GetDivergenceName());
|
||||
}
|
||||
|
||||
if (filter->GetComputeVorticity())
|
||||
{
|
||||
vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>> vorticity;
|
||||
vtkm::worklet::DispatcherMapField<vtkm::worklet::Vorticity, DeviceAdapter> dispatcher;
|
||||
dispatcher.Invoke(inField, vorticity);
|
||||
|
||||
add_field(result, vorticity, filter->GetVorticityName());
|
||||
}
|
||||
|
||||
if (filter->GetComputeQCriterion())
|
||||
{
|
||||
vtkm::cont::ArrayHandle<T> qc;
|
||||
vtkm::worklet::DispatcherMapField<vtkm::worklet::QCriterion, DeviceAdapter> dispatcher;
|
||||
dispatcher.Invoke(inField, qc);
|
||||
|
||||
add_field(result, qc, filter->GetQCriterionName());
|
||||
}
|
||||
}
|
||||
|
||||
template <typename T, typename S, typename DeviceAdapter>
|
||||
inline void add_extra_vec_fields(const vtkm::cont::ArrayHandle<T, S>&,
|
||||
const vtkm::filter::Gradient* const,
|
||||
vtkm::filter::ResultField&,
|
||||
const DeviceAdapter&)
|
||||
{
|
||||
//not a vector array handle so add nothing
|
||||
}
|
||||
|
||||
} //namespace
|
||||
|
||||
namespace vtkm
|
||||
@ -100,6 +54,8 @@ Gradient::Gradient()
|
||||
: ComputePointGradient(false)
|
||||
, ComputeVorticity(false)
|
||||
, ComputeQCriterion(false)
|
||||
, StoreGradient(true)
|
||||
, GradientsName("Gradients")
|
||||
, DivergenceName("Divergence")
|
||||
, VorticityName("Vorticity")
|
||||
, QCriterionName("QCriterion")
|
||||
@ -130,11 +86,16 @@ inline vtkm::filter::ResultField Gradient::DoExecute(
|
||||
std::string outputName = this->GetOutputFieldName();
|
||||
if (outputName.empty())
|
||||
{
|
||||
outputName = "Gradients";
|
||||
outputName = this->GradientsName;
|
||||
}
|
||||
|
||||
//todo: we need to ask the policy what storage type we should be using
|
||||
//If the input is implicit, we should know what to fall back to
|
||||
vtkm::worklet::GradientOutputFields<T> gradientfields(this->GetComputeGradient(),
|
||||
this->GetComputeDivergence(),
|
||||
this->GetComputeVorticity(),
|
||||
this->GetComputeQCriterion());
|
||||
|
||||
vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>> outArray;
|
||||
if (this->ComputePointGradient)
|
||||
{
|
||||
@ -142,6 +103,7 @@ inline vtkm::filter::ResultField Gradient::DoExecute(
|
||||
outArray = gradient.Run(vtkm::filter::ApplyPolicy(cells, policy),
|
||||
vtkm::filter::ApplyPolicy(coords, policy),
|
||||
inField,
|
||||
gradientfields,
|
||||
adapter);
|
||||
}
|
||||
else
|
||||
@ -150,17 +112,29 @@ inline vtkm::filter::ResultField Gradient::DoExecute(
|
||||
outArray = gradient.Run(vtkm::filter::ApplyPolicy(cells, policy),
|
||||
vtkm::filter::ApplyPolicy(coords, policy),
|
||||
inField,
|
||||
gradientfields,
|
||||
adapter);
|
||||
}
|
||||
VTKM_CONSTEXPR bool isVector = std::is_same<typename vtkm::VecTraits<T>::HasMultipleComponents,
|
||||
vtkm::VecTraitsTagMultipleComponents>::value;
|
||||
|
||||
vtkm::cont::Field::AssociationEnum fieldAssociation(this->ComputePointGradient
|
||||
? vtkm::cont::Field::ASSOC_POINTS
|
||||
: vtkm::cont::Field::ASSOC_CELL_SET);
|
||||
vtkm::filter::ResultField result(input, outArray, outputName, fieldAssociation, cells.GetName());
|
||||
|
||||
//Add the vorticity and qcriterion fields if they are enabled to the result
|
||||
add_extra_vec_fields(outArray, this, result, adapter);
|
||||
|
||||
if (this->GetComputeDivergence() && isVector)
|
||||
{
|
||||
add_field(result, gradientfields.Divergence, this->GetDivergenceName());
|
||||
}
|
||||
if (this->GetComputeVorticity() && isVector)
|
||||
{
|
||||
add_field(result, gradientfields.Vorticity, this->GetVorticityName());
|
||||
}
|
||||
if (this->GetComputeQCriterion() && isVector)
|
||||
{
|
||||
add_field(result, gradientfields.QCriterion, this->GetQCriterionName());
|
||||
}
|
||||
return result;
|
||||
}
|
||||
}
|
||||
|
||||
@ -66,6 +66,7 @@ set(headers
|
||||
#-----------------------------------------------------------------------------
|
||||
add_subdirectory(internal)
|
||||
add_subdirectory(contourtree)
|
||||
add_subdirectory(gradient)
|
||||
add_subdirectory(splatkernels)
|
||||
add_subdirectory(tetrahedralize)
|
||||
add_subdirectory(triangulate)
|
||||
|
||||
@ -21,192 +21,35 @@
|
||||
#ifndef vtk_m_worklet_Gradient_h
|
||||
#define vtk_m_worklet_Gradient_h
|
||||
|
||||
#include <vtkm/CellTraits.h>
|
||||
#include <vtkm/VecFromPortal.h>
|
||||
#include <vtkm/VecFromPortalPermute.h>
|
||||
#include <vtkm/exec/CellDerivative.h>
|
||||
#include <vtkm/exec/ExecutionWholeArray.h>
|
||||
#include <vtkm/exec/ParametricCoordinates.h>
|
||||
#include <vtkm/worklet/DispatcherMapTopology.h>
|
||||
#include <vtkm/worklet/WorkletMapField.h>
|
||||
#include <vtkm/worklet/WorkletMapTopology.h>
|
||||
|
||||
#include <vtkm/exec/CellDerivative.h>
|
||||
#include <vtkm/worklet/gradient/CellGradient.h>
|
||||
#include <vtkm/worklet/gradient/Divergence.h>
|
||||
#include <vtkm/worklet/gradient/GradientOutput.h>
|
||||
#include <vtkm/worklet/gradient/PointGradient.h>
|
||||
#include <vtkm/worklet/gradient/QCriterion.h>
|
||||
#include <vtkm/worklet/gradient/Vorticity.h>
|
||||
|
||||
namespace vtkm
|
||||
{
|
||||
namespace worklet
|
||||
{
|
||||
|
||||
template <typename T>
|
||||
struct GradientOutputFields;
|
||||
|
||||
namespace gradient
|
||||
{
|
||||
|
||||
template <typename T>
|
||||
struct GradientInType : vtkm::ListTagBase<T>
|
||||
{
|
||||
};
|
||||
|
||||
template <typename T>
|
||||
struct GradientOutType : vtkm::ListTagBase<vtkm::Vec<T, 3>>
|
||||
{
|
||||
};
|
||||
|
||||
struct GradientVecOutTypes : vtkm::ListTagBase<vtkm::Vec<vtkm::Vec<vtkm::Float32, 3>, 3>,
|
||||
vtkm::Vec<vtkm::Vec<vtkm::Float64, 3>, 3>>
|
||||
{
|
||||
};
|
||||
|
||||
template <typename T>
|
||||
struct CellGradient : vtkm::worklet::WorkletMapPointToCell
|
||||
{
|
||||
typedef void ControlSignature(CellSetIn,
|
||||
FieldInPoint<Vec3> pointCoordinates,
|
||||
FieldInPoint<GradientInType<T>> inputField,
|
||||
FieldOutCell<GradientOutType<T>> outputField);
|
||||
|
||||
typedef void ExecutionSignature(CellShape, PointCount, _2, _3, _4);
|
||||
typedef _1 InputDomain;
|
||||
|
||||
template <typename CellTagType,
|
||||
typename PointCoordVecType,
|
||||
typename FieldInVecType,
|
||||
typename FieldOutType>
|
||||
VTKM_EXEC void operator()(CellTagType shape,
|
||||
vtkm::IdComponent pointCount,
|
||||
const PointCoordVecType& wCoords,
|
||||
const FieldInVecType& field,
|
||||
FieldOutType& outputField) const
|
||||
{
|
||||
vtkm::Vec<vtkm::FloatDefault, 3> center =
|
||||
vtkm::exec::ParametricCoordinatesCenter(pointCount, shape, *this);
|
||||
|
||||
outputField = vtkm::exec::CellDerivative(field, wCoords, center, shape, *this);
|
||||
}
|
||||
};
|
||||
|
||||
template <typename T>
|
||||
struct PointGradient : public vtkm::worklet::WorkletMapCellToPoint
|
||||
{
|
||||
typedef void ControlSignature(CellSetIn,
|
||||
WholeCellSetIn<Point, Cell>,
|
||||
WholeArrayIn<Vec3> pointCoordinates,
|
||||
WholeArrayIn<GradientInType<T>> inputField,
|
||||
FieldOutPoint<GradientOutType<T>> 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
|
||||
{
|
||||
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 BaseGradientType = typename vtkm::BaseComponent<ValueType>::Type;
|
||||
const BaseGradientType invNumCells =
|
||||
static_cast<BaseGradientType>(1.) / static_cast<BaseGradientType>(numCells);
|
||||
using OutValueType = typename FieldOutType::ComponentType;
|
||||
outputField[0] = static_cast<OutValueType>(gradient[0] * invNumCells);
|
||||
outputField[1] = static_cast<OutValueType>(gradient[1] * invNumCells);
|
||||
outputField[2] = static_cast<OutValueType>(gradient[2] * invNumCells);
|
||||
}
|
||||
|
||||
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());
|
||||
}
|
||||
};
|
||||
|
||||
//-----------------------------------------------------------------------------
|
||||
template <typename CoordinateSystem, typename T, typename S, typename Device>
|
||||
struct DeducedPointGrad
|
||||
{
|
||||
DeducedPointGrad(const CoordinateSystem& coords,
|
||||
const vtkm::cont::ArrayHandle<T, S>& field,
|
||||
vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>>* result)
|
||||
GradientOutputFields<T>* result)
|
||||
: Points(&coords)
|
||||
, Field(&field)
|
||||
, Result(result)
|
||||
@ -226,7 +69,7 @@ struct DeducedPointGrad
|
||||
|
||||
const CoordinateSystem* const Points;
|
||||
const vtkm::cont::ArrayHandle<T, S>* const Field;
|
||||
vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>>* Result;
|
||||
GradientOutputFields<T>* Result;
|
||||
|
||||
private:
|
||||
void operator=(const DeducedPointGrad<CoordinateSystem, T, S, Device>&) = delete;
|
||||
@ -234,6 +77,94 @@ private:
|
||||
|
||||
} //namespace gradient
|
||||
|
||||
template <typename T>
|
||||
struct GradientOutputFields : public vtkm::exec::ExecutionObjectBase
|
||||
{
|
||||
|
||||
using ValueType = T;
|
||||
using BaseTType = typename vtkm::BaseComponent<T>::Type;
|
||||
|
||||
template <typename DeviceAdapter>
|
||||
struct ExecutionTypes
|
||||
{
|
||||
using Portal = vtkm::exec::GradientOutput<T, DeviceAdapter>;
|
||||
};
|
||||
|
||||
GradientOutputFields()
|
||||
: Gradient()
|
||||
, Divergence()
|
||||
, Vorticity()
|
||||
, QCriterion()
|
||||
, StoreGradient(true)
|
||||
, ComputeDivergence(false)
|
||||
, ComputeVorticity(false)
|
||||
, ComputeQCriterion(false)
|
||||
{
|
||||
}
|
||||
|
||||
GradientOutputFields(bool store, bool divergence, bool vorticity, bool qc)
|
||||
: Gradient()
|
||||
, Divergence()
|
||||
, Vorticity()
|
||||
, QCriterion()
|
||||
, StoreGradient(store)
|
||||
, ComputeDivergence(divergence)
|
||||
, ComputeVorticity(vorticity)
|
||||
, ComputeQCriterion(qc)
|
||||
{
|
||||
}
|
||||
|
||||
/// Add divergence field to the output data.
|
||||
/// The input array must have 3 components in order to compute this.
|
||||
/// The default is off.
|
||||
void SetComputeDivergence(bool enable) { ComputeDivergence = enable; }
|
||||
bool GetComputeDivergence() const { return ComputeDivergence; }
|
||||
|
||||
/// Add voriticity/curl field to the output data.
|
||||
/// The input array must have 3 components in order to compute this.
|
||||
/// The default is off.
|
||||
void SetComputeVorticity(bool enable) { ComputeVorticity = enable; }
|
||||
bool GetComputeVorticity() const { return ComputeVorticity; }
|
||||
|
||||
/// Add Q-criterion field to the output data.
|
||||
/// The input array must have 3 components in order to compute this.
|
||||
/// The default is off.
|
||||
void SetComputeQCriterion(bool enable) { ComputeQCriterion = enable; }
|
||||
bool GetComputeQCriterion() const { return ComputeQCriterion; }
|
||||
|
||||
/// Add gradient field to the output data.
|
||||
/// The input array must have 3 components in order to disable this.
|
||||
/// The default is on.
|
||||
void SetComputeGradient(bool enable) { StoreGradient = enable; }
|
||||
bool GetComputeGradient() const { return StoreGradient; }
|
||||
|
||||
//todo fix this for scalar
|
||||
template <typename DeviceAdapter>
|
||||
vtkm::exec::GradientOutput<T, DeviceAdapter> PrepareForOutput(vtkm::Id size, DeviceAdapter)
|
||||
{
|
||||
vtkm::exec::GradientOutput<T, DeviceAdapter> portal(this->StoreGradient,
|
||||
this->ComputeDivergence,
|
||||
this->ComputeVorticity,
|
||||
this->ComputeQCriterion,
|
||||
this->Gradient,
|
||||
this->Divergence,
|
||||
this->Vorticity,
|
||||
this->QCriterion,
|
||||
size);
|
||||
return portal;
|
||||
}
|
||||
|
||||
vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>> Gradient;
|
||||
vtkm::cont::ArrayHandle<BaseTType> Divergence;
|
||||
vtkm::cont::ArrayHandle<vtkm::Vec<BaseTType, 3>> Vorticity;
|
||||
vtkm::cont::ArrayHandle<BaseTType> QCriterion;
|
||||
|
||||
private:
|
||||
bool StoreGradient;
|
||||
bool ComputeDivergence;
|
||||
bool ComputeVorticity;
|
||||
bool ComputeQCriterion;
|
||||
};
|
||||
class PointGradient
|
||||
{
|
||||
public:
|
||||
@ -246,15 +177,29 @@ public:
|
||||
const CoordinateSystem& coords,
|
||||
const vtkm::cont::ArrayHandle<T, S>& field,
|
||||
DeviceAdapter device)
|
||||
{
|
||||
vtkm::worklet::GradientOutputFields<T> extraOutput(true, false, false, false);
|
||||
return this->Run(cells, coords, field, extraOutput, device);
|
||||
}
|
||||
|
||||
template <typename CellSetType,
|
||||
typename CoordinateSystem,
|
||||
typename T,
|
||||
typename S,
|
||||
typename DeviceAdapter>
|
||||
vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>> Run(const CellSetType& cells,
|
||||
const CoordinateSystem& coords,
|
||||
const vtkm::cont::ArrayHandle<T, S>& field,
|
||||
GradientOutputFields<T>& extraOutput,
|
||||
DeviceAdapter)
|
||||
{
|
||||
//we are using cast and call here as we pass the cells twice to the invoke
|
||||
//and want the type resolved once before hand instead of twice
|
||||
//by the dispatcher ( that will cost more in time and binary size )
|
||||
|
||||
vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>> output;
|
||||
gradient::DeducedPointGrad<CoordinateSystem, T, S, DeviceAdapter> func(coords, field, &output);
|
||||
gradient::DeducedPointGrad<CoordinateSystem, T, S, DeviceAdapter> func(
|
||||
coords, field, &extraOutput);
|
||||
vtkm::cont::CastAndCall(cells, func);
|
||||
return output;
|
||||
return extraOutput.Gradient;
|
||||
}
|
||||
};
|
||||
|
||||
@ -271,66 +216,30 @@ public:
|
||||
const vtkm::cont::ArrayHandle<T, S>& field,
|
||||
DeviceAdapter device)
|
||||
{
|
||||
vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>> output;
|
||||
vtkm::worklet::DispatcherMapTopology<vtkm::worklet::gradient::CellGradient<T>, DeviceAdapter>
|
||||
dispatcher;
|
||||
dispatcher.Invoke(cells, coords, field, output);
|
||||
return output;
|
||||
vtkm::worklet::GradientOutputFields<T> extra(true, false, false, false);
|
||||
return this->Run(cells, coords, field, extra, device);
|
||||
}
|
||||
};
|
||||
|
||||
struct Divergence : public vtkm::worklet::WorkletMapField
|
||||
{
|
||||
typedef void ControlSignature(FieldIn<gradient::GradientVecOutTypes> input,
|
||||
FieldOut<Scalar> output);
|
||||
typedef void ExecutionSignature(_1, _2);
|
||||
typedef _1 InputDomain;
|
||||
|
||||
template <typename InputType, typename OutputType>
|
||||
VTKM_EXEC void operator()(const InputType& input, OutputType& divergence) const
|
||||
template <typename CellSetType,
|
||||
typename CoordinateSystem,
|
||||
typename T,
|
||||
typename S,
|
||||
typename DeviceAdapter>
|
||||
vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>> Run(const CellSetType& cells,
|
||||
const CoordinateSystem& coords,
|
||||
const vtkm::cont::ArrayHandle<T, S>& field,
|
||||
GradientOutputFields<T>& extraOutput,
|
||||
DeviceAdapter)
|
||||
{
|
||||
divergence = input[0][0] + input[1][1] + input[2][2];
|
||||
}
|
||||
};
|
||||
using DispatcherType =
|
||||
vtkm::worklet::DispatcherMapTopology<vtkm::worklet::gradient::CellGradient<T>, DeviceAdapter>;
|
||||
|
||||
struct Vorticity : public vtkm::worklet::WorkletMapField
|
||||
{
|
||||
typedef void ControlSignature(FieldIn<gradient::GradientVecOutTypes> input,
|
||||
FieldOut<Vec3> output);
|
||||
typedef void ExecutionSignature(_1, _2);
|
||||
typedef _1 InputDomain;
|
||||
vtkm::worklet::gradient::CellGradient<T> worklet;
|
||||
DispatcherType dispatcher(worklet);
|
||||
|
||||
template <typename InputType, typename OutputType>
|
||||
VTKM_EXEC void operator()(const InputType& input, OutputType& vorticity) const
|
||||
{
|
||||
vorticity[0] = input[2][1] - input[1][2];
|
||||
vorticity[1] = input[0][2] - input[2][0];
|
||||
vorticity[2] = input[1][0] - input[0][1];
|
||||
}
|
||||
};
|
||||
|
||||
struct QCriterion : public vtkm::worklet::WorkletMapField
|
||||
{
|
||||
typedef void ControlSignature(FieldIn<gradient::GradientVecOutTypes> input,
|
||||
FieldOut<Scalar> output);
|
||||
typedef void ExecutionSignature(_1, _2);
|
||||
typedef _1 InputDomain;
|
||||
|
||||
template <typename InputType, typename OutputType>
|
||||
VTKM_EXEC void operator()(const InputType& input, OutputType& qcriterion) const
|
||||
{
|
||||
OutputType t1 = ((input[2][1] - input[1][2]) * (input[2][1] - input[1][2]) +
|
||||
(input[1][0] - input[0][1]) * (input[1][0] - input[0][1]) +
|
||||
(input[0][2] - input[2][0]) * (input[0][2] - input[2][0])) /
|
||||
2.0f;
|
||||
OutputType t2 = input[0][0] * input[0][0] + input[1][1] * input[1][1] +
|
||||
input[2][2] * input[2][2] +
|
||||
((input[1][0] + input[0][1]) * (input[1][0] + input[0][1]) +
|
||||
(input[2][0] + input[0][2]) * (input[2][0] + input[0][2]) +
|
||||
(input[2][1] + input[1][2]) * (input[2][1] + input[1][2])) /
|
||||
2.0f;
|
||||
|
||||
qcriterion = (t1 - t2) / 2.0f;
|
||||
dispatcher.Invoke(cells, coords, field, extraOutput);
|
||||
return extraOutput.Gradient;
|
||||
}
|
||||
};
|
||||
}
|
||||
|
||||
31
vtkm/worklet/gradient/CMakeLists.txt
Normal file
31
vtkm/worklet/gradient/CMakeLists.txt
Normal file
@ -0,0 +1,31 @@
|
||||
##============================================================================
|
||||
## 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 2016 Sandia Corporation.
|
||||
## Copyright 2016 UT-Battelle, LLC.
|
||||
## Copyright 2016 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.
|
||||
##============================================================================
|
||||
|
||||
set(headers
|
||||
CellGradient.h
|
||||
Divergence.h
|
||||
GradientOutput.h
|
||||
PointGradient.h
|
||||
QCriterion.h
|
||||
Vorticity.h
|
||||
)
|
||||
|
||||
#-----------------------------------------------------------------------------
|
||||
vtkm_declare_headers(${headers})
|
||||
73
vtkm/worklet/gradient/CellGradient.h
Normal file
73
vtkm/worklet/gradient/CellGradient.h
Normal file
@ -0,0 +1,73 @@
|
||||
//============================================================================
|
||||
// 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_worklet_gradient_CellGradient_h
|
||||
#define vtk_m_worklet_gradient_CellGradient_h
|
||||
|
||||
#include <vtkm/exec/CellDerivative.h>
|
||||
#include <vtkm/exec/ParametricCoordinates.h>
|
||||
#include <vtkm/worklet/WorkletMapTopology.h>
|
||||
|
||||
#include <vtkm/worklet/gradient/GradientOutput.h>
|
||||
|
||||
namespace vtkm
|
||||
{
|
||||
namespace worklet
|
||||
{
|
||||
namespace gradient
|
||||
{
|
||||
|
||||
template <typename T>
|
||||
struct CellGradientInType : vtkm::ListTagBase<T>
|
||||
{
|
||||
};
|
||||
|
||||
template <typename T>
|
||||
struct CellGradient : vtkm::worklet::WorkletMapPointToCell
|
||||
{
|
||||
typedef void ControlSignature(CellSetIn,
|
||||
FieldInPoint<Vec3> pointCoordinates,
|
||||
FieldInPoint<CellGradientInType<T>> inputField,
|
||||
GradientOutputs outputFields);
|
||||
|
||||
typedef void ExecutionSignature(CellShape, PointCount, _2, _3, _4);
|
||||
typedef _1 InputDomain;
|
||||
|
||||
template <typename CellTagType,
|
||||
typename PointCoordVecType,
|
||||
typename FieldInVecType,
|
||||
typename GradientOutType>
|
||||
VTKM_EXEC void operator()(CellTagType shape,
|
||||
vtkm::IdComponent pointCount,
|
||||
const PointCoordVecType& wCoords,
|
||||
const FieldInVecType& field,
|
||||
GradientOutType& outputGradient) const
|
||||
{
|
||||
vtkm::Vec<vtkm::FloatDefault, 3> center =
|
||||
vtkm::exec::ParametricCoordinatesCenter(pointCount, shape, *this);
|
||||
|
||||
outputGradient = vtkm::exec::CellDerivative(field, wCoords, center, shape, *this);
|
||||
}
|
||||
};
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#endif
|
||||
54
vtkm/worklet/gradient/Divergence.h
Normal file
54
vtkm/worklet/gradient/Divergence.h
Normal file
@ -0,0 +1,54 @@
|
||||
//============================================================================
|
||||
// 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_worklet_gradient_Divergence_h
|
||||
#define vtk_m_worklet_gradient_Divergence_h
|
||||
|
||||
#include <vtkm/worklet/WorkletMapField.h>
|
||||
|
||||
namespace vtkm
|
||||
{
|
||||
namespace worklet
|
||||
{
|
||||
namespace gradient
|
||||
{
|
||||
|
||||
struct DivergenceTypes : vtkm::ListTagBase<vtkm::Vec<vtkm::Vec<vtkm::Float32, 3>, 3>,
|
||||
vtkm::Vec<vtkm::Vec<vtkm::Float64, 3>, 3>>
|
||||
{
|
||||
};
|
||||
|
||||
|
||||
struct Divergence : public vtkm::worklet::WorkletMapField
|
||||
{
|
||||
typedef void ControlSignature(FieldIn<DivergenceTypes> input, FieldOut<Scalar> output);
|
||||
typedef void ExecutionSignature(_1, _2);
|
||||
typedef _1 InputDomain;
|
||||
|
||||
template <typename InputType, typename OutputType>
|
||||
VTKM_EXEC void operator()(const InputType& input, OutputType& divergence) const
|
||||
{
|
||||
divergence = input[0][0] + input[1][1] + input[2][2];
|
||||
}
|
||||
};
|
||||
}
|
||||
}
|
||||
}
|
||||
#endif
|
||||
275
vtkm/worklet/gradient/GradientOutput.h
Normal file
275
vtkm/worklet/gradient/GradientOutput.h
Normal file
@ -0,0 +1,275 @@
|
||||
//============================================================================
|
||||
// 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_worklet_gradient_GradientOutput_h
|
||||
#define vtk_m_worklet_gradient_GradientOutput_h
|
||||
|
||||
#include <vtkm/BaseComponent.h>
|
||||
|
||||
#include <vtkm/cont/arg/TransportTagArrayOut.h>
|
||||
#include <vtkm/cont/arg/TransportTagExecObject.h>
|
||||
|
||||
#include <vtkm/exec/ExecutionObjectBase.h>
|
||||
#include <vtkm/exec/arg/FetchTagArrayDirectOut.h>
|
||||
|
||||
#include <vtkm/worklet/gradient/Divergence.h>
|
||||
#include <vtkm/worklet/gradient/QCriterion.h>
|
||||
#include <vtkm/worklet/gradient/Vorticity.h>
|
||||
|
||||
namespace vtkm
|
||||
{
|
||||
namespace exec
|
||||
{
|
||||
|
||||
template <typename T, typename DeviceAdapter>
|
||||
struct GradientScalarOutput : public vtkm::exec::ExecutionObjectBase
|
||||
{
|
||||
using ValueType = vtkm::Vec<T, 3>;
|
||||
using BaseTType = typename vtkm::BaseComponent<T>::Type;
|
||||
|
||||
struct PortalTypes
|
||||
{
|
||||
using HandleType = vtkm::cont::ArrayHandle<ValueType>;
|
||||
using ExecutionTypes = typename HandleType::template ExecutionTypes<DeviceAdapter>;
|
||||
using Portal = typename ExecutionTypes::Portal;
|
||||
};
|
||||
|
||||
GradientScalarOutput() = default;
|
||||
|
||||
GradientScalarOutput(bool,
|
||||
bool,
|
||||
bool,
|
||||
bool,
|
||||
vtkm::cont::ArrayHandle<ValueType>& gradient,
|
||||
vtkm::cont::ArrayHandle<BaseTType>&,
|
||||
vtkm::cont::ArrayHandle<vtkm::Vec<BaseTType, 3>>&,
|
||||
vtkm::cont::ArrayHandle<BaseTType>&,
|
||||
vtkm::Id size)
|
||||
{
|
||||
this->GradientPortal = gradient.PrepareForOutput(size, DeviceAdapter());
|
||||
}
|
||||
|
||||
VTKM_SUPPRESS_EXEC_WARNINGS
|
||||
VTKM_EXEC
|
||||
void Set(vtkm::Id index, const vtkm::Vec<T, 3>& value) const
|
||||
{
|
||||
this->GradientPortal.Set(index, value);
|
||||
}
|
||||
|
||||
typename PortalTypes::Portal GradientPortal;
|
||||
};
|
||||
|
||||
template <typename T, typename DeviceAdapter>
|
||||
struct GradientVecOutput : public vtkm::exec::ExecutionObjectBase
|
||||
{
|
||||
using ValueType = vtkm::Vec<T, 3>;
|
||||
using BaseTType = typename vtkm::BaseComponent<T>::Type;
|
||||
|
||||
template <typename FieldType>
|
||||
struct PortalTypes
|
||||
{
|
||||
using HandleType = vtkm::cont::ArrayHandle<FieldType>;
|
||||
using ExecutionTypes = typename HandleType::template ExecutionTypes<DeviceAdapter>;
|
||||
using Portal = typename ExecutionTypes::Portal;
|
||||
};
|
||||
|
||||
GradientVecOutput() = default;
|
||||
|
||||
GradientVecOutput(bool g,
|
||||
bool d,
|
||||
bool v,
|
||||
bool q,
|
||||
vtkm::cont::ArrayHandle<ValueType>& gradient,
|
||||
vtkm::cont::ArrayHandle<BaseTType>& divergence,
|
||||
vtkm::cont::ArrayHandle<vtkm::Vec<BaseTType, 3>>& vorticity,
|
||||
vtkm::cont::ArrayHandle<BaseTType>& qcriterion,
|
||||
vtkm::Id size)
|
||||
{
|
||||
this->SetGradient = g;
|
||||
this->SetDivergence = d;
|
||||
this->SetVorticity = v;
|
||||
this->SetQCriterion = q;
|
||||
|
||||
DeviceAdapter device;
|
||||
if (g)
|
||||
{
|
||||
this->GradientPortal = gradient.PrepareForOutput(size, device);
|
||||
}
|
||||
if (d)
|
||||
{
|
||||
this->DivergencePortal = divergence.PrepareForOutput(size, device);
|
||||
}
|
||||
if (v)
|
||||
{
|
||||
this->VorticityPortal = vorticity.PrepareForOutput(size, device);
|
||||
}
|
||||
if (q)
|
||||
{
|
||||
this->QCriterionPortal = qcriterion.PrepareForOutput(size, device);
|
||||
}
|
||||
}
|
||||
|
||||
VTKM_SUPPRESS_EXEC_WARNINGS
|
||||
VTKM_EXEC
|
||||
void Set(vtkm::Id index, const vtkm::Vec<T, 3>& value) const
|
||||
{
|
||||
if (this->SetGradient)
|
||||
{
|
||||
this->GradientPortal.Set(index, value);
|
||||
}
|
||||
if (this->SetDivergence)
|
||||
{
|
||||
vtkm::worklet::gradient::Divergence divergence;
|
||||
BaseTType output;
|
||||
divergence(value, output);
|
||||
this->DivergencePortal.Set(index, output);
|
||||
}
|
||||
if (this->SetVorticity)
|
||||
{
|
||||
vtkm::worklet::gradient::Vorticity vorticity;
|
||||
T output;
|
||||
vorticity(value, output);
|
||||
this->VorticityPortal.Set(index, output);
|
||||
}
|
||||
if (this->SetQCriterion)
|
||||
{
|
||||
vtkm::worklet::gradient::QCriterion qc;
|
||||
BaseTType output;
|
||||
qc(value, output);
|
||||
this->QCriterionPortal.Set(index, output);
|
||||
}
|
||||
}
|
||||
|
||||
bool SetGradient;
|
||||
bool SetDivergence;
|
||||
bool SetVorticity;
|
||||
bool SetQCriterion;
|
||||
|
||||
typename PortalTypes<ValueType>::Portal GradientPortal;
|
||||
typename PortalTypes<BaseTType>::Portal DivergencePortal;
|
||||
typename PortalTypes<vtkm::Vec<BaseTType, 3>>::Portal VorticityPortal;
|
||||
typename PortalTypes<BaseTType>::Portal QCriterionPortal;
|
||||
};
|
||||
|
||||
template <typename T, typename DeviceAdapter>
|
||||
struct GradientOutput : public GradientScalarOutput<T, DeviceAdapter>
|
||||
{
|
||||
#if defined(VTKM_MSVC) && (_MSC_VER == 1800) // workaround for VS2013
|
||||
template <typename... Args>
|
||||
GradientOutput(Args&&... args)
|
||||
: GradientScalarOutput<T, DeviceAdapter>::GradientScalarOutput(std::forward<Args>(args)...)
|
||||
{
|
||||
}
|
||||
#else
|
||||
using GradientScalarOutput<T, DeviceAdapter>::GradientScalarOutput;
|
||||
#endif
|
||||
};
|
||||
|
||||
template <typename DeviceAdapter>
|
||||
struct GradientOutput<vtkm::Vec<vtkm::Float32, 3>, DeviceAdapter>
|
||||
: public GradientVecOutput<vtkm::Vec<vtkm::Float32, 3>, DeviceAdapter>
|
||||
{
|
||||
|
||||
#if defined(VTKM_MSVC) && (_MSC_VER == 1800) // workaround for VS2013
|
||||
template <typename... Args>
|
||||
GradientOutput(Args&&... args)
|
||||
: GradientVecOutput<vtkm::Vec<vtkm::Float32, 3>, DeviceAdapter>::GradientVecOutput(
|
||||
std::forward<Args>(args)...)
|
||||
{
|
||||
}
|
||||
#else
|
||||
using GradientVecOutput<vtkm::Vec<vtkm::Float32, 3>, DeviceAdapter>::GradientVecOutput;
|
||||
#endif
|
||||
};
|
||||
|
||||
template <typename DeviceAdapter>
|
||||
struct GradientOutput<vtkm::Vec<vtkm::Float64, 3>, DeviceAdapter>
|
||||
: public GradientVecOutput<vtkm::Vec<vtkm::Float64, 3>, DeviceAdapter>
|
||||
{
|
||||
|
||||
#if defined(VTKM_MSVC) && (_MSC_VER == 1800) // workaround for VS2013
|
||||
template <typename... Args>
|
||||
GradientOutput(Args&&... args)
|
||||
: GradientVecOutput<vtkm::Vec<vtkm::Float64, 3>, DeviceAdapter>::GradientVecOutput(
|
||||
std::forward<Args>(args)...)
|
||||
{
|
||||
}
|
||||
#else
|
||||
using GradientVecOutput<vtkm::Vec<vtkm::Float64, 3>, DeviceAdapter>::GradientVecOutput;
|
||||
#endif
|
||||
};
|
||||
}
|
||||
} // namespace vtkm::exec
|
||||
|
||||
|
||||
namespace vtkm
|
||||
{
|
||||
namespace cont
|
||||
{
|
||||
namespace arg
|
||||
{
|
||||
|
||||
/// \brief \c Transport tag for output arrays.
|
||||
///
|
||||
/// \c TransportTagArrayOut is a tag used with the \c Transport class to
|
||||
/// transport \c ArrayHandle objects for output data.
|
||||
///
|
||||
struct TransportTagGradientOut
|
||||
{
|
||||
};
|
||||
|
||||
template <typename ContObjectType, typename Device>
|
||||
struct Transport<vtkm::cont::arg::TransportTagGradientOut, ContObjectType, Device>
|
||||
{
|
||||
using ExecObjectType = vtkm::exec::GradientOutput<typename ContObjectType::ValueType, Device>;
|
||||
|
||||
template <typename InputDomainType>
|
||||
VTKM_CONT ExecObjectType operator()(ContObjectType object,
|
||||
const InputDomainType& vtkmNotUsed(inputDomain),
|
||||
vtkm::Id vtkmNotUsed(inputRange),
|
||||
vtkm::Id outputRange) const
|
||||
{
|
||||
return object.PrepareForOutput(outputRange, Device());
|
||||
}
|
||||
};
|
||||
}
|
||||
}
|
||||
} // namespace vtkm::cont::arg
|
||||
|
||||
|
||||
namespace vtkm
|
||||
{
|
||||
namespace worklet
|
||||
{
|
||||
namespace gradient
|
||||
{
|
||||
|
||||
|
||||
struct GradientOutputs : vtkm::cont::arg::ControlSignatureTagBase
|
||||
{
|
||||
using TypeCheckTag = vtkm::cont::arg::TypeCheckTagExecObject;
|
||||
using TransportTag = vtkm::cont::arg::TransportTagGradientOut;
|
||||
using FetchTag = vtkm::exec::arg::FetchTagArrayDirectOut;
|
||||
};
|
||||
}
|
||||
}
|
||||
} // namespace vtkm::worklet::gradient
|
||||
|
||||
#endif
|
||||
164
vtkm/worklet/gradient/PointGradient.h
Normal file
164
vtkm/worklet/gradient/PointGradient.h
Normal file
@ -0,0 +1,164 @@
|
||||
//============================================================================
|
||||
// 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_worklet_gradient_PointGradient_h
|
||||
#define vtk_m_worklet_gradient_PointGradient_h
|
||||
|
||||
#include <vtkm/exec/CellDerivative.h>
|
||||
#include <vtkm/exec/ParametricCoordinates.h>
|
||||
#include <vtkm/worklet/WorkletMapTopology.h>
|
||||
|
||||
#include <vtkm/worklet/gradient/GradientOutput.h>
|
||||
|
||||
|
||||
namespace vtkm
|
||||
{
|
||||
namespace worklet
|
||||
{
|
||||
namespace gradient
|
||||
{
|
||||
|
||||
template <typename T>
|
||||
struct PointGradientInType : vtkm::ListTagBase<T>
|
||||
{
|
||||
};
|
||||
|
||||
template <typename T>
|
||||
struct PointGradient : public vtkm::worklet::WorkletMapCellToPoint
|
||||
{
|
||||
typedef void ControlSignature(CellSetIn,
|
||||
WholeCellSetIn<Point, Cell>,
|
||||
WholeArrayIn<Vec3> pointCoordinates,
|
||||
WholeArrayIn<PointGradientInType<T>> inputField,
|
||||
GradientOutputs outputFields);
|
||||
|
||||
typedef void ExecutionSignature(CellCount, CellIndices, WorkIndex, _2, _3, _4, _5);
|
||||
typedef _1 InputDomain;
|
||||
|
||||
template <typename FromIndexType,
|
||||
typename CellSetInType,
|
||||
typename WholeCoordinatesIn,
|
||||
typename WholeFieldIn,
|
||||
typename GradientOutType>
|
||||
VTKM_EXEC void operator()(const vtkm::IdComponent& numCells,
|
||||
const FromIndexType& cellIds,
|
||||
const vtkm::Id& pointId,
|
||||
const CellSetInType& geometry,
|
||||
const WholeCoordinatesIn& pointCoordinates,
|
||||
const WholeFieldIn& inputField,
|
||||
GradientOutType& outputGradient) 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 BaseGradientType = typename vtkm::BaseComponent<ValueType>::Type;
|
||||
const BaseGradientType invNumCells =
|
||||
static_cast<BaseGradientType>(1.) / static_cast<BaseGradientType>(numCells);
|
||||
|
||||
gradient[0] = gradient[0] * invNumCells;
|
||||
gradient[1] = gradient[1] * invNumCells;
|
||||
gradient[2] = gradient[2] * invNumCells;
|
||||
outputGradient = gradient;
|
||||
}
|
||||
|
||||
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());
|
||||
}
|
||||
};
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#endif
|
||||
60
vtkm/worklet/gradient/QCriterion.h
Normal file
60
vtkm/worklet/gradient/QCriterion.h
Normal file
@ -0,0 +1,60 @@
|
||||
//============================================================================
|
||||
// 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_worklet_gradient_QCriterion_h
|
||||
#define vtk_m_worklet_gradient_QCriterion_h
|
||||
|
||||
#include <vtkm/worklet/WorkletMapField.h>
|
||||
|
||||
namespace vtkm
|
||||
{
|
||||
namespace worklet
|
||||
{
|
||||
namespace gradient
|
||||
{
|
||||
|
||||
struct QCriterionTypes : vtkm::ListTagBase<vtkm::Vec<vtkm::Vec<vtkm::Float32, 3>, 3>,
|
||||
vtkm::Vec<vtkm::Vec<vtkm::Float64, 3>, 3>>
|
||||
{
|
||||
};
|
||||
|
||||
struct QCriterion : public vtkm::worklet::WorkletMapField
|
||||
{
|
||||
typedef void ControlSignature(FieldIn<QCriterionTypes> input, FieldOut<Scalar> output);
|
||||
typedef void ExecutionSignature(_1, _2);
|
||||
typedef _1 InputDomain;
|
||||
|
||||
template <typename InputType, typename OutputType>
|
||||
VTKM_EXEC void operator()(const InputType& input, OutputType& qcriterion) const
|
||||
{
|
||||
const vtkm::Vec<OutputType, 3> v(
|
||||
input[2][1] - input[1][2], input[0][2] - input[2][0], input[1][0] - input[0][1]);
|
||||
const vtkm::Vec<OutputType, 3> s(
|
||||
input[2][1] + input[1][2], input[0][2] + input[2][0], input[1][0] + input[0][1]);
|
||||
const vtkm::Vec<OutputType, 3> d(input[0][0], input[1][1], input[2][2]);
|
||||
|
||||
//compute QCriterion
|
||||
qcriterion = ((vtkm::dot(v, v) / 2.0f) - (vtkm::dot(d, d) + (vtkm::dot(s, s) / 2.0f))) / 2.0f;
|
||||
}
|
||||
};
|
||||
}
|
||||
}
|
||||
}
|
||||
#endif
|
||||
57
vtkm/worklet/gradient/Vorticity.h
Normal file
57
vtkm/worklet/gradient/Vorticity.h
Normal file
@ -0,0 +1,57 @@
|
||||
//============================================================================
|
||||
// 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_worklet_gradient_Vorticity_h
|
||||
#define vtk_m_worklet_gradient_Vorticity_h
|
||||
|
||||
#include <vtkm/worklet/WorkletMapField.h>
|
||||
|
||||
namespace vtkm
|
||||
{
|
||||
namespace worklet
|
||||
{
|
||||
namespace gradient
|
||||
{
|
||||
|
||||
struct VorticityTypes : vtkm::ListTagBase<vtkm::Vec<vtkm::Vec<vtkm::Float32, 3>, 3>,
|
||||
vtkm::Vec<vtkm::Vec<vtkm::Float64, 3>, 3>>
|
||||
{
|
||||
};
|
||||
|
||||
|
||||
struct Vorticity : public vtkm::worklet::WorkletMapField
|
||||
{
|
||||
typedef void ControlSignature(FieldIn<VorticityTypes> input, FieldOut<Vec3> output);
|
||||
typedef void ExecutionSignature(_1, _2);
|
||||
typedef _1 InputDomain;
|
||||
|
||||
template <typename InputType, typename OutputType>
|
||||
VTKM_EXEC void operator()(const InputType& input, OutputType& vorticity) const
|
||||
{
|
||||
vorticity[0] = input[2][1] - input[1][2];
|
||||
vorticity[1] = input[0][2] - input[2][0];
|
||||
vorticity[2] = input[1][0] - input[0][1];
|
||||
}
|
||||
};
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#endif
|
||||
@ -69,10 +69,10 @@ void TestCellGradientUniform3D()
|
||||
gradient.Run(dataSet.GetCellSet(), dataSet.GetCoordinateSystem(), input, DeviceAdapter());
|
||||
|
||||
vtkm::Vec<vtkm::Float32, 3> expected[4] = {
|
||||
{ 10.025, 30.075, 60.125 },
|
||||
{ 10.025, 30.075, 60.125 },
|
||||
{ 10.025, 30.075, 60.175 },
|
||||
{ 10.025, 30.075, 60.175 },
|
||||
{ 10.025f, 30.075f, 60.125f },
|
||||
{ 10.025f, 30.075f, 60.125f },
|
||||
{ 10.025f, 30.075f, 60.175f },
|
||||
{ 10.025f, 30.075f, 60.175f },
|
||||
};
|
||||
for (int i = 0; i < 4; ++i)
|
||||
{
|
||||
@ -84,8 +84,9 @@ void TestCellGradientUniform3D()
|
||||
template <typename DeviceAdapter>
|
||||
void TestCellGradientUniform3DWithVectorField()
|
||||
{
|
||||
std::cout << "Testing CellGradient Worklet with a vector field on 3D strucutred data"
|
||||
<< std::endl;
|
||||
std::cout
|
||||
<< "Testing CellGradient and QCriterion Worklet with a vector field on 3D strucutred data"
|
||||
<< std::endl;
|
||||
vtkm::cont::testing::MakeTestDataSet testDataSet;
|
||||
vtkm::cont::DataSet dataSet = testDataSet.Make3DUniformDataSet0();
|
||||
|
||||
@ -102,9 +103,24 @@ void TestCellGradientUniform3DWithVectorField()
|
||||
|
||||
//we need to add Vec3 array to the dataset
|
||||
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Vec<vtkm::Float64, 3>, 3>> result;
|
||||
|
||||
vtkm::worklet::GradientOutputFields<vtkm::Vec<vtkm::Float64, 3>> extraOutput;
|
||||
extraOutput.SetComputeDivergence(false);
|
||||
extraOutput.SetComputeVorticity(false);
|
||||
extraOutput.SetComputeQCriterion(true);
|
||||
|
||||
vtkm::worklet::CellGradient gradient;
|
||||
result =
|
||||
gradient.Run(dataSet.GetCellSet(), dataSet.GetCoordinateSystem(), input, DeviceAdapter());
|
||||
result = gradient.Run(
|
||||
dataSet.GetCellSet(), dataSet.GetCoordinateSystem(), input, extraOutput, DeviceAdapter());
|
||||
|
||||
VTKM_TEST_ASSERT((extraOutput.Gradient.GetNumberOfValues() == 4),
|
||||
"Gradient field should be generated");
|
||||
VTKM_TEST_ASSERT((extraOutput.Divergence.GetNumberOfValues() == 0),
|
||||
"Divergence field shouldn't be generated");
|
||||
VTKM_TEST_ASSERT((extraOutput.Vorticity.GetNumberOfValues() == 0),
|
||||
"Vorticity field shouldn't be generated");
|
||||
VTKM_TEST_ASSERT((extraOutput.QCriterion.GetNumberOfValues() == 4),
|
||||
"QCriterion field should be generated");
|
||||
|
||||
vtkm::Vec<vtkm::Vec<vtkm::Float64, 3>, 3> expected[4] = {
|
||||
{ { 10.025, 10.025, 10.025 }, { 30.075, 30.075, 30.075 }, { 60.125, 60.125, 60.125 } },
|
||||
@ -123,6 +139,83 @@ void TestCellGradientUniform3DWithVectorField()
|
||||
"Wrong result for vec field CellGradient worklet on 3D uniform data");
|
||||
VTKM_TEST_ASSERT(test_equal(e[2], r[2]),
|
||||
"Wrong result for vec field CellGradient worklet on 3D uniform data");
|
||||
|
||||
const vtkm::Vec<vtkm::Float64, 3> v(e[2][1] - e[1][2], e[0][2] - e[2][0], e[1][0] - e[0][1]);
|
||||
const vtkm::Vec<vtkm::Float64, 3> s(e[2][1] + e[1][2], e[0][2] + e[2][0], e[1][0] + e[0][1]);
|
||||
const vtkm::Vec<vtkm::Float64, 3> d(e[0][0], e[1][1], e[2][2]);
|
||||
|
||||
//compute QCriterion
|
||||
vtkm::Float64 qcriterion =
|
||||
((vtkm::dot(v, v) / 2.0f) - (vtkm::dot(d, d) + (vtkm::dot(s, s) / 2.0f))) / 2.0f;
|
||||
|
||||
vtkm::Float64 q = extraOutput.QCriterion.GetPortalConstControl().Get(i);
|
||||
|
||||
VTKM_TEST_ASSERT(
|
||||
test_equal(qcriterion, q),
|
||||
"Wrong result for QCriterion field of CellGradient worklet on 3D uniform data");
|
||||
}
|
||||
}
|
||||
|
||||
template <typename DeviceAdapter>
|
||||
void TestCellGradientUniform3DWithVectorField2()
|
||||
{
|
||||
std::cout << "Testing CellGradient Worklet with a vector field on 3D strucutred data" << std::endl
|
||||
<< "Disabling Gradient computation and enabling Divergence, and Vorticity" << 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);
|
||||
|
||||
vtkm::worklet::GradientOutputFields<vtkm::Vec<vtkm::Float64, 3>> extraOutput;
|
||||
extraOutput.SetComputeGradient(false);
|
||||
extraOutput.SetComputeDivergence(true);
|
||||
extraOutput.SetComputeVorticity(true);
|
||||
extraOutput.SetComputeQCriterion(false);
|
||||
|
||||
vtkm::worklet::CellGradient gradient;
|
||||
auto result = gradient.Run(
|
||||
dataSet.GetCellSet(), dataSet.GetCoordinateSystem(), input, extraOutput, DeviceAdapter());
|
||||
|
||||
//Verify that the result is 0 size
|
||||
VTKM_TEST_ASSERT((result.GetNumberOfValues() == 0), "Gradient field shouldn't be generated");
|
||||
//Verify that the extra arrays are the correct size
|
||||
VTKM_TEST_ASSERT((extraOutput.Gradient.GetNumberOfValues() == 0),
|
||||
"Gradient field shouldn't be generated");
|
||||
VTKM_TEST_ASSERT((extraOutput.Divergence.GetNumberOfValues() == 4),
|
||||
"Divergence field should be generated");
|
||||
VTKM_TEST_ASSERT((extraOutput.Vorticity.GetNumberOfValues() == 4),
|
||||
"Vorticity field should be generated");
|
||||
VTKM_TEST_ASSERT((extraOutput.QCriterion.GetNumberOfValues() == 0),
|
||||
"QCriterion field shouldn't be generated");
|
||||
|
||||
//Verify the contents of the other arrays
|
||||
vtkm::Vec<vtkm::Vec<vtkm::Float64, 3>, 3> expected_gradients[4] = {
|
||||
{ { 10.025, 10.025, 10.025 }, { 30.075, 30.075, 30.075 }, { 60.125, 60.125, 60.125 } },
|
||||
{ { 10.025, 10.025, 10.025 }, { 30.075, 30.075, 30.075 }, { 60.125, 60.125, 60.125 } },
|
||||
{ { 10.025, 10.025, 10.025 }, { 30.075, 30.075, 30.075 }, { 60.175, 60.175, 60.175 } },
|
||||
{ { 10.025, 10.025, 10.025 }, { 30.075, 30.075, 30.075 }, { 60.175, 60.175, 60.175 } }
|
||||
};
|
||||
|
||||
for (int i = 0; i < 4; ++i)
|
||||
{
|
||||
vtkm::Vec<vtkm::Vec<vtkm::Float64, 3>, 3> eg = expected_gradients[i];
|
||||
|
||||
vtkm::Float64 d = extraOutput.Divergence.GetPortalConstControl().Get(i);
|
||||
VTKM_TEST_ASSERT(test_equal((eg[0][0] + eg[1][1] + eg[2][2]), d),
|
||||
"Wrong result for Divergence on 3D uniform data");
|
||||
|
||||
vtkm::Vec<vtkm::Float64, 3> ev(eg[2][1] - eg[1][2], eg[0][2] - eg[2][0], eg[1][0] - eg[0][1]);
|
||||
vtkm::Vec<vtkm::Float64, 3> v = extraOutput.Vorticity.GetPortalConstControl().Get(i);
|
||||
VTKM_TEST_ASSERT(test_equal(ev, v), "Wrong result for Vorticity on 3D uniform data");
|
||||
}
|
||||
}
|
||||
|
||||
@ -156,6 +249,7 @@ void TestCellGradient()
|
||||
TestCellGradientUniform2D<DeviceAdapter>();
|
||||
TestCellGradientUniform3D<DeviceAdapter>();
|
||||
TestCellGradientUniform3DWithVectorField<DeviceAdapter>();
|
||||
TestCellGradientUniform3DWithVectorField2<DeviceAdapter>();
|
||||
TestCellGradientExplicit<DeviceAdapter>();
|
||||
}
|
||||
}
|
||||
|
||||
@ -66,7 +66,10 @@ void TestPointGradientUniform3D()
|
||||
gradient.Run(dataSet.GetCellSet(), dataSet.GetCoordinateSystem(), fieldArray, DeviceAdapter());
|
||||
|
||||
vtkm::Vec<vtkm::Float32, 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 },
|
||||
{ 10.0f, 30.f, 60.1f },
|
||||
{ 10.0f, 30.1f, 60.1f },
|
||||
{ 10.0f, 30.1f, 60.2f },
|
||||
{ 10.1f, 30.f, 60.2f },
|
||||
};
|
||||
for (int i = 0; i < 4; ++i)
|
||||
{
|
||||
@ -119,6 +122,84 @@ void TestPointGradientUniform3DWithVectorField()
|
||||
}
|
||||
}
|
||||
|
||||
template <typename DeviceAdapter>
|
||||
void TestPointGradientUniform3DWithVectorField2()
|
||||
{
|
||||
std::cout << "Testing PointGradient Worklet with a vector field on 3D strucutred data"
|
||||
<< std::endl
|
||||
<< "Disabling Gradient computation and enabling Divergence, Vorticity, and QCriterion"
|
||||
<< 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);
|
||||
|
||||
vtkm::worklet::GradientOutputFields<vtkm::Vec<vtkm::Float64, 3>> extraOutput;
|
||||
extraOutput.SetComputeGradient(false);
|
||||
extraOutput.SetComputeDivergence(true);
|
||||
extraOutput.SetComputeVorticity(true);
|
||||
extraOutput.SetComputeQCriterion(true);
|
||||
|
||||
vtkm::worklet::PointGradient gradient;
|
||||
auto result = gradient.Run(
|
||||
dataSet.GetCellSet(), dataSet.GetCoordinateSystem(), input, extraOutput, DeviceAdapter());
|
||||
|
||||
//Verify that the result is 0 size
|
||||
VTKM_TEST_ASSERT((result.GetNumberOfValues() == 0), "Gradient field shouldn't be generated");
|
||||
//Verify that the extra arrays are the correct size
|
||||
VTKM_TEST_ASSERT((extraOutput.Gradient.GetNumberOfValues() == 0),
|
||||
"Gradient field shouldn't be generated");
|
||||
VTKM_TEST_ASSERT((extraOutput.Divergence.GetNumberOfValues() == nVerts),
|
||||
"Divergence field should be generated");
|
||||
VTKM_TEST_ASSERT((extraOutput.Vorticity.GetNumberOfValues() == nVerts),
|
||||
"Vorticity field should be generated");
|
||||
VTKM_TEST_ASSERT((extraOutput.QCriterion.GetNumberOfValues() == nVerts),
|
||||
"QCriterion field should be generated");
|
||||
|
||||
vtkm::Vec<vtkm::Vec<vtkm::Float64, 3>, 3> expected_gradients[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> eg = expected_gradients[i];
|
||||
|
||||
vtkm::Float64 d = extraOutput.Divergence.GetPortalConstControl().Get(i);
|
||||
|
||||
VTKM_TEST_ASSERT(test_equal((eg[0][0] + eg[1][1] + eg[2][2]), d),
|
||||
"Wrong result for Divergence on 3D uniform data");
|
||||
|
||||
vtkm::Vec<vtkm::Float64, 3> ev(eg[2][1] - eg[1][2], eg[0][2] - eg[2][0], eg[1][0] - eg[0][1]);
|
||||
vtkm::Vec<vtkm::Float64, 3> v = extraOutput.Vorticity.GetPortalConstControl().Get(i);
|
||||
VTKM_TEST_ASSERT(test_equal(ev, v), "Wrong result for Vorticity on 3D uniform data");
|
||||
|
||||
const vtkm::Vec<vtkm::Float64, 3> es(
|
||||
eg[2][1] + eg[1][2], eg[0][2] + eg[2][0], eg[1][0] + eg[0][1]);
|
||||
const vtkm::Vec<vtkm::Float64, 3> ed(eg[0][0], eg[1][1], eg[2][2]);
|
||||
|
||||
//compute QCriterion
|
||||
vtkm::Float64 qcriterion =
|
||||
((vtkm::dot(ev, ev) / 2.0f) - (vtkm::dot(ed, ed) + (vtkm::dot(es, es) / 2.0f))) / 2.0f;
|
||||
|
||||
vtkm::Float64 q = extraOutput.QCriterion.GetPortalConstControl().Get(i);
|
||||
|
||||
VTKM_TEST_ASSERT(
|
||||
test_equal(qcriterion, q),
|
||||
"Wrong result for QCriterion field of PointGradient worklet on 3D uniform data");
|
||||
}
|
||||
}
|
||||
|
||||
template <typename DeviceAdapter>
|
||||
void TestPointGradientExplicit()
|
||||
{
|
||||
@ -148,6 +229,7 @@ void TestPointGradient()
|
||||
TestPointGradientUniform2D<DeviceAdapter>();
|
||||
TestPointGradientUniform3D<DeviceAdapter>();
|
||||
TestPointGradientUniform3DWithVectorField<DeviceAdapter>();
|
||||
TestPointGradientUniform3DWithVectorField2<DeviceAdapter>();
|
||||
TestPointGradientExplicit<DeviceAdapter>();
|
||||
}
|
||||
}
|
||||
|
||||
Loading…
Reference in New Issue
Block a user