Migrate Gradient filter

This commit is contained in:
Li-Ta Lo 2022-01-30 12:18:22 -07:00
parent 4fa5f2f29c
commit 59cc8cbaeb
31 changed files with 330 additions and 944 deletions

@ -29,7 +29,6 @@
#include <vtkm/filter/CellAverage.h>
#include <vtkm/filter/FieldSelection.h>
#include <vtkm/filter/Gradient.h>
#include <vtkm/filter/PointAverage.h>
#include <vtkm/filter/PolicyBase.h>
#include <vtkm/filter/Tetrahedralize.h>
@ -42,6 +41,7 @@
#include <vtkm/filter/entity_extraction/ExternalFaces.h>
#include <vtkm/filter/entity_extraction/Threshold.h>
#include <vtkm/filter/entity_extraction/ThresholdPoints.h>
#include <vtkm/filter/vector_calculus/Gradient.h>
#include <vtkm/io/VTKDataSetReader.h>
@ -117,7 +117,7 @@ void BenchGradient(::benchmark::State& state, int options)
{
const vtkm::cont::DeviceAdapterId device = Config.Device;
vtkm::filter::Gradient filter;
vtkm::filter::vector_calculus::Gradient filter;
if (options & ScalarInput)
{

@ -25,12 +25,12 @@
#include <vtkm/cont/internal/OptionParser.h>
#include <vtkm/filter/Gradient.h>
#include <vtkm/filter/Streamline.h>
#include <vtkm/filter/Tetrahedralize.h>
#include <vtkm/filter/Tube.h>
#include <vtkm/filter/contour/Contour.h>
#include <vtkm/filter/contour/Slice.h>
#include <vtkm/filter/vector_calculus/Gradient.h>
#include <vtkm/rendering/Actor.h>
#include <vtkm/rendering/CanvasRayTracer.h>
@ -119,7 +119,7 @@ void BuildInputDataSet(uint32_t cycle, bool isStructured, bool isMultiBlock, vtk
}
// Generate Perln Noise Gradient point vector field
vtkm::filter::Gradient gradientFilter;
vtkm::filter::vector_calculus::Gradient gradientFilter;
gradientFilter.SetActiveField(PointScalarsName, vtkm::cont::Field::Association::POINTS);
gradientFilter.SetComputePointGradient(true);
gradientFilter.SetOutputFieldName(PointVectorsName);

@ -15,7 +15,7 @@
#include <vtkm/cont/openmp/DeviceAdapterOpenMP.h>
#include <vtkm/cont/tbb/DeviceAdapterTBB.h>
#include <vtkm/filter/Gradient.h>
#include <vtkm/filter/vector_calculus/Gradient.h>
namespace
@ -206,7 +206,7 @@ inline VTKM_CONT vtkm::cont::PartitionedDataSet MultiDeviceGradient::PrepareForE
//Step 3. Construct the filter we want to run on each partition
vtkm::filter::Gradient gradient;
vtkm::filter::vector_calculus::Gradient gradient;
gradient.SetComputePointGradient(this->GetComputePointGradient());
gradient.SetActiveField(this->GetActiveFieldName());
@ -218,7 +218,7 @@ inline VTKM_CONT vtkm::cont::PartitionedDataSet MultiDeviceGradient::PrepareForE
vtkm::cont::DataSet input = *partition;
this->Queue.push( //build a lambda that is the work to do
[=]() {
vtkm::filter::Gradient perThreadGrad = gradient;
vtkm::filter::vector_calculus::Gradient perThreadGrad = gradient;
vtkm::cont::DataSet result = perThreadGrad.Execute(input);
outPtr->ReplacePartition(0, result);
@ -237,7 +237,7 @@ inline VTKM_CONT vtkm::cont::PartitionedDataSet MultiDeviceGradient::PrepareForE
//blocking manner
this->Queue.push( //build a lambda that is the work to do
[=]() {
vtkm::filter::Gradient perThreadGrad = gradient;
vtkm::filter::vector_calculus::Gradient perThreadGrad = gradient;
vtkm::cont::DataSet result = perThreadGrad.Execute(input);
outPtr->ReplacePartition(index, result);

@ -8,8 +8,6 @@
## PURPOSE. See the above copyright notice for more information.
##============================================================================
vtkm_add_instantiations(GradientInstantiations FILTER Gradient)
set(deprecated_headers
CellSetConnectivity.h
CleanGrid.h
@ -25,6 +23,7 @@ set(deprecated_headers
ExtractStructured.h
GenerateIds.h
GhostCellRemove.h
Gradient.h
Histogram.h
ImageConnectivity.h
Mask.h
@ -164,19 +163,6 @@ set(extra_sources_device
particleadvection/ParticleMessenger.cxx
)
set(gradient_headers
Gradient.h
)
set(gradient_header_template_sources
Gradient.hxx
)
set(gradient_sources_device
${GradientInstantiations}
)
set(core_headers
NewFilter.h
NewFilterField.h
@ -218,19 +204,10 @@ vtkm_library(
USE_VTKM_JOB_POOL
)
vtkm_library(
NAME vtkm_filter_gradient
TEMPLATE_SOURCES ${gradient_header_template_sources}
HEADERS ${gradient_headers}
DEVICE_SOURCES ${gradient_sources_device}
USE_VTKM_JOB_POOL
)
set_target_properties(
vtkm_filter_common
vtkm_filter_core
vtkm_filter_extra
vtkm_filter_gradient
PROPERTIES
UNITY_BUILD ON
UNITY_BUILD_MODE BATCH
@ -239,15 +216,13 @@ set_target_properties(
target_link_libraries(vtkm_filter_common PUBLIC vtkm_filter_core vtkm_worklet) # TODO: deprecate vtkm_filter_common
target_link_libraries(vtkm_filter_core PUBLIC vtkm_cont vtkm_worklet)
target_link_libraries(vtkm_filter_extra PUBLIC vtkm_filter_common)
target_link_libraries(vtkm_filter_gradient PUBLIC vtkm_filter_common)
if (VTKm_ENABLE_MPI)
target_link_libraries(vtkm_filter_common PUBLIC MPI::MPI_CXX)
target_link_libraries(vtkm_filter_extra PUBLIC MPI::MPI_CXX)
target_link_libraries(vtkm_filter_gradient PUBLIC MPI::MPI_CXX)
endif()
target_link_libraries(vtkm_filter PUBLIC INTERFACE
vtkm_filter_extra
vtkm_filter_gradient
vtkm_filter_common
vtkm_filter_core
)

@ -7,239 +7,33 @@
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#ifndef vtk_m_filter_Gradient_h
#define vtk_m_filter_Gradient_h
#include <vtkm/filter/vtkm_filter_gradient_export.h>
#include <vtkm/filter/FilterField.h>
#include <vtkm/filter/Instantiations.h>
#include <vtkm/cont/ArrayHandleSOA.h>
#include <vtkm/Deprecated.h>
#include <vtkm/filter/vector_calculus/Gradient.h>
namespace vtkm
{
namespace filter
{
/// \brief A general filter for gradient estimation.
/// Estimates the gradient of a point field in a data set. The created gradient array
/// can be determined at either each point location or at the center of each cell.
///
/// The default for the filter is output as cell centered gradients.
/// To enable point based gradient computation enable \c SetComputePointGradient
///
/// Note: If no explicit name for the output field is provided the filter will
/// default to "Gradients"
class VTKM_FILTER_GRADIENT_EXPORT Gradient : public vtkm::filter::FilterField<Gradient>
VTKM_DEPRECATED(1.8,
"Use vtkm/filter/vector_calculus/Gradient.h instead of vtkm/filter/Gradient.h.")
inline void Gradient_deprecated() {}
inline void Gradient_deprecated_warning()
{
public:
using SupportedTypes = vtkm::List<vtkm::Float32, vtkm::Float64, vtkm::Vec3f_32, vtkm::Vec3f_64>;
Gradient_deprecated();
}
/// When this flag is on (default is off), the gradient filter will provide a
/// point based gradients, which are significantly more costly since for each
/// point we need to compute the gradient of each cell that uses it.
void SetComputePointGradient(bool enable) { ComputePointGradient = enable; }
bool GetComputePointGradient() const { return ComputePointGradient; }
/// Add divergence field to the output data. The name of the array
/// will be Divergence and will be a cell field unless \c ComputePointGradient
/// is enabled. 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 name of the array
/// will be Vorticity and will be a cell field unless \c ComputePointGradient
/// is enabled. 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 name of the array
/// will be QCriterion and will be a cell field unless \c ComputePointGradient
/// is enabled. 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 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; }
/// Make the vector gradient output format be in FORTRAN Column-major order.
/// This is only used when the input field is a vector field ( 3 components ).
/// Enabling column-major is important if integrating with other projects
/// such as VTK.
/// Default: Row Order
void SetColumnMajorOrdering() { RowOrdering = false; }
/// Make the vector gradient output format be in C Row-major order.
/// This is only used when the input field is a vector field ( 3 components ).
/// Default: Row Order
void SetRowMajorOrdering() { RowOrdering = true; }
void SetDivergenceName(const std::string& name) { this->DivergenceName = name; }
const std::string& GetDivergenceName() const { return this->DivergenceName; }
void SetVorticityName(const std::string& name) { this->VorticityName = name; }
const std::string& GetVorticityName() const { return this->VorticityName; }
void SetQCriterionName(const std::string& name) { this->QCriterionName = name; }
const std::string& GetQCriterionName() const { return this->QCriterionName; }
template <typename T, typename StorageType, typename DerivedPolicy>
vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet& input,
const vtkm::cont::ArrayHandle<T, StorageType>& field,
const vtkm::filter::FieldMetadata& fieldMeta,
vtkm::filter::PolicyBase<DerivedPolicy> policy);
VTKM_CONT
Filter* Clone() const override
{
Gradient* clone = new Gradient();
clone->CopyStateFrom(this);
return clone;
}
VTKM_CONT
bool CanThread() const override { return true; }
protected:
VTKM_CONT
void CopyStateFrom(const Gradient* gradient)
{
this->FilterField<Gradient>::CopyStateFrom(gradient);
this->ComputePointGradient = gradient->ComputePointGradient;
this->ComputeDivergence = gradient->ComputeDivergence;
this->ComputeVorticity = gradient->ComputeVorticity;
this->ComputeQCriterion = gradient->ComputeQCriterion;
this->StoreGradient = gradient->StoreGradient;
this->RowOrdering = gradient->RowOrdering;
}
private:
bool ComputePointGradient = false;
bool ComputeDivergence = false;
bool ComputeVorticity = false;
bool ComputeQCriterion = false;
bool StoreGradient = true;
bool RowOrdering = true;
std::string DivergenceName = "Divergence";
std::string GradientsName = "Gradients";
std::string QCriterionName = "QCriterion";
std::string VorticityName = "Vorticity";
class VTKM_DEPRECATED(1.8, "Use vtkm::filter::vector_calculus::Gradient.") Gradient
: public vtkm::filter::vector_calculus::Gradient
{
using vector_calculus::Gradient::Gradient;
};
#ifndef vtkm_filter_Gradient_cxx
VTKM_INSTANTIATION_BEGIN
extern template VTKM_FILTER_GRADIENT_TEMPLATE_EXPORT vtkm::cont::DataSet Gradient::DoExecute(
const vtkm::cont::DataSet&,
const vtkm::cont::ArrayHandle<vtkm::Float32>&,
const vtkm::filter::FieldMetadata&,
vtkm::filter::PolicyBase<vtkm::filter::PolicyDefault>);
VTKM_INSTANTIATION_END
VTKM_INSTANTIATION_BEGIN
extern template VTKM_FILTER_GRADIENT_TEMPLATE_EXPORT vtkm::cont::DataSet Gradient::DoExecute(
const vtkm::cont::DataSet&,
const vtkm::cont::ArrayHandle<vtkm::Float64>&,
const vtkm::filter::FieldMetadata&,
vtkm::filter::PolicyBase<vtkm::filter::PolicyDefault>);
VTKM_INSTANTIATION_END
VTKM_INSTANTIATION_BEGIN
extern template VTKM_FILTER_GRADIENT_TEMPLATE_EXPORT vtkm::cont::DataSet Gradient::DoExecute(
const vtkm::cont::DataSet&,
const vtkm::cont::ArrayHandle<vtkm::Vec3f_32>&,
const vtkm::filter::FieldMetadata&,
vtkm::filter::PolicyBase<vtkm::filter::PolicyDefault>);
VTKM_INSTANTIATION_END
VTKM_INSTANTIATION_BEGIN
extern template VTKM_FILTER_GRADIENT_TEMPLATE_EXPORT vtkm::cont::DataSet Gradient::DoExecute(
const vtkm::cont::DataSet&,
const vtkm::cont::ArrayHandle<vtkm::Vec3f_32, vtkm::cont::StorageTagSOA>&,
const vtkm::filter::FieldMetadata&,
vtkm::filter::PolicyBase<vtkm::filter::PolicyDefault>);
VTKM_INSTANTIATION_END
VTKM_INSTANTIATION_BEGIN
extern template VTKM_FILTER_GRADIENT_TEMPLATE_EXPORT vtkm::cont::DataSet Gradient::DoExecute(
const vtkm::cont::DataSet&,
const vtkm::cont::ArrayHandle<vtkm::Vec3f_64>&,
const vtkm::filter::FieldMetadata&,
vtkm::filter::PolicyBase<vtkm::filter::PolicyDefault>);
VTKM_INSTANTIATION_END
VTKM_INSTANTIATION_BEGIN
extern template VTKM_FILTER_GRADIENT_TEMPLATE_EXPORT vtkm::cont::DataSet Gradient::DoExecute(
const vtkm::cont::DataSet&,
const vtkm::cont::ArrayHandle<vtkm::Vec3f_64, vtkm::cont::StorageTagSOA>&,
const vtkm::filter::FieldMetadata&,
vtkm::filter::PolicyBase<vtkm::filter::PolicyDefault>);
VTKM_INSTANTIATION_END
VTKM_INSTANTIATION_BEGIN
extern template VTKM_FILTER_GRADIENT_TEMPLATE_EXPORT vtkm::cont::DataSet Gradient::DoExecute(
const vtkm::cont::DataSet&,
const vtkm::cont::ArrayHandle<vtkm::Vec3f, vtkm::cont::StorageTagUniformPoints>&,
const vtkm::filter::FieldMetadata&,
vtkm::filter::PolicyBase<vtkm::filter::PolicyDefault>);
VTKM_INSTANTIATION_END
VTKM_INSTANTIATION_BEGIN
extern template VTKM_FILTER_GRADIENT_TEMPLATE_EXPORT vtkm::cont::DataSet Gradient::DoExecute(
const vtkm::cont::DataSet&,
const vtkm::cont::ArrayHandle<
vtkm::Vec3f_32,
vtkm::cont::StorageTagCartesianProduct<vtkm::cont::StorageTagBasic,
vtkm::cont::StorageTagBasic,
vtkm::cont::StorageTagBasic>>&,
const vtkm::filter::FieldMetadata&,
vtkm::filter::PolicyBase<vtkm::filter::PolicyDefault>);
VTKM_INSTANTIATION_END
VTKM_INSTANTIATION_BEGIN
extern template VTKM_FILTER_GRADIENT_TEMPLATE_EXPORT vtkm::cont::DataSet Gradient::DoExecute(
const vtkm::cont::DataSet&,
const vtkm::cont::ArrayHandle<
vtkm::Vec3f_64,
vtkm::cont::StorageTagCartesianProduct<vtkm::cont::StorageTagBasic,
vtkm::cont::StorageTagBasic,
vtkm::cont::StorageTagBasic>>&,
const vtkm::filter::FieldMetadata&,
vtkm::filter::PolicyBase<vtkm::filter::PolicyDefault>);
VTKM_INSTANTIATION_END
VTKM_INSTANTIATION_BEGIN
#ifdef VTKM_ADD_XGC_DEFAULT_TYPES
extern template VTKM_FILTER_GRADIENT_TEMPLATE_EXPORT vtkm::cont::DataSet Gradient::DoExecute(
const vtkm::cont::DataSet&,
const vtkm::cont::ArrayHandle<vtkm::Vec<float, 3>, vtkm::cont::StorageTagXGCCoordinates>&,
const vtkm::filter::FieldMetadata&,
vtkm::filter::PolicyBase<vtkm::filter::PolicyDefault>);
#endif
VTKM_INSTANTIATION_END
VTKM_INSTANTIATION_BEGIN
#ifdef VTKM_ADD_XGC_DEFAULT_TYPES
extern template VTKM_FILTER_GRADIENT_TEMPLATE_EXPORT vtkm::cont::DataSet Gradient::DoExecute(
const vtkm::cont::DataSet&,
const vtkm::cont::ArrayHandle<vtkm::Vec<double, 3>, vtkm::cont::StorageTagXGCCoordinates>&,
const vtkm::filter::FieldMetadata&,
vtkm::filter::PolicyBase<vtkm::filter::PolicyDefault>);
#endif
VTKM_INSTANTIATION_END
#endif //vtkm_filter_Gradient_cxx
}
} // namespace vtkm::filter
#endif // vtk_m_filter_Gradient_h
#endif //vtk_m_filter_Gradient_h

@ -1,129 +0,0 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
//
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#ifndef vtk_m_filter_Gradient_hxx
#define vtk_m_filter_Gradient_hxx
#include <vtkm/cont/ErrorFilterExecution.h>
#include <vtkm/cont/UnknownCellSet.h>
#include <vtkm/worklet/Gradient.h>
namespace
{
//-----------------------------------------------------------------------------
template <typename T, typename S>
inline void transpose_3x3(vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Vec<T, 3>, 3>, S>& field)
{
vtkm::worklet::gradient::Transpose3x3<T> transpose;
transpose.Run(field);
}
//-----------------------------------------------------------------------------
template <typename T, typename S>
inline void transpose_3x3(vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>, S>&)
{ //This is not a 3x3 matrix so no transpose needed
}
} //namespace
namespace vtkm
{
namespace filter
{
//-----------------------------------------------------------------------------
template <typename T, typename StorageType, typename DerivedPolicy>
inline vtkm::cont::DataSet Gradient::DoExecute(
const vtkm::cont::DataSet& input,
const vtkm::cont::ArrayHandle<T, StorageType>& inField,
const vtkm::filter::FieldMetadata& fieldMetadata,
vtkm::filter::PolicyBase<DerivedPolicy> policy)
{
if (!fieldMetadata.IsPointField())
{
throw vtkm::cont::ErrorFilterExecution("Point field expected.");
}
constexpr bool isVector = std::is_same<typename vtkm::VecTraits<T>::HasMultipleComponents,
vtkm::VecTraitsTagMultipleComponents>::value;
if (GetComputeQCriterion() && !isVector)
{
throw vtkm::cont::ErrorFilterExecution("scalar gradients can't generate qcriterion");
}
if (GetComputeVorticity() && !isVector)
{
throw vtkm::cont::ErrorFilterExecution("scalar gradients can't generate vorticity");
}
const vtkm::cont::UnknownCellSet& cells = input.GetCellSet();
const vtkm::cont::CoordinateSystem& coords =
input.GetCoordinateSystem(this->GetActiveCoordinateSystemIndex());
std::string outputName = this->GetOutputFieldName();
if (outputName.empty())
{
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)
{
vtkm::worklet::PointGradient gradient;
outArray = gradient.Run(
vtkm::filter::ApplyPolicyCellSet(cells, policy, *this), coords, inField, gradientfields);
}
else
{
vtkm::worklet::CellGradient gradient;
outArray = gradient.Run(
vtkm::filter::ApplyPolicyCellSet(cells, policy, *this), coords, inField, gradientfields);
}
if (!this->RowOrdering)
{
transpose_3x3(outArray);
}
vtkm::cont::Field::Association fieldAssociation(this->ComputePointGradient
? vtkm::cont::Field::Association::POINTS
: vtkm::cont::Field::Association::CELL_SET);
vtkm::cont::DataSet result;
result.CopyStructure(input);
result.AddField(vtkm::cont::Field{ outputName, fieldAssociation, outArray });
if (this->GetComputeDivergence() && isVector)
{
result.AddField(
vtkm::cont::Field{ this->GetDivergenceName(), fieldAssociation, gradientfields.Divergence });
}
if (this->GetComputeVorticity() && isVector)
{
result.AddField(
vtkm::cont::Field{ this->GetVorticityName(), fieldAssociation, gradientfields.Vorticity });
}
if (this->GetComputeQCriterion() && isVector)
{
result.AddField(
vtkm::cont::Field{ this->GetQCriterionName(), fieldAssociation, gradientfields.QCriterion });
}
return result;
}
}
} // namespace vtkm::filter
#endif

@ -18,7 +18,7 @@
#include <vtkm/filter/contour/worklet/FlyingEdgesTables.h>
#include <vtkm/VectorAnalysis.h>
#include <vtkm/worklet/gradient/StructuredPointGradient.h>
#include <vtkm/filter/vector_calculus/worklet/gradient/StructuredPointGradient.h>
namespace vtkm
{

@ -30,9 +30,9 @@
#include <vtkm/filter/contour/worklet/CommonState.h>
#include <vtkm/filter/contour/worklet/MarchingCellTables.h>
#include <vtkm/filter/vector_calculus/worklet/gradient/PointGradient.h>
#include <vtkm/filter/vector_calculus/worklet/gradient/StructuredPointGradient.h>
#include <vtkm/worklet/WorkletReduceByKey.h>
#include <vtkm/worklet/gradient/PointGradient.h>
#include <vtkm/worklet/gradient/StructuredPointGradient.h>
namespace vtkm
{

@ -22,8 +22,6 @@ set(unit_tests
UnitTestFieldMetadata.cxx
UnitTestFieldSelection.cxx
UnitTestFieldToColors.cxx
UnitTestGradientExplicit.cxx
UnitTestGradientUniform.cxx
UnitTestGhostCellClassify.cxx
UnitTestImageDifferenceFilter.cxx
UnitTestImageMedianFilter.cxx

@ -13,10 +13,10 @@
#include <vtkm/cont/testing/MakeTestDataSet.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/filter/Gradient.h>
#include <vtkm/filter/clean_grid/CleanGrid.h>
#include <vtkm/filter/contour/ClipWithField.h>
#include <vtkm/filter/contour/Contour.h>
#include <vtkm/filter/vector_calculus/Gradient.h>
#include <vtkm/io/VTKDataSetReader.h>
#include <vtkm/source/Tangle.h>
@ -147,7 +147,7 @@ void TestMultiBlockFilter()
results.clear();
for (const auto doThreading : flags)
{
vtkm::filter::Gradient grad;
vtkm::filter::vector_calculus::Gradient grad;
grad.SetRunMultiThreadedFilter(doThreading);
grad.SetComputePointGradient(true);
grad.SetActiveField("tangle");

@ -10,10 +10,12 @@
set(vector_calculus_headers
CrossProduct.h
DotProduct.h
Gradient.h
)
set(vector_calculus_sources_device
CrossProduct.cxx
DotProduct.cxx
Gradient.cxx
)
vtkm_library(
@ -26,6 +28,12 @@ vtkm_library(
target_link_libraries(vtkm_filter_vector_calculus PUBLIC vtkm_worklet vtkm_filter_core)
target_link_libraries(vtkm_filter PUBLIC INTERFACE vtkm_filter_vector_calculus)
if (VTKm_ENABLE_MPI)
# TODO: is this necessary?
# target_link_libraries(vtkm_filter_gradient PUBLIC MPI::MPI_CXX)
endif()
add_subdirectory(worklet)
#-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -
if (VTKm_ENABLE_TESTING)
add_subdirectory(testing)

@ -0,0 +1,139 @@
//============================================================================
// 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.
//============================================================================
#include <vtkm/cont/ErrorFilterExecution.h>
#include <vtkm/cont/UnknownCellSet.h>
#include <vtkm/filter/vector_calculus/Gradient.h>
#include <vtkm/filter/vector_calculus/worklet/Gradient.h>
namespace
{
//-----------------------------------------------------------------------------
template <typename T, typename S>
inline void transpose_3x3(vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Vec<T, 3>, 3>, S>& field)
{
vtkm::worklet::gradient::Transpose3x3<T> transpose;
transpose.Run(field);
}
//-----------------------------------------------------------------------------
template <typename T, typename S>
inline void transpose_3x3(vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>, S>&)
{ //This is not a 3x3 matrix so no transpose needed
}
} //namespace
namespace vtkm
{
namespace filter
{
namespace vector_calculus
{
//-----------------------------------------------------------------------------
vtkm::cont::DataSet Gradient::DoExecute(const vtkm::cont::DataSet& inputDataSet)
{
auto field = this->GetFieldFromDataSet(inputDataSet);
if (!field.IsFieldPoint())
{
throw vtkm::cont::ErrorFilterExecution("Point field expected.");
}
const bool isVector = field.GetData().GetNumberOfComponents() == 3;
if (GetComputeQCriterion() && !isVector)
{
throw vtkm::cont::ErrorFilterExecution("scalar gradients can't generate qcriterion");
}
if (GetComputeVorticity() && !isVector)
{
throw vtkm::cont::ErrorFilterExecution("scalar gradients can't generate vorticity");
}
const vtkm::cont::UnknownCellSet& inputCellSet = inputDataSet.GetCellSet();
const vtkm::cont::CoordinateSystem& coords =
inputDataSet.GetCoordinateSystem(this->GetActiveCoordinateSystemIndex());
vtkm::cont::UnknownArrayHandle gradientArray;
vtkm::cont::UnknownArrayHandle divergenceArray;
vtkm::cont::UnknownArrayHandle vorticityArray;
vtkm::cont::UnknownArrayHandle qcriterionArray;
// TODO: there are a humungous number of (weak) symbols in the .o file. Investigate if
// they are all legit.
auto resolveType = [&, this](const auto& concrete) {
// use std::decay to remove const ref from the decltype of concrete.
using T = typename std::decay_t<decltype(concrete)>::ValueType;
vtkm::worklet::GradientOutputFields<T> gradientfields(this->GetComputeGradient(),
this->GetComputeDivergence(),
this->GetComputeVorticity(),
this->GetComputeQCriterion());
vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>> result;
if (this->ComputePointGradient)
{
vtkm::worklet::PointGradient gradient;
result = gradient.Run(inputCellSet, coords, concrete, gradientfields);
}
else
{
vtkm::worklet::CellGradient gradient;
result = gradient.Run(inputCellSet, coords, concrete, gradientfields);
}
if (!this->RowOrdering)
{
transpose_3x3(result);
}
gradientArray = result;
divergenceArray = gradientfields.Divergence;
vorticityArray = gradientfields.Vorticity;
qcriterionArray = gradientfields.QCriterion;
};
// TODO: do we need to deal with vtkm::cont::StorageTagXGCCoordinates? or it is already part of
// VTKM_DEFAULT_STORAGE_LIST when enabled.
field.GetData().CastAndCallForTypesWithFloatFallback<SupportedTypes, VTKM_DEFAULT_STORAGE_LIST>(
resolveType);
// This copies the CellSet and Fields to be passed from inputDataSet to outputDataSet
vtkm::cont::DataSet outputDataSet = this->CreateResult(inputDataSet);
std::string outputName = this->GetOutputFieldName();
if (outputName.empty())
{
outputName = this->GradientsName;
}
vtkm::cont::Field::Association fieldAssociation(this->ComputePointGradient
? vtkm::cont::Field::Association::POINTS
: vtkm::cont::Field::Association::CELL_SET);
outputDataSet.AddField(vtkm::cont::Field{ outputName, fieldAssociation, gradientArray });
if (this->GetComputeDivergence() && isVector)
{
outputDataSet.AddField(
vtkm::cont::Field{ this->GetDivergenceName(), fieldAssociation, divergenceArray });
}
if (this->GetComputeVorticity() && isVector)
{
outputDataSet.AddField(
vtkm::cont::Field{ this->GetVorticityName(), fieldAssociation, vorticityArray });
}
if (this->GetComputeQCriterion() && isVector)
{
outputDataSet.AddField(
vtkm::cont::Field{ this->GetQCriterionName(), fieldAssociation, qcriterionArray });
}
return outputDataSet;
}
}
}
} // namespace vtkm::filter

@ -0,0 +1,112 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
//
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#ifndef vtk_m_filter_Gradient_h
#define vtk_m_filter_Gradient_h
#include <vtkm/filter/NewFilterField.h>
#include <vtkm/filter/vector_calculus/vtkm_filter_vector_calculus_export.h>
namespace vtkm
{
namespace filter
{
namespace vector_calculus
{
/// \brief A general filter for gradient estimation.
/// Estimates the gradient of a point field in a data set. The created gradient array
/// can be determined at either each point location or at the center of each cell.
///
/// The default for the filter is output as cell centered gradients.
/// To enable point based gradient computation enable \c SetComputePointGradient
///
/// Note: If no explicit name for the output field is provided the filter will
/// default to "Gradients"
class VTKM_FILTER_VECTOR_CALCULUS_EXPORT Gradient : public vtkm::filter::NewFilterField
{
public:
using SupportedTypes = vtkm::List<vtkm::Float32, vtkm::Float64, vtkm::Vec3f_32, vtkm::Vec3f_64>;
/// When this flag is on (default is off), the gradient filter will provide a
/// point based gradients, which are significantly more costly since for each
/// point we need to compute the gradient of each cell that uses it.
void SetComputePointGradient(bool enable) { ComputePointGradient = enable; }
bool GetComputePointGradient() const { return ComputePointGradient; }
/// Add divergence field to the output data. The name of the array
/// will be Divergence and will be a cell field unless \c ComputePointGradient
/// is enabled. 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 name of the array
/// will be Vorticity and will be a cell field unless \c ComputePointGradient
/// is enabled. 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 name of the array
/// will be QCriterion and will be a cell field unless \c ComputePointGradient
/// is enabled. 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 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; }
/// Make the vector gradient output format be in FORTRAN Column-major order.
/// This is only used when the input field is a vector field ( 3 components ).
/// Enabling column-major is important if integrating with other projects
/// such as VTK.
/// Default: Row Order
void SetColumnMajorOrdering() { RowOrdering = false; }
/// Make the vector gradient output format be in C Row-major order.
/// This is only used when the input field is a vector field ( 3 components ).
/// Default: Row Order
void SetRowMajorOrdering() { RowOrdering = true; }
void SetDivergenceName(const std::string& name) { this->DivergenceName = name; }
const std::string& GetDivergenceName() const { return this->DivergenceName; }
void SetVorticityName(const std::string& name) { this->VorticityName = name; }
const std::string& GetVorticityName() const { return this->VorticityName; }
void SetQCriterionName(const std::string& name) { this->QCriterionName = name; }
const std::string& GetQCriterionName() const { return this->QCriterionName; }
private:
vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet& inputDataSet) override;
bool ComputePointGradient = false;
bool ComputeDivergence = false;
bool ComputeVorticity = false;
bool ComputeQCriterion = false;
bool StoreGradient = true;
bool RowOrdering = true;
std::string DivergenceName = "Divergence";
std::string GradientsName = "Gradients";
std::string QCriterionName = "QCriterion";
std::string VorticityName = "Vorticity";
};
} // namespace vector_calculus
} // namespace filter
} // namespace vtkm::filter
#endif // vtk_m_filter_Gradient_h

@ -11,6 +11,8 @@
set(unit_tests
UnitTestCrossProductFilter.cxx
UnitTestDotProductFilter.cxx
UnitTestGradientExplicit.cxx
UnitTestGradientUniform.cxx
)
set(libraries

@ -8,7 +8,7 @@
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#include <vtkm/filter/Gradient.h>
#include <vtkm/filter/vector_calculus/Gradient.h>
#include <vtkm/cont/testing/MakeTestDataSet.h>
#include <vtkm/cont/testing/Testing.h>
@ -23,7 +23,7 @@ void TestCellGradientExplicit()
vtkm::cont::testing::MakeTestDataSet testDataSet;
vtkm::cont::DataSet dataSet = testDataSet.Make3DExplicitDataSet0();
vtkm::filter::Gradient gradient;
vtkm::filter::vector_calculus::Gradient gradient;
gradient.SetOutputFieldName("gradient");
gradient.SetActiveField("pointvar");
@ -48,7 +48,7 @@ void TestPointGradientExplicit()
vtkm::cont::testing::MakeTestDataSet testDataSet;
vtkm::cont::DataSet dataSet = testDataSet.Make3DExplicitDataSet0();
vtkm::filter::Gradient gradient;
vtkm::filter::vector_calculus::Gradient gradient;
gradient.SetComputePointGradient(true);
gradient.SetOutputFieldName("gradient");
gradient.SetActiveField("pointvar");

@ -8,8 +8,9 @@
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#include <vtkm/filter/Gradient.h>
#include <vtkm/filter/vector_calculus/Gradient.h>
#include <vtkm/cont/ErrorFilterExecution.h>
#include <vtkm/cont/testing/MakeTestDataSet.h>
#include <vtkm/cont/testing/Testing.h>
@ -23,7 +24,7 @@ void TestCellGradientUniform3D()
vtkm::cont::testing::MakeTestDataSet testDataSet;
vtkm::cont::DataSet dataSet = testDataSet.Make3DUniformDataSet0();
vtkm::filter::Gradient gradient;
vtkm::filter::vector_calculus::Gradient gradient;
gradient.SetOutputFieldName("Gradient");
gradient.SetComputeVorticity(true); //this won't work as we have a scalar field
@ -65,7 +66,7 @@ void TestCellGradientUniform3DWithVectorField()
dataSet.AddPointField("vec_pointvar", input);
//we need to add Vec3 array to the dataset
vtkm::filter::Gradient gradient;
vtkm::filter::vector_calculus::Gradient gradient;
gradient.SetOutputFieldName("vec_gradient");
gradient.SetComputeVorticity(true);
gradient.SetComputeQCriterion(true);
@ -124,7 +125,7 @@ void TestPointGradientUniform3DWithVectorField()
dataSet.AddPointField("vec_pointvar", input);
//we need to add Vec3 array to the dataset
vtkm::filter::Gradient gradient;
vtkm::filter::vector_calculus::Gradient gradient;
gradient.SetComputePointGradient(true);
gradient.SetOutputFieldName("vec_gradient");
gradient.SetActiveField("vec_pointvar");

@ -0,0 +1,18 @@
##============================================================================
## 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.
##============================================================================
set(headers
Gradient.h
)
add_subdirectory(gradient)
#-----------------------------------------------------------------------------
vtkm_declare_headers(${headers})

@ -14,14 +14,14 @@
#include <vtkm/worklet/DispatcherMapTopology.h>
#include <vtkm/worklet/DispatcherPointNeighborhood.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/StructuredPointGradient.h>
#include <vtkm/worklet/gradient/Transpose.h>
#include <vtkm/worklet/gradient/Vorticity.h>
#include <vtkm/filter/vector_calculus/worklet/gradient/CellGradient.h>
#include <vtkm/filter/vector_calculus/worklet/gradient/Divergence.h>
#include <vtkm/filter/vector_calculus/worklet/gradient/GradientOutput.h>
#include <vtkm/filter/vector_calculus/worklet/gradient/PointGradient.h>
#include <vtkm/filter/vector_calculus/worklet/gradient/QCriterion.h>
#include <vtkm/filter/vector_calculus/worklet/gradient/StructuredPointGradient.h>
#include <vtkm/filter/vector_calculus/worklet/gradient/Transpose.h>
#include <vtkm/filter/vector_calculus/worklet/gradient/Vorticity.h>
namespace vtkm
{

@ -15,7 +15,7 @@
#include <vtkm/exec/ParametricCoordinates.h>
#include <vtkm/worklet/WorkletMapTopology.h>
#include <vtkm/worklet/gradient/GradientOutput.h>
#include <vtkm/filter/vector_calculus/worklet/gradient/GradientOutput.h>
namespace vtkm
{

@ -19,9 +19,9 @@
#include <vtkm/cont/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>
#include <vtkm/filter/vector_calculus/worklet/gradient/Divergence.h>
#include <vtkm/filter/vector_calculus/worklet/gradient/QCriterion.h>
#include <vtkm/filter/vector_calculus/worklet/gradient/Vorticity.h>
namespace vtkm
{

@ -16,7 +16,7 @@
#include <vtkm/worklet/WorkletMapTopology.h>
#include <utility>
#include <vtkm/worklet/gradient/GradientOutput.h>
#include <vtkm/filter/vector_calculus/worklet/gradient/GradientOutput.h>
namespace vtkm

@ -11,8 +11,8 @@
#ifndef vtk_m_worklet_gradient_StructuredPointGradient_h
#define vtk_m_worklet_gradient_StructuredPointGradient_h
#include <vtkm/filter/vector_calculus/worklet/gradient/GradientOutput.h>
#include <vtkm/worklet/WorkletPointNeighborhood.h>
#include <vtkm/worklet/gradient/GradientOutput.h>
namespace vtkm

@ -25,7 +25,6 @@ set(headers
DispatcherPointNeighborhood.h
DispatcherReduceByKey.h
FieldStatistics.h
Gradient.h
ImageDifference.h
KdTree3D.h # Deprecated
KernelSplatter.h
@ -111,7 +110,6 @@ add_subdirectory(contourtree)
add_subdirectory(contourtree_augmented)
add_subdirectory(contourtree_distributed)
add_subdirectory(cosmotools)
add_subdirectory(gradient)
add_subdirectory(histogram)
add_subdirectory(lcs)
add_subdirectory(mir)

@ -18,7 +18,6 @@ set(unit_tests
UnitTestBoundingIntervalHierarchy.cxx
UnitTestCellAverage.cxx
UnitTestCellDeepCopy.cxx
UnitTestCellGradient.cxx
UnitTestCellMeasure.cxx
UnitTestContourTreeUniform.cxx
UnitTestContourTreeUniformAugmented.cxx
@ -36,7 +35,6 @@ set(unit_tests
UnitTestOrientNormals.cxx
UnitTestParticleAdvection.cxx
UnitTestPointElevation.cxx
UnitTestPointGradient.cxx
UnitTestPointTransform.cxx
UnitTestProbe.cxx
UnitTestScalarsToColors.cxx

@ -1,248 +0,0 @@
//============================================================================
// 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.
//============================================================================
#include <vtkm/worklet/Gradient.h>
#include <vtkm/cont/testing/MakeTestDataSet.h>
#include <vtkm/cont/testing/Testing.h>
namespace
{
void TestCellGradientUniform2D()
{
std::cout << "Testing CellGradient Worklet on 2D structured data" << std::endl;
vtkm::cont::testing::MakeTestDataSet testDataSet;
vtkm::cont::DataSet dataSet = testDataSet.Make2DUniformDataSet0();
vtkm::cont::ArrayHandle<vtkm::Float32> input;
vtkm::cont::ArrayHandle<vtkm::Vec3f_32> result;
dataSet.GetField("pointvar").GetData().AsArrayHandle(input);
vtkm::worklet::CellGradient gradient;
result = gradient.Run(dataSet.GetCellSet(), dataSet.GetCoordinateSystem(), input);
vtkm::Vec3f_32 expected[2] = { { 10, 30, 0 }, { 10, 30, 0 } };
for (int i = 0; i < 2; ++i)
{
VTKM_TEST_ASSERT(test_equal(result.ReadPortal().Get(i), expected[i]),
"Wrong result for CellGradient worklet on 2D uniform data");
}
}
void TestCellGradientUniform3D()
{
std::cout << "Testing CellGradient Worklet on 3D structured data" << std::endl;
vtkm::cont::testing::MakeTestDataSet testDataSet;
vtkm::cont::DataSet dataSet = testDataSet.Make3DUniformDataSet0();
vtkm::cont::ArrayHandle<vtkm::Float32> input;
vtkm::cont::ArrayHandle<vtkm::Vec3f_32> result;
dataSet.GetField("pointvar").GetData().AsArrayHandle(input);
vtkm::worklet::CellGradient gradient;
result = gradient.Run(dataSet.GetCellSet(), dataSet.GetCoordinateSystem(), input);
vtkm::Vec3f_32 expected[4] = {
{ 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)
{
VTKM_TEST_ASSERT(test_equal(result.ReadPortal().Get(i), expected[i]),
"Wrong result for CellGradient worklet on 3D uniform data");
}
}
void TestCellGradientUniform3DWithVectorField()
{
std::cout
<< "Testing CellGradient and QCriterion Worklet with a vector field on 3D structured 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::Vec3f_64> 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::Vec3f_64> input =
vtkm::cont::make_ArrayHandle(vec, vtkm::CopyFlag::Off);
//we need to add Vec3 array to the dataset
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Vec3f_64, 3>> result;
vtkm::worklet::GradientOutputFields<vtkm::Vec3f_64> extraOutput;
extraOutput.SetComputeDivergence(false);
extraOutput.SetComputeVorticity(false);
extraOutput.SetComputeQCriterion(true);
vtkm::worklet::CellGradient gradient;
result = gradient.Run(dataSet.GetCellSet(), dataSet.GetCoordinateSystem(), input, extraOutput);
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::Vec3f_64, 3> expected[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::Vec3f_64, 3> e = expected[i];
vtkm::Vec<vtkm::Vec3f_64, 3> r = result.ReadPortal().Get(i);
VTKM_TEST_ASSERT(test_equal(e[0], r[0]),
"Wrong result for vec field CellGradient worklet on 3D uniform data");
VTKM_TEST_ASSERT(test_equal(e[1], r[1]),
"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::Vec3f_64 v(e[1][2] - e[2][1], e[2][0] - e[0][2], e[0][1] - e[1][0]);
const vtkm::Vec3f_64 s(e[1][2] + e[2][1], e[2][0] + e[0][2], e[0][1] + e[1][0]);
const vtkm::Vec3f_64 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.ReadPortal().Get(i);
std::cout << "qcriterion expected: " << qcriterion << std::endl;
std::cout << "qcriterion actual: " << q << std::endl;
VTKM_TEST_ASSERT(
test_equal(qcriterion, q),
"Wrong result for QCriterion field of CellGradient worklet on 3D uniform data");
}
}
void TestCellGradientUniform3DWithVectorField2()
{
std::cout << "Testing CellGradient Worklet with a vector field on 3D structured 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::Vec3f_64> 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::Vec3f_64> input =
vtkm::cont::make_ArrayHandle(vec, vtkm::CopyFlag::Off);
vtkm::worklet::GradientOutputFields<vtkm::Vec3f_64> 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);
//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::Vec3f_64, 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 } }
};
auto vorticityPortal = extraOutput.Vorticity.ReadPortal();
auto divergencePortal = extraOutput.Divergence.ReadPortal();
for (int i = 0; i < 4; ++i)
{
vtkm::Vec<vtkm::Vec3f_64, 3> eg = expected_gradients[i];
vtkm::Float64 d = divergencePortal.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::Vec3f_64 ev(eg[1][2] - eg[2][1], eg[2][0] - eg[0][2], eg[0][1] - eg[1][0]);
vtkm::Vec3f_64 v = vorticityPortal.Get(i);
VTKM_TEST_ASSERT(test_equal(ev, v), "Wrong result for Vorticity on 3D uniform data");
}
}
void TestCellGradientExplicit()
{
std::cout << "Testing CellGradient Worklet on Explicit data" << std::endl;
vtkm::cont::testing::MakeTestDataSet testDataSet;
vtkm::cont::DataSet dataSet = testDataSet.Make3DExplicitDataSet0();
vtkm::cont::ArrayHandle<vtkm::Float32> input;
vtkm::cont::ArrayHandle<vtkm::Vec3f_32> result;
dataSet.GetField("pointvar").GetData().AsArrayHandle(input);
vtkm::worklet::CellGradient gradient;
result = gradient.Run(dataSet.GetCellSet(), dataSet.GetCoordinateSystem(), input);
vtkm::Vec3f_32 expected[2] = { { 10.f, 10.1f, 0.0f }, { 10.f, 10.1f, -0.0f } };
auto resultPortal = result.ReadPortal();
for (int i = 0; i < 2; ++i)
{
VTKM_TEST_ASSERT(test_equal(resultPortal.Get(i), expected[i]),
"Wrong result for CellGradient worklet on 3D explicit data");
}
}
void TestCellGradient()
{
TestCellGradientUniform2D();
TestCellGradientUniform3D();
TestCellGradientUniform3DWithVectorField();
TestCellGradientUniform3DWithVectorField2();
TestCellGradientExplicit();
}
}
int UnitTestCellGradient(int argc, char* argv[])
{
return vtkm::cont::testing::Testing::Run(TestCellGradient, argc, argv);
}

@ -1,280 +0,0 @@
//============================================================================
// 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.
//============================================================================
#include <vtkm/worklet/DispatcherMapTopology.h>
#include <vtkm/worklet/Gradient.h>
#include <vtkm/cont/testing/MakeTestDataSet.h>
#include <vtkm/cont/testing/Testing.h>
namespace
{
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::Float32> fieldArray;
dataSet.GetField("pointvar").GetData().AsArrayHandle(fieldArray);
vtkm::worklet::PointGradient gradient;
auto result = gradient.Run(dataSet.GetCellSet(), dataSet.GetCoordinateSystem(), fieldArray);
vtkm::Vec3f_32 expected[2] = { { 10, 30, 0 }, { 10, 30, 0 } };
for (int i = 0; i < 2; ++i)
{
VTKM_TEST_ASSERT(test_equal(result.ReadPortal().Get(i), expected[i]),
"Wrong result for PointGradient worklet on 2D uniform data",
"\nExpected ",
expected[i],
"\nGot ",
result.ReadPortal().Get(i),
"\n");
}
}
void TestPointGradientUniform3D()
{
std::cout << "Testing PointGradient Worklet on 3D structured data" << std::endl;
vtkm::cont::testing::MakeTestDataSet testDataSet;
vtkm::cont::DataSet dataSet = testDataSet.Make3DUniformDataSet0();
vtkm::cont::ArrayHandle<vtkm::Float32> fieldArray;
dataSet.GetField("pointvar").GetData().AsArrayHandle(fieldArray);
vtkm::worklet::PointGradient gradient;
auto result = gradient.Run(dataSet.GetCellSet(), dataSet.GetCoordinateSystem(), fieldArray);
vtkm::Vec3f_32 expected[4] = {
{ 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)
{
VTKM_TEST_ASSERT(test_equal(result.ReadPortal().Get(i), expected[i]),
"Wrong result for PointGradient worklet on 3D uniform data",
"\nExpected ",
expected[i],
"\nGot ",
result.ReadPortal().Get(i),
"\n");
}
}
void TestPointGradientUniform3DWithVectorField()
{
std::cout << "Testing PointGradient Worklet with a vector field on 3D structured 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::Vec3f_64> 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::Vec3f_64> input =
vtkm::cont::make_ArrayHandle(vec, vtkm::CopyFlag::On);
vtkm::worklet::PointGradient gradient;
auto result = gradient.Run(dataSet.GetCellSet(), dataSet.GetCoordinateSystem(), input);
vtkm::Vec<vtkm::Vec3f_64, 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::Vec3f_64, 3> e = expected[i];
vtkm::Vec<vtkm::Vec3f_64, 3> r = result.ReadPortal().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 TestPointGradientUniform3DWithVectorField2()
{
std::cout << "Testing PointGradient Worklet with a vector field on 3D structured 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::Vec3f_64> 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::Vec3f_64> input =
vtkm::cont::make_ArrayHandle(vec, vtkm::CopyFlag::On);
vtkm::worklet::GradientOutputFields<vtkm::Vec3f_64> 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);
//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::Vec3f_64, 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::Vec3f_64, 3> eg = expected_gradients[i];
vtkm::Float64 d = extraOutput.Divergence.ReadPortal().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::Vec3f_64 ev(eg[1][2] - eg[2][1], eg[2][0] - eg[0][2], eg[0][1] - eg[1][0]);
vtkm::Vec3f_64 v = extraOutput.Vorticity.ReadPortal().Get(i);
VTKM_TEST_ASSERT(test_equal(ev, v), "Wrong result for Vorticity on 3D uniform data");
const vtkm::Vec3f_64 es(eg[1][2] + eg[2][1], eg[2][0] + eg[0][2], eg[0][1] + eg[1][0]);
const vtkm::Vec3f_64 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.ReadPortal().Get(i);
VTKM_TEST_ASSERT(
test_equal(qcriterion, q),
"Wrong result for QCriterion field of PointGradient worklet on 3D uniform data");
}
}
void TestPointGradientExplicit3D()
{
std::cout << "Testing PointGradient Worklet on Explicit 3D data" << std::endl;
vtkm::cont::testing::MakeTestDataSet testDataSet;
vtkm::cont::DataSet dataSet = testDataSet.Make3DExplicitDataSet5();
vtkm::cont::ArrayHandle<vtkm::Float32> fieldArray;
dataSet.GetField("pointvar").GetData().AsArrayHandle(fieldArray);
vtkm::worklet::PointGradient gradient;
auto result = gradient.Run(dataSet.GetCellSet(), dataSet.GetCoordinateSystem(), fieldArray);
//vtkm::cont::printSummary_ArrayHandle(result, std::cout, true);
const int nVerts = 11;
vtkm::Vec3f_32 expected[nVerts] = {
{ 10.0f, 40.2f, 30.1f }, { 27.4f, 40.1f, 10.1f }, { 17.425f, 40.0f, 10.1f },
{ -10.0f, 40.1f, 30.1f }, { 9.9f, -0.0500011f, 30.0f }, { 16.2125f, -4.55f, 10.0f },
{ 6.2f, -4.6f, 10.0f }, { -10.1f, -0.0999985f, 30.0f }, { 22.5125f, -4.575f, 10.025f },
{ 1.0f, -40.3f, 30.0f }, { 0.6f, -49.2f, 10.0f }
};
for (int i = 0; i < nVerts; ++i)
{
VTKM_TEST_ASSERT(test_equal(result.ReadPortal().Get(i), expected[i]),
"Wrong result for PointGradient worklet on 3D explicit data",
"\nExpected ",
expected[i],
"\nGot ",
result.ReadPortal().Get(i),
"\n");
}
}
void TestPointGradientExplicit2D()
{
std::cout << "Testing PointGradient Worklet on Explicit 2D data" << std::endl;
vtkm::cont::testing::MakeTestDataSet testDataSet;
vtkm::cont::DataSet dataSet = testDataSet.Make2DExplicitDataSet0();
vtkm::cont::ArrayHandle<vtkm::Float32> fieldArray;
dataSet.GetField("pointvar").GetData().AsArrayHandle(fieldArray);
vtkm::worklet::PointGradient gradient;
auto result = gradient.Run(dataSet.GetCellSet(), dataSet.GetCoordinateSystem(), fieldArray);
//vtkm::cont::printSummary_ArrayHandle(result, std::cout, true);
const int nVerts = 16;
vtkm::Vec3f_32 expected[nVerts] = {
{ -22.0f, -7.0f, 0.0f }, { -25.5f, -7.0f, 0.0f }, { -30.5f, 7.0f, 0.0f },
{ -32.0f, 16.0f, 0.0f }, { -23.0f, -42.0f, 0.0f }, { -23.25f, -17.0f, 0.0f },
{ -20.6667f, 1.33333f, 0.0f }, { -23.0f, 14.0f, 0.0f }, { -8.0f, -42.0f, 0.0f },
{ 2.91546f, -24.8357f, 0.0f }, { -0.140736f, -7.71853f, 0.0f }, { -5.0f, 12.0f, 0.0f },
{ 31.8803f, 1.0f, 0.0f }, { -44.8148f, 20.5f, 0.0f }, { 38.5653f, 5.86938f, 0.0f },
{ 26.3967f, 86.7934f, 0.0f }
};
for (int i = 0; i < nVerts; ++i)
{
VTKM_TEST_ASSERT(test_equal(result.ReadPortal().Get(i), expected[i]),
"Wrong result for PointGradient worklet on 2D explicit data",
"\nExpected ",
expected[i],
"\nGot ",
result.ReadPortal().Get(i),
"\n");
}
}
void TestPointGradient()
{
TestPointGradientUniform2D();
TestPointGradientUniform3D();
TestPointGradientUniform3DWithVectorField();
TestPointGradientUniform3DWithVectorField2();
TestPointGradientExplicit2D();
TestPointGradientExplicit3D();
}
}
int UnitTestPointGradient(int argc, char* argv[])
{
return vtkm::cont::testing::Testing::Run(TestPointGradient, argc, argv);
}