Shuenhoy/vtk-m
1
0
Fork 0
mirror of https://gitlab.kitware.com/vtk/vtk-m synced 2024-09-16 17:22:55 +00:00
Code Issues 2 Actions Packages Projects Releases Wiki Activity

Merge branch 'master' of https://gitlab.kitware.com/vtk/vtk-m into sl_block_duplication

Browse Source
This commit is contained in:
Dave Pugmire 2023-04-14 12:54:00 -04:00
commit fe35202de8
29 changed files with 1768 additions and 155 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
.gitlab/ci/crusher.yml
View File

@ -24,8 +24,8 @@
DefApps/default
craype-accel-amd-gfx90a
rocm/4.5.0
gcc/10
cmake/3.22
gcc/12
cmake/3.23
git
git-lfs
ninja
@ -79,7 +79,7 @@ test:crusher_gcc_hip:
dependencies:
- build:crusher_gcc_hip
variables:
SCHEDULER_PARAMETERS: "-ACSC331 -pbatch -t120 --nice=0 -c32 --gpus=4 -N 1"
SCHEDULER_PARAMETERS: "-ACSC331_crusher -pbatch -t120 --nice=0 -c32 --gpus=4 -N 1"
# We need this to skip ctest_submit from being run inside a jsrun job
GITLAB_CI_EMULATION: 1
# Tests errors to address due to different env in Spock

9
CONTRIBUTING.md
View File

@ -32,14 +32,19 @@ Before you begin, perform initial setup:
This will prompt for your GitLab user name and configure a remote
called `gitlab` to refer to it.
5. (Optional but highly recommended.)
5. (Optional but highly recommended.)
[Disabling git-lfs in your fork] is needed to add/modify git-lfs files.
Find the setting to disable git-lfs in your fork through your fork web UI:
Settings/General/Project Features/Git LFS; set it to off; and _save changes_.
6. (Optional but highly recommended.)
[Register with the VTK-m dashboard] on Kitware's CDash instance to
better know how your code performs in regression tests. After
registering and signing in, click on "All Dashboards" link in the upper
left corner, scroll down and click "Subscribe to this project" on the
right of VTK-m.
6. (Optional but highly recommended.)
7. (Optional but highly recommended.)
[Sign up for the VTK-m mailing list] to communicate with other
developers and users.

9
docs/changelog/clip-remove-unused-points.md Normal file
View File

@ -0,0 +1,9 @@
# Clip now doesn't copy unused points from the input to the output
Previously, clip would just copy all the points and point data from the input to the output,
and only append the new points. This would affect the bounds computation of the result.
If the caller wanted to remove the unused points, they had to run the CleanGrid filter
on the result.
With this change, clip now keeps track of which inputs are actually part of the output
and copies only those points.

7
docs/changelog/continuous-scatterplot.md Normal file
View File

@ -0,0 +1,7 @@
# Continuous Scatterplot filter
This new filter, designed for bi-variate analysis, computes the continuous scatter-plot of a 3D mesh for two given scalar point fields.
The continuous scatter-plot is an extension of the discrete scatter-plot for continuous bi-variate analysis. The output dataset consists of triangle-shaped cells, whose coordinates on the 2D plane represent respectively the values of both scalar fields. The point centered scalar field generated on the triangular mesh quantifies the density of values in the data domain.
This VTK-m implementation is based on the algorithm presented in the paper "Continuous Scatterplots" by S. Bachthaler and D. Weiskopf.

9
examples/contour_tree_augmented/ContourTreeApp.cxx
View File

@ -756,14 +756,13 @@ int main(int argc, char* argv[])
//----Isovalue seleciton start
if (numLevels > 0) // if compute isovalues
{
// Get the data values for computing the explicit branch decomposition
#ifdef WITH_MPI
// Get the data values for computing the explicit branch decomposition
vtkm::cont::ArrayHandle<ValueType> dataField;
result.GetPartitions()[0].GetField(0).GetData().AsArrayHandle(dataField);
#ifdef WITH_MPI
result.GetPartitions()[0].GetField("values").GetData().AsArrayHandle(dataField);
bool dataFieldIsSorted = true;
#else
vtkm::cont::ArrayHandle<ValueType> dataField;
useDataSet.GetField(0).GetData().AsArrayHandle(dataField);
useDataSet.GetField("values").GetData().AsArrayHandle(dataField);
bool dataFieldIsSorted = false;
#endif

4
vtkm/filter/contour/testing/UnitTestClipWithFieldFilter.cxx
View File

@ -70,8 +70,8 @@ void TestClipExplicit()
vtkm::cont::ArrayHandle<vtkm::Float32> resultArrayHandle;
temp.AsArrayHandle(resultArrayHandle);
vtkm::Float32 expected[7] = { 1, 2, 1, 0, 0.5, 0.5, 0.5 };
for (int i = 0; i < 7; ++i)
vtkm::Float32 expected[6] = { 1, 2, 1, 0.5, 0.5, 0.5 };
for (int i = 0; i < 6; ++i)
{
VTKM_TEST_ASSERT(test_equal(resultArrayHandle.ReadPortal().Get(i), expected[i]),
"Wrong result for Clip fliter on triangle explicit data");

12
vtkm/filter/contour/testing/UnitTestClipWithImplicitFunctionFilter.cxx
View File

@ -71,11 +71,11 @@ void TestClipStructured(vtkm::Float64 offset)
vtkm::cont::ArrayHandle<vtkm::Float32> resultArrayHandle;
temp.AsArrayHandle(resultArrayHandle);
VTKM_TEST_ASSERT(resultArrayHandle.GetNumberOfValues() == 13,
VTKM_TEST_ASSERT(resultArrayHandle.GetNumberOfValues() == 12,
"Wrong number of points in the output dataset");
vtkm::Float32 expected[13] = { 1, 1, 1, 1, 0, 1, 1, 1, 1, 0.25, 0.25, 0.25, 0.25 };
for (int i = 0; i < 13; ++i)
vtkm::Float32 expected[12] = { 1, 1, 1, 1, 1, 1, 1, 1, 0.25, 0.25, 0.25, 0.25 };
for (int i = 0; i < 12; ++i)
{
VTKM_TEST_ASSERT(test_equal(resultArrayHandle.ReadPortal().Get(i), expected[i]),
"Wrong result for ClipWithImplicitFunction fliter on sturctured quads data");
@ -107,11 +107,11 @@ void TestClipStructuredInverted()
vtkm::cont::ArrayHandle<vtkm::Float32> resultArrayHandle;
temp.AsArrayHandle(resultArrayHandle);
VTKM_TEST_ASSERT(resultArrayHandle.GetNumberOfValues() == 13,
VTKM_TEST_ASSERT(resultArrayHandle.GetNumberOfValues() == 5,
"Wrong number of points in the output dataset");
vtkm::Float32 expected[13] = { 1, 1, 1, 1, 0, 1, 1, 1, 1, 0.25, 0.25, 0.25, 0.25 };
for (int i = 0; i < 13; ++i)
vtkm::Float32 expected[5] = { 0, 0.25, 0.25, 0.25, 0.25 };
for (int i = 0; i < 5; ++i)
{
VTKM_TEST_ASSERT(test_equal(resultArrayHandle.ReadPortal().Get(i), expected[i]),
"Wrong result for ClipWithImplicitFunction fliter on sturctured quads data");

130
vtkm/filter/contour/worklet/Clip.h
View File

@ -10,10 +10,8 @@
#ifndef vtkm_m_worklet_Clip_h
#define vtkm_m_worklet_Clip_h
#include <vtkm/filter/clean_grid/worklet/RemoveUnusedPoints.h>
#include <vtkm/filter/contour/worklet/clip/ClipTables.h>
#include <vtkm/worklet/DispatcherMapField.h>
#include <vtkm/worklet/DispatcherMapTopology.h>
#include <vtkm/worklet/DispatcherReduceByKey.h>
#include <vtkm/worklet/Keys.h>
#include <vtkm/worklet/WorkletMapField.h>
#include <vtkm/worklet/WorkletMapTopology.h>
@ -26,6 +24,7 @@
#include <vtkm/cont/CellSetExplicit.h>
#include <vtkm/cont/ConvertNumComponentsToOffsets.h>
#include <vtkm/cont/CoordinateSystem.h>
#include <vtkm/cont/Invoker.h>
#include <vtkm/cont/Timer.h>
#include <vtkm/cont/UnknownArrayHandle.h>
@ -327,6 +326,7 @@ public:
FieldInCell clipStats,
ExecObject clipTables,
ExecObject connectivityObject,
WholeArrayOut pointsOnlyConnectivityIndices,
WholeArrayOut edgePointReverseConnectivity,
WholeArrayOut edgePointInterpolation,
WholeArrayOut inCellReverseConnectivity,
@ -351,7 +351,8 @@ public:
_11,
_12,
_13,
_14);
_14,
_15);
template <typename CellShapeTag,
typename PointVecType,
@ -368,6 +369,7 @@ public:
const ClipStats& clipStats,
const internal::ClipTables::DevicePortal<DeviceAdapter>& clippingData,
ConnectivityObject& connectivityObject,
IdArrayType& pointsOnlyConnectivityIndices,
IdArrayType& edgePointReverseConnectivity,
EdgeInterpolationPortalType& edgePointInterpolation,
IdArrayType& inCellReverseConnectivity,
@ -390,6 +392,8 @@ public:
// Start Indices to keep track of interpolation points for new cell.
vtkm::Id inCellInterpPointIndex = clipStats.NumberOfInCellInterpPoints;
vtkm::Id inCellEdgeInterpIndex = clipStats.NumberOfInCellEdgeIndices;
// Start index of connectivityPointsOnly
vtkm::Id pointsOnlyConnectivityIndicesIndex = connectivityIndex - edgeIndex - inCellIndex;
// Iterate over the shapes for the current cell and begin to fill connectivity.
vtkm::Id numberOfCells = clippingData.ValueAt(clipIndex++);
@ -473,6 +477,8 @@ public:
}
else if (entry >= 100) // existing vertex
{
pointsOnlyConnectivityIndices.Set(pointsOnlyConnectivityIndicesIndex++,
connectivityIndex);
connectivityObject.SetConnectivity(connectivityIndex++, points[entry - 100]);
}
else // case of a new edge point
@ -589,14 +595,15 @@ public:
vtkm::Float64 value,
bool invert)
{
vtkm::cont::Invoker invoke;
// Create the required output fields.
vtkm::cont::ArrayHandle<ClipStats> clipStats;
vtkm::cont::ArrayHandle<vtkm::Id> clipTableIndices;
ComputeStats statsWorklet(value, invert);
//Send this CellSet to process
vtkm::worklet::DispatcherMapTopology<ComputeStats> statsDispatcher(statsWorklet);
statsDispatcher.Invoke(cellSet, scalars, this->ClipTablesInstance, clipStats, clipTableIndices);
ComputeStats statsWorklet(value, invert);
invoke(statsWorklet, cellSet, scalars, this->ClipTablesInstance, clipStats, clipTableIndices);
ClipStats zero;
vtkm::cont::ArrayHandle<ClipStats> cellSetStats;
@ -612,6 +619,10 @@ public:
shapes, numberOfIndices, connectivity, offsets, total);
//Begin Process of Constructing the new CellSet.
vtkm::cont::ArrayHandle<vtkm::Id> pointsOnlyConnectivityIndices;
pointsOnlyConnectivityIndices.Allocate(total.NumberOfIndices - total.NumberOfEdgeIndices -
total.NumberOfInCellIndices);
vtkm::cont::ArrayHandle<vtkm::Id> edgePointReverseConnectivity;
edgePointReverseConnectivity.Allocate(total.NumberOfEdgeIndices);
vtkm::cont::ArrayHandle<EdgeInterpolation> edgeInterpolation;
@ -628,25 +639,53 @@ public:
this->InCellInterpolationInfo.Allocate(total.NumberOfInCellInterpPoints);
this->CellMapOutputToInput.Allocate(total.NumberOfCells);
GenerateCellSet cellSetWorklet(value);
//Send this CellSet to process
vtkm::worklet::DispatcherMapTopology<GenerateCellSet> cellSetDispatcher(cellSetWorklet);
cellSetDispatcher.Invoke(cellSet,
scalars,
clipTableIndices,
cellSetStats,
this->ClipTablesInstance,
connectivityObject,
edgePointReverseConnectivity,
edgeInterpolation,
cellPointReverseConnectivity,
cellPointEdgeReverseConnectivity,
cellPointEdgeInterpolation,
this->InCellInterpolationKeys,
this->InCellInterpolationInfo,
this->CellMapOutputToInput);
GenerateCellSet cellSetWorklet(value);
invoke(cellSetWorklet,
cellSet,
scalars,
clipTableIndices,
cellSetStats,
this->ClipTablesInstance,
connectivityObject,
pointsOnlyConnectivityIndices,
edgePointReverseConnectivity,
edgeInterpolation,
cellPointReverseConnectivity,
cellPointEdgeReverseConnectivity,
cellPointEdgeInterpolation,
this->InCellInterpolationKeys,
this->InCellInterpolationInfo,
this->CellMapOutputToInput);
this->InterpolationKeysBuilt = false;
clipTableIndices.ReleaseResources();
cellSetStats.ReleaseResources();
// extract only the used points from the input
{
vtkm::cont::ArrayHandle<vtkm::IdComponent> pointMask;
pointMask.AllocateAndFill(scalars.GetNumberOfValues(), 0);
auto pointsOnlyConnectivity =
vtkm::cont::make_ArrayHandlePermutation(pointsOnlyConnectivityIndices, connectivity);
invoke(
vtkm::worklet::RemoveUnusedPoints::GeneratePointMask{}, pointsOnlyConnectivity, pointMask);
vtkm::worklet::ScatterCounting scatter(pointMask, true);
auto pointMapInputToOutput = scatter.GetInputToOutputMap();
this->PointMapOutputToInput = scatter.GetOutputToInputMap();
pointMask.ReleaseResources();
invoke(vtkm::worklet::RemoveUnusedPoints::TransformPointIndices{},
pointsOnlyConnectivity,
pointMapInputToOutput,
pointsOnlyConnectivity);
pointsOnlyConnectivityIndices.ReleaseResources();
}
// Get unique EdgeInterpolation : unique edge points.
// LowerBound for edgeInterpolation : get index into new edge points array.
// LowerBound for cellPointEdgeInterpolation : get index into new edge points array.
@ -660,35 +699,37 @@ public:
edgeInterpolation,
edgeInterpolationIndexToUnique,
EdgeInterpolation::LessThanOp());
edgeInterpolation.ReleaseResources();
vtkm::cont::ArrayHandle<vtkm::Id> cellInterpolationIndexToUnique;
vtkm::cont::Algorithm::LowerBounds(this->EdgePointsInterpolation,
cellPointEdgeInterpolation,
cellInterpolationIndexToUnique,
EdgeInterpolation::LessThanOp());
cellPointEdgeInterpolation.ReleaseResources();
this->EdgePointsOffset = scalars.GetNumberOfValues();
this->EdgePointsOffset = this->PointMapOutputToInput.GetNumberOfValues();
this->InCellPointsOffset =
this->EdgePointsOffset + this->EdgePointsInterpolation.GetNumberOfValues();
// Scatter these values into the connectivity array,
// scatter indices are given in reverse connectivity.
ScatterEdgeConnectivity scatterEdgePointConnectivity(this->EdgePointsOffset);
vtkm::worklet::DispatcherMapField<ScatterEdgeConnectivity> scatterEdgeDispatcher(
scatterEdgePointConnectivity);
scatterEdgeDispatcher.Invoke(
edgeInterpolationIndexToUnique, edgePointReverseConnectivity, connectivity);
scatterEdgeDispatcher.Invoke(cellInterpolationIndexToUnique,
cellPointEdgeReverseConnectivity,
this->InCellInterpolationInfo);
invoke(scatterEdgePointConnectivity,
edgeInterpolationIndexToUnique,
edgePointReverseConnectivity,
connectivity);
invoke(scatterEdgePointConnectivity,
cellInterpolationIndexToUnique,
cellPointEdgeReverseConnectivity,
this->InCellInterpolationInfo);
// Add offset in connectivity of all new in-cell points.
ScatterInCellConnectivity scatterInCellPointConnectivity(this->InCellPointsOffset);
vtkm::worklet::DispatcherMapField<ScatterInCellConnectivity> scatterInCellDispatcher(
scatterInCellPointConnectivity);
scatterInCellDispatcher.Invoke(cellPointReverseConnectivity, connectivity);
invoke(scatterInCellPointConnectivity, cellPointReverseConnectivity, connectivity);
vtkm::cont::CellSetExplicit<> output;
vtkm::Id numberOfPoints = scalars.GetNumberOfValues() +
vtkm::Id numberOfPoints = this->PointMapOutputToInput.GetNumberOfValues() +
this->EdgePointsInterpolation.GetNumberOfValues() + total.NumberOfInCellPoints;
vtkm::cont::ConvertNumComponentsToOffsets(numberOfIndices, offsets);
@ -830,21 +871,25 @@ public:
this->InterpolationKeys.BuildArrays(this->InCellInterpolationKeys, KeysSortType::Unstable);
}
vtkm::Id numberOfOriginalValues = input.GetNumberOfValues();
vtkm::Id numberOfVertexPoints = this->PointMapOutputToInput.GetNumberOfValues();
vtkm::Id numberOfEdgePoints = this->EdgePointsInterpolation.GetNumberOfValues();
vtkm::Id numberOfInCellPoints = this->InterpolationKeys.GetUniqueKeys().GetNumberOfValues();
// Copy over the original values. They are still part of the output. (Unused points are
// not culled. Use CleanGrid for that.)
output.Allocate(numberOfOriginalValues + numberOfEdgePoints + numberOfInCellPoints);
vtkm::cont::Algorithm::CopySubRange(input, 0, numberOfOriginalValues, output);
output.Allocate(numberOfVertexPoints + numberOfEdgePoints + numberOfInCellPoints);
// Copy over the original values that are still part of the output.
vtkm::cont::Algorithm::CopySubRange(
vtkm::cont::make_ArrayHandlePermutation(this->PointMapOutputToInput, input),
0,
numberOfVertexPoints,
output);
// Interpolate all new points that lie on edges of the input mesh.
vtkm::cont::Invoker invoke;
invoke(PerformEdgeInterpolations{},
this->EdgePointsInterpolation,
input,
vtkm::cont::make_ArrayHandleView(output, numberOfOriginalValues, numberOfEdgePoints));
vtkm::cont::make_ArrayHandleView(output, numberOfVertexPoints, numberOfEdgePoints));
// Perform a gather on the output to get all the required values for calculation of centroids
// using the interpolation info array.
@ -854,7 +899,7 @@ public:
this->InterpolationKeys,
toReduceValues,
vtkm::cont::make_ArrayHandleView(
output, numberOfOriginalValues + numberOfEdgePoints, numberOfInCellPoints));
output, numberOfVertexPoints + numberOfEdgePoints, numberOfInCellPoints));
}
vtkm::cont::ArrayHandle<vtkm::Id> GetCellMapOutputToInput() const
@ -868,6 +913,7 @@ private:
vtkm::cont::ArrayHandle<vtkm::Id> InCellInterpolationKeys;
vtkm::cont::ArrayHandle<vtkm::Id> InCellInterpolationInfo;
vtkm::cont::ArrayHandle<vtkm::Id> CellMapOutputToInput;
vtkm::cont::ArrayHandle<vtkm::Id> PointMapOutputToInput;
vtkm::Id EdgePointsOffset;
vtkm::Id InCellPointsOffset;
vtkm::worklet::Keys<vtkm::Id> InterpolationKeys;

2
vtkm/filter/density_estimate/CMakeLists.txt
View File

@ -8,6 +8,7 @@
## PURPOSE. See the above copyright notice for more information.
##============================================================================
set(density_estimate_headers
ContinuousScatterPlot.h
Entropy.h
Histogram.h
NDEntropy.h
@ -19,6 +20,7 @@ set(density_estimate_headers
)
set(density_estimate_sources_device
ContinuousScatterPlot.cxx
Entropy.cxx
Histogram.cxx
NDEntropy.cxx

82
vtkm/filter/density_estimate/ContinuousScatterPlot.cxx Normal file
View File

@ -0,0 +1,82 @@
//============================================================================
// 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/ArrayCopy.h>
#include <vtkm/cont/CoordinateSystem.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/ErrorFilterExecution.h>
#include <vtkm/filter/density_estimate/ContinuousScatterPlot.h>
#include <vtkm/filter/density_estimate/worklet/ContinuousScatterPlot.h>
#include <vtkm/filter/geometry_refinement/Tetrahedralize.h>
namespace vtkm
{
namespace filter
{
namespace density_estimate
{
VTKM_CONT vtkm::cont::DataSet ContinuousScatterPlot::DoExecute(const vtkm::cont::DataSet& input)
{
// This algorithm only operate on tetra cells, we need to apply the tetrahedralize filter first.
auto tetrahedralizeFilter = vtkm::filter::geometry_refinement::Tetrahedralize();
auto tetraInput = tetrahedralizeFilter.Execute(input);
vtkm::cont::CellSetSingleType<> tetraCellSet;
tetraInput.GetCellSet().AsCellSet(tetraCellSet);
vtkm::cont::Field scalarField1 = input.GetField(GetActiveFieldName(0));
vtkm::cont::Field scalarField2 = input.GetField(GetActiveFieldName(1));
if (!(scalarField1.IsPointField() && scalarField2.IsPointField()))
{
throw vtkm::cont::ErrorFilterExecution("Point fields expected.");
}
const auto& coords = input.GetCoordinateSystem().GetDataAsMultiplexer();
vtkm::cont::CoordinateSystem activeCoordSystem = input.GetCoordinateSystem();
vtkm::cont::DataSet scatterplotDataSet;
vtkm::worklet::ContinuousScatterPlot worklet;
auto resolveFieldType = [&](const auto& resolvedScalar) {
using FieldType = typename std::decay_t<decltype(resolvedScalar)>::ValueType;
vtkm::cont::CellSetSingleType<> scatterplotCellSet;
vtkm::cont::ArrayHandle<FieldType> density;
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>> newCoords;
// Both fields need to resolve to the same type in order to perform calculations
vtkm::cont::ArrayHandle<FieldType> resolvedScalar2;
vtkm::cont::ArrayCopyShallowIfPossible(scalarField2.GetData(), resolvedScalar2);
worklet.Run(tetraCellSet,
coords,
newCoords,
density,
resolvedScalar,
resolvedScalar2,
scatterplotCellSet);
// Populate the new dataset representing the continuous scatterplot
// Using the density field and coordinates calculated by the worklet
activeCoordSystem = vtkm::cont::CoordinateSystem(activeCoordSystem.GetName(), newCoords);
scatterplotDataSet.AddCoordinateSystem(activeCoordSystem);
scatterplotDataSet.SetCellSet(scatterplotCellSet);
scatterplotDataSet.AddPointField(this->GetOutputFieldName(), density);
};
this->CastAndCallScalarField(scalarField1.GetData(), resolveFieldType);
return scatterplotDataSet;
}
} // namespace density_estimate
} // namespace filter
} // namespace vtkm

60
vtkm/filter/density_estimate/ContinuousScatterPlot.h Normal file
View 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.
//============================================================================
#ifndef vtk_m_filter_density_estimate_ContinuousScatterPlot_h
#define vtk_m_filter_density_estimate_ContinuousScatterPlot_h
#include <vtkm/filter/FilterField.h>
#include <vtkm/filter/density_estimate/vtkm_filter_density_estimate_export.h>
namespace vtkm
{
namespace filter
{
namespace density_estimate
{
/// \brief Constructs the continuous scatterplot for two given scalar point fields of a mesh.
///
/// The continuous scatterplot is an extension of the discrete scatterplot for continuous bi-variate analysis.
/// This filter outputs an ExplicitDataSet of triangle-shaped cells, whose coordinates on the 2D plane represent respectively
/// the values of both scalar fields. Triangles' points are associated with a scalar field, representing the
/// density of values in the data domain. The filter tetrahedralizes the input dataset before operating.
///
/// If both fields provided don't have the same floating point precision, the output will
/// have the precision of the first one of the pair.
///
/// This implementation is based on the algorithm presented in the publication :
///
/// S. Bachthaler and D. Weiskopf, "Continuous Scatterplots"
/// in IEEE Transactions on Visualization and Computer Graphics,
/// vol. 14, no. 6, pp. 1428-1435, Nov.-Dec. 2008
/// doi: 10.1109/TVCG.2008.119.
class VTKM_FILTER_DENSITY_ESTIMATE_EXPORT ContinuousScatterPlot : public vtkm::filter::FilterField
{
public:
VTKM_CONT ContinuousScatterPlot() { this->SetOutputFieldName("density"); }
/// Select both point fields to use when running the filter.
/// Replaces setting each one individually using `SetActiveField` on indices 0 and 1.
VTKM_CONT
void SetActiveFieldsPair(const std::string& fieldName1, const std::string& fieldName2)
{
SetActiveField(0, fieldName1, vtkm::cont::Field::Association::Points);
SetActiveField(1, fieldName2, vtkm::cont::Field::Association::Points);
};
private:
VTKM_CONT vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet& input) override;
};
} // namespace density_estimate
} // namespace filter
} // namespace vtm
#endif //vtk_m_filter_density_estimate_ContinuousScatterPlot_h

1
vtkm/filter/density_estimate/testing/CMakeLists.txt
View File

@ -9,6 +9,7 @@
##============================================================================
set(unit_tests
UnitTestContinuousScatterPlot.cxx
UnitTestEntropyFilter.cxx
UnitTestHistogramFilter.cxx
UnitTestNDEntropyFilter.cxx

690
vtkm/filter/density_estimate/testing/UnitTestContinuousScatterPlot.cxx Normal file
View File

@ -0,0 +1,690 @@
//============================================================================
// 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/ArrayHandle.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/DataSetBuilderExplicit.h>
#include <vtkm/cont/ErrorFilterExecution.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/filter/density_estimate/ContinuousScatterPlot.h>
namespace
{
std::vector<vtkm::Vec3f> tetraCoords{ vtkm::Vec3f(0.0f, 0.0f, 0.0f),
vtkm::Vec3f(2.0f, 0.0f, 0.0f),
vtkm::Vec3f(2.0f, 2.0f, 0.0f),
vtkm::Vec3f(1.0f, 0.0f, 2.0f) };
std::vector<vtkm::UInt8> tetraShape{ vtkm::CELL_SHAPE_TETRA };
std::vector<vtkm::IdComponent> tetraIndex{ 4 };
std::vector<vtkm::Id> tetraConnectivity{ 0, 1, 2, 3 };
std::vector<vtkm::Vec3f> multiCoords{
vtkm::Vec3f(0.0f, 0.0f, 0.0f), vtkm::Vec3f(2.0f, 0.0f, 0.0f), vtkm::Vec3f(2.0f, 2.0f, 0.0f),
vtkm::Vec3f(1.0f, 0.0f, 2.0f), vtkm::Vec3f(0.0f, 0.0f, 0.0f), vtkm::Vec3f(2.0f, 0.0f, 0.0f),
vtkm::Vec3f(2.0f, 2.0f, 0.0f), vtkm::Vec3f(1.0f, 0.0f, 2.0f), vtkm::Vec3f(0.0f, 0.0f, 0.0f),
vtkm::Vec3f(2.0f, 0.0f, 0.0f), vtkm::Vec3f(2.0f, 2.0f, 0.0f), vtkm::Vec3f(1.0f, 0.0f, 2.0f)
};
std::vector<vtkm::UInt8> multiShapes{ vtkm::CELL_SHAPE_TETRA,
vtkm::CELL_SHAPE_TETRA,
vtkm::CELL_SHAPE_TETRA };
std::vector<vtkm::IdComponent> multiIndices{ 4, 4, 4 };
std::vector<vtkm::Id> multiConnectivity{ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 };
template <typename FieldType1, typename FieldType2>
vtkm::cont::DataSet makeDataSet(const std::vector<vtkm::Vec3f>& ds_coords,
const std::vector<vtkm::UInt8>& ds_shapes,
const std::vector<vtkm::IdComponent>& ds_indices,
const std::vector<vtkm::Id>& ds_connectivity,
const FieldType1& scalar1,
const FieldType2& scalar2)
{
vtkm::cont::DataSetBuilderExplicit dsb;
vtkm::cont::DataSet ds = dsb.Create(ds_coords, ds_shapes, ds_indices, ds_connectivity);
ds.AddPointField("scalar1", scalar1, 4);
ds.AddPointField("scalar2", scalar2, 4);
return ds;
}
vtkm::cont::DataSet executeFilter(const vtkm::cont::DataSet ds)
{
auto continuousSCP = vtkm::filter::density_estimate::ContinuousScatterPlot();
continuousSCP.SetActiveFieldsPair("scalar1", "scalar2");
return continuousSCP.Execute(ds);
}
template <typename PositionsPortalType, typename FieldType1, typename FieldType2>
void testCoordinates(const PositionsPortalType& positionsP,
const FieldType1& scalar1,
const FieldType2& scalar2,
const vtkm::IdComponent numberOfPoints)
{
for (vtkm::IdComponent i = 0; i < numberOfPoints; i++)
{
VTKM_TEST_ASSERT(test_equal(positionsP.Get(i)[0], scalar1[i]), "Wrong point coordinates");
VTKM_TEST_ASSERT(test_equal(positionsP.Get(i)[1], scalar2[i]), "Wrong point coordinates");
VTKM_TEST_ASSERT(test_equal(positionsP.Get(i)[2], 0.0f),
"Z coordinate value in the scatterplot should always be null");
}
}
template <typename DensityArrayType, typename FieldType>
void testDensity(const DensityArrayType& density,
const vtkm::IdComponent& centerId,
const FieldType& centerDensity)
{
for (vtkm::IdComponent i = 0; i < density.GetNumberOfValues(); i++)
{
if (i == centerId)
{
VTKM_TEST_ASSERT(test_equal(density.Get(i), centerDensity),
"Wrong density in the middle point of the cell");
}
else
{
VTKM_TEST_ASSERT(test_equal(density.Get(i), 0.0f),
"Density on the edge of the tetrahedron should be null");
}
}
}
template <typename CellSetType>
void testShapes(const CellSetType& cellSet)
{
for (vtkm::IdComponent i = 0; i < cellSet.GetNumberOfCells(); i++)
{
VTKM_TEST_ASSERT(cellSet.GetCellShape(i) == vtkm::CELL_SHAPE_TRIANGLE);
}
}
template <typename CellSetType>
void testConnectivity(const CellSetType& cellSet,
const vtkm::cont::ArrayHandle<vtkm::IdComponent>& expectedConnectivityArray)
{
VTKM_TEST_ASSERT(
test_equal_ArrayHandles(
cellSet.GetConnectivityArray(vtkm::TopologyElementTagCell(), vtkm::TopologyElementTagPoint()),
expectedConnectivityArray),
"Wrong connectivity");
}
void TestSingleTetraProjectionQuadConvex()
{
constexpr vtkm::FloatDefault scalar1[4] = {
0.0f,
1.0f,
0.0f,
-2.0f,
};
constexpr vtkm::FloatDefault scalar2[4] = {
-1.0f,
0.0f,
2.0f,
0.0f,
};
vtkm::cont::DataSet ds =
makeDataSet(tetraCoords, tetraShape, tetraIndex, tetraConnectivity, scalar1, scalar2);
auto scatterPlot = executeFilter(ds);
VTKM_TEST_ASSERT(test_equal(scatterPlot.GetNumberOfCells(), 4),
"Wrong number of projected triangles in the continuous scatter plot");
VTKM_TEST_ASSERT(test_equal(scatterPlot.GetNumberOfPoints(), 5),
"Wrong number of projected points in the continuous scatter plot");
// Test point positions
auto positions = scatterPlot.GetCoordinateSystem().GetDataAsMultiplexer();
auto positionsP = positions.ReadPortal();
testCoordinates(positionsP, scalar1, scalar2, 4);
// Test intersection point coordinates
VTKM_TEST_ASSERT(test_equal(positionsP.Get(4), vtkm::Vec3f{ 0.0f, 0.0f, 0.0f }),
"Wrong intersection point coordinates");
// Test for triangle shapes
vtkm::cont::CellSetSingleType<> cellSet;
scatterPlot.GetCellSet().AsCellSet(cellSet);
testShapes(cellSet);
// Test connectivity
auto expectedConnectivityArray =
vtkm::cont::make_ArrayHandle<vtkm::IdComponent>({ 0, 1, 4, 1, 2, 4, 2, 3, 4, 3, 0, 4 });
testConnectivity(cellSet, expectedConnectivityArray);
// Test density values
auto density = scatterPlot.GetField("density")
.GetData()
.AsArrayHandle<vtkm::cont::ArrayHandle<vtkm::FloatDefault>>()
.ReadPortal();
testDensity(density, 4, 0.888889f);
}
void TestSingleTetraProjectionQuadSelfIntersect()
{
constexpr vtkm::FloatDefault scalar1[4] = {
0.0f,
0.0f,
1.0f,
-2.0f,
};
constexpr vtkm::FloatDefault scalar2[4] = {
-1.0f,
2.0f,
0.0f,
0.0f,
};
vtkm::cont::DataSet ds =
makeDataSet(tetraCoords, tetraShape, tetraIndex, tetraConnectivity, scalar1, scalar2);
auto scatterPlot = executeFilter(ds);
VTKM_TEST_ASSERT(test_equal(scatterPlot.GetNumberOfCells(), 4),
"Wrong number of projected triangles in the continuous scatter plot");
// Test point positions
auto positions = scatterPlot.GetCoordinateSystem().GetDataAsMultiplexer();
auto positionsP = positions.ReadPortal();
testCoordinates(positionsP, scalar1, scalar2, 4);
// Test connectivity
vtkm::cont::CellSetSingleType<> cellSet;
scatterPlot.GetCellSet().AsCellSet(cellSet);
auto expectedConnectivityArray =
vtkm::cont::make_ArrayHandle<vtkm::IdComponent>({ 0, 2, 4, 2, 1, 4, 1, 3, 4, 3, 0, 4 });
testConnectivity(cellSet, expectedConnectivityArray);
}
void TestSingleTetraProjectionQuadInverseOrder()
{
constexpr vtkm::FloatDefault scalar1[4] = {
-2.0f,
0.0f,
1.0f,
0.0f,
};
constexpr vtkm::FloatDefault scalar2[4] = {
0.0f,
2.0f,
0.0f,
-1.0f,
};
vtkm::cont::DataSet ds =
makeDataSet(tetraCoords, tetraShape, tetraIndex, tetraConnectivity, scalar1, scalar2);
auto scatterPlot = executeFilter(ds);
VTKM_TEST_ASSERT(test_equal(scatterPlot.GetNumberOfCells(), 4),
"Wrong number of projected triangles in the continuous scatter plot");
// Test connectivity
vtkm::cont::CellSetSingleType<> cellSet;
scatterPlot.GetCellSet().AsCellSet(cellSet);
auto expectedConnectivityArray = vtkm::cont::make_ArrayHandle<vtkm::IdComponent>(
{ 0,
1,
4,
1,
2,
4,
2,
3,
4,
3,
0,
4 }); // Inverting the order of points should not change connectivity
testConnectivity(cellSet, expectedConnectivityArray);
}
void TestSingleTetraProjectionQuadSelfIntersectSecond()
{
constexpr vtkm::FloatDefault scalar1[4] = {
0.0f,
1.0f,
-2.0f,
0.0f,
};
constexpr vtkm::FloatDefault scalar2[4] = {
-1.0f,
0.0f,
0.0f,
2.0f,
};
vtkm::cont::DataSet ds =
makeDataSet(tetraCoords, tetraShape, tetraIndex, tetraConnectivity, scalar1, scalar2);
auto scatterPlot = executeFilter(ds);
VTKM_TEST_ASSERT(test_equal(scatterPlot.GetNumberOfCells(), 4),
"Wrong number of projected triangles in the continuous scatter plot");
// Test connectivity
vtkm::cont::CellSetSingleType<> cellSet;
scatterPlot.GetCellSet().AsCellSet(cellSet);
auto expectedConnectivityArray =
vtkm::cont::make_ArrayHandle<vtkm::IdComponent>({ 0, 2, 4, 2, 3, 4, 3, 1, 4, 1, 0, 4 });
testConnectivity(cellSet, expectedConnectivityArray);
}
void TestSingleTetra_projection_triangle_point0Inside()
{
constexpr vtkm::FloatDefault scalar1[4] = {
3.0f,
3.0f,
4.0f,
1.0f,
};
constexpr vtkm::FloatDefault scalar2[4] = {
1.0f,
0.0f,
2.0f,
2.0f,
};
vtkm::cont::DataSet ds =
makeDataSet(tetraCoords, tetraShape, tetraIndex, tetraConnectivity, scalar1, scalar2);
auto scatterPlot = executeFilter(ds);
VTKM_TEST_ASSERT(test_equal(scatterPlot.GetNumberOfCells(), 3),
"Wrong number of projected triangles in the continuous scatter plot");
VTKM_TEST_ASSERT(test_equal(scatterPlot.GetNumberOfPoints(), 4),
"Wrong number of projected points in the continuous scatter plot");
// Test point positions
auto positions = scatterPlot.GetCoordinateSystem().GetDataAsMultiplexer();
auto positionsP = positions.ReadPortal();
testCoordinates(positionsP, scalar1, scalar2, 3);
// Test for triangle shapes
vtkm::cont::CellSetSingleType<> cellSet;
scatterPlot.GetCellSet().AsCellSet(cellSet);
testShapes(cellSet);
// Test connectivity
auto expectedConnectivityArray =
vtkm::cont::make_ArrayHandle<vtkm::IdComponent>({ 1, 2, 0, 2, 3, 0, 3, 1, 0 });
testConnectivity(cellSet, expectedConnectivityArray);
// Test density values
auto density = scatterPlot.GetField("density")
.GetData()
.AsArrayHandle<vtkm::cont::ArrayHandle<vtkm::FloatDefault>>()
.ReadPortal();
testDensity(density, 0, 1.33333f);
}
void TestSingleTetra_projection_triangle_point1Inside()
{
constexpr vtkm::FloatDefault scalar1[4] = {
3.0f,
3.0f,
4.0f,
1.0f,
};
constexpr vtkm::FloatDefault scalar2[4] = {
0.0f,
1.0f,
2.0f,
2.0f,
};
vtkm::cont::DataSet ds =
makeDataSet(tetraCoords, tetraShape, tetraIndex, tetraConnectivity, scalar1, scalar2);
auto scatterPlot = executeFilter(ds);
// Test connectivity
vtkm::cont::CellSetSingleType<> cellSet;
scatterPlot.GetCellSet().AsCellSet(cellSet);
auto expectedConnectivityArray =
vtkm::cont::make_ArrayHandle<vtkm::IdComponent>({ 0, 2, 1, 2, 3, 1, 3, 0, 1 });
testConnectivity(cellSet, expectedConnectivityArray);
// Test density values
auto density = scatterPlot.GetField("density")
.GetData()
.AsArrayHandle<vtkm::cont::ArrayHandle<vtkm::FloatDefault>>()
.ReadPortal();
testDensity(density, 1, 1.33333f);
}
void TestSingleTetra_projection_triangle_point2Inside()
{
constexpr vtkm::FloatDefault scalar1[4] = {
3.0f,
4.0f,
3.0f,
1.0f,
};
constexpr vtkm::FloatDefault scalar2[4] = {
0.0f,
2.0f,
1.0f,
2.0f,
};
vtkm::cont::DataSet ds =
makeDataSet(tetraCoords, tetraShape, tetraIndex, tetraConnectivity, scalar1, scalar2);
auto scatterPlot = executeFilter(ds);
// Test connectivity
vtkm::cont::CellSetSingleType<> cellSet;
scatterPlot.GetCellSet().AsCellSet(cellSet);
auto expectedConnectivityArray =
vtkm::cont::make_ArrayHandle<vtkm::IdComponent>({ 0, 1, 2, 1, 3, 2, 3, 0, 2 });
testConnectivity(cellSet, expectedConnectivityArray);
// Test density values
auto density = scatterPlot.GetField("density")
.GetData()
.AsArrayHandle<vtkm::cont::ArrayHandle<vtkm::FloatDefault>>()
.ReadPortal();
testDensity(density, 2, 1.33333f);
}
void TestSingleTetra_projection_triangle_point3Inside()
{
constexpr vtkm::FloatDefault scalar1[4] = {
3.0f,
4.0f,
1.0f,
3.0f,
};
constexpr vtkm::FloatDefault scalar2[4] = {
0.0f,
2.0f,
2.0f,
1.0f,
};
vtkm::cont::DataSet ds =
makeDataSet(tetraCoords, tetraShape, tetraIndex, tetraConnectivity, scalar1, scalar2);
auto scatterPlot = executeFilter(ds);
// Test connectivity
vtkm::cont::CellSetSingleType<> cellSet;
scatterPlot.GetCellSet().AsCellSet(cellSet);
auto expectedConnectivityArray =
vtkm::cont::make_ArrayHandle<vtkm::IdComponent>({ 0, 1, 3, 1, 2, 3, 2, 0, 3 });
testConnectivity(cellSet, expectedConnectivityArray);
// Test density values
auto density = scatterPlot.GetField("density")
.GetData()
.AsArrayHandle<vtkm::cont::ArrayHandle<vtkm::FloatDefault>>()
.ReadPortal();
testDensity(density, 3, 1.33333f);
}
void TestNullSpatialVolume()
{
constexpr vtkm::FloatDefault scalar1[4] = {
3.0f,
3.0f,
4.0f,
1.0f,
};
constexpr vtkm::FloatDefault scalar2[4] = {
1.0f,
0.0f,
2.0f,
2.0f,
};
std::vector<vtkm::Vec3f> nullCoords{ vtkm::Vec3f(1.0f, 1.0f, 1.0f),
vtkm::Vec3f(1.0f, 1.0f, 1.0f),
vtkm::Vec3f(1.0f, 1.0f, 1.0f),
vtkm::Vec3f(1.0f, 1.0f, 1.0f) };
vtkm::cont::DataSet ds =
makeDataSet(nullCoords, tetraShape, tetraIndex, tetraConnectivity, scalar1, scalar2);
auto scatterPlot = executeFilter(ds);
// Test density values
auto density = scatterPlot.GetField("density")
.GetData()
.AsArrayHandle<vtkm::cont::ArrayHandle<vtkm::FloatDefault>>()
.ReadPortal();
VTKM_TEST_ASSERT(vtkm::IsInf(density.Get(0)), "Density should be infinite for a null volume");
}
void TestNullDataSurface()
{
constexpr vtkm::FloatDefault scalar1[4] = {
0.0f,
1.0f,
3.0f,
2.0f,
};
constexpr vtkm::FloatDefault scalar2[4] = {
0.0f,
1.0f,
3.0f,
2.0f,
};
vtkm::cont::DataSet ds =
makeDataSet(tetraCoords, tetraShape, tetraIndex, tetraConnectivity, scalar1, scalar2);
auto scatterPlot = executeFilter(ds);
auto density = scatterPlot.GetField("density")
.GetData()
.AsArrayHandle<vtkm::cont::ArrayHandle<vtkm::FloatDefault>>()
.ReadPortal();
VTKM_TEST_ASSERT(vtkm::IsInf(density.Get(4)), "Density should be infinite for a null volume");
}
void TestMultipleTetra()
{
constexpr vtkm::FloatDefault multiscalar1[12] = { 3.0f, 3.0f, 4.0f, 1.0f, 0.0f, 1.0f,
0.0f, -2.0f, 3.0f, 3.0f, 4.0f, 1.0f };
constexpr vtkm::FloatDefault multiscalar2[12] = { 1.0f, 0.0f, 2.0f, 2.0f, -1.0f, 0.0f,
2.0f, 0.0f, 1.0f, 0.0f, 2.0f, 2.0f };
vtkm::cont::DataSetBuilderExplicit dsb;
vtkm::cont::DataSet ds = dsb.Create(multiCoords, multiShapes, multiIndices, multiConnectivity);
ds.AddPointField("scalar1", multiscalar1, 12);
ds.AddPointField("scalar2", multiscalar2, 12);
// Filtering
auto scatterPlot = executeFilter(ds);
VTKM_TEST_ASSERT(test_equal(scatterPlot.GetNumberOfCells(), 10),
"Wrong number of projected triangles in the continuous scatter plot");
VTKM_TEST_ASSERT(test_equal(scatterPlot.GetNumberOfPoints(), 13),
"Wrong number of projected points in the continuous scatter plot");
vtkm::cont::CellSetSingleType<> cellSet;
scatterPlot.GetCellSet().AsCellSet(cellSet);
testShapes(cellSet);
}
void TestNonTetra()
{
std::vector<vtkm::Vec3f> wedgeCoords{
vtkm::Vec3f(0.0f, 0.0f, 0.0f), vtkm::Vec3f(2.0f, 0.0f, 0.0f), vtkm::Vec3f(2.0f, 4.0f, 0.0f),
vtkm::Vec3f(0.0f, 4.0f, 0.0f), vtkm::Vec3f(1.0f, 0.0f, 3.0f), vtkm::Vec3f(1.0f, 4.0f, 3.0f),
};
constexpr vtkm::FloatDefault scalar1[6] = { 0.0f, 3.0f, 3.0f, 2.0f, 2.0f, 1.0f };
constexpr vtkm::FloatDefault scalar2[6] = { 0.0f, 1.0f, 3.0f, 2.0f, 0.0f, 1.0f };
std::vector<vtkm::UInt8> w_shape{ vtkm::CELL_SHAPE_WEDGE };
std::vector<vtkm::IdComponent> w_indices{ 6 };
std::vector<vtkm::Id> w_connectivity{ 0, 1, 2, 3, 4, 5 };
vtkm::cont::DataSetBuilderExplicit dsb;
vtkm::cont::DataSet ds = dsb.Create(wedgeCoords, w_shape, w_indices, w_connectivity);
ds.AddPointField("scalar1", scalar1, 6);
ds.AddPointField("scalar2", scalar2, 6);
auto scatterPlot = executeFilter(ds);
VTKM_TEST_ASSERT(test_equal(scatterPlot.GetNumberOfCells(), 12),
"Wrong number of projected triangles in the continuous scatter plot");
VTKM_TEST_ASSERT(test_equal(scatterPlot.GetNumberOfPoints(), 15),
"Wrong number of projected points in the continuous scatter plot");
}
void TestNonPointFields()
{
constexpr vtkm::FloatDefault cellField1[1] = { 1.0f };
constexpr vtkm::FloatDefault cellField2[1] = {
0.0f,
};
vtkm::cont::DataSetBuilderExplicit dsb;
vtkm::cont::DataSet ds = dsb.Create(tetraCoords, tetraShape, tetraIndex, tetraConnectivity);
ds.AddCellField("scalar1", cellField1, 1);
ds.AddCellField("scalar2", cellField2, 1);
auto continuousSCP = vtkm::filter::density_estimate::ContinuousScatterPlot();
continuousSCP.SetActiveField(0, "scalar1", vtkm::cont::Field::Association::Cells);
continuousSCP.SetActiveField(1, "scalar2", vtkm::cont::Field::Association::Cells);
try
{
auto scatterPlot = continuousSCP.Execute(ds);
VTKM_TEST_FAIL(
"Filter execution was not aborted after providing active fields not associated with points");
}
catch (const vtkm::cont::ErrorFilterExecution&)
{
std::cout << "Execution successfully aborted" << std::endl;
}
}
void TestDataTypes()
{
constexpr vtkm::Float32 scalar1_f32[4] = {
-2.0f,
0.0f,
1.0f,
0.0f,
};
constexpr vtkm::Float32 scalar2_f32[4] = {
0.0f,
2.0f,
0.0f,
-1.0f,
};
constexpr vtkm::Float64 scalar1_f64[4] = {
-2.0,
0.0,
1.0,
0.0,
};
constexpr vtkm::Float64 scalar2_f64[4] = {
0.0,
2.0,
0.0,
-1.0,
};
// The filter should run whatever float precision we use for both fields
vtkm::cont::DataSet ds =
makeDataSet(tetraCoords, tetraShape, tetraIndex, tetraConnectivity, scalar1_f32, scalar2_f32);
auto scatterPlot = executeFilter(ds);
auto density = scatterPlot.GetField("density")
.GetData()
.AsArrayHandle<vtkm::cont::ArrayHandle<vtkm::Float32>>()
.ReadPortal();
testDensity(density, 4, 0.888889f);
vtkm::cont::DataSet ds2 =
makeDataSet(tetraCoords, tetraShape, tetraIndex, tetraConnectivity, scalar1_f32, scalar2_f64);
auto scatterPlot2 = executeFilter(ds);
auto density2 = scatterPlot.GetField("density")
.GetData()
.AsArrayHandle<vtkm::cont::ArrayHandle<vtkm::Float32>>()
.ReadPortal();
testDensity(density2, 4, 0.888889f);
vtkm::cont::DataSet ds3 =
makeDataSet(tetraCoords, tetraShape, tetraIndex, tetraConnectivity, scalar1_f64, scalar2_f32);
auto scatterPlot3 = executeFilter(ds);
auto density3 = scatterPlot.GetField("density")
.GetData()
.AsArrayHandle<vtkm::cont::ArrayHandle<vtkm::Float32>>()
.ReadPortal();
testDensity(density3, 4, 0.888889f);
vtkm::cont::DataSet ds4 =
makeDataSet(tetraCoords, tetraShape, tetraIndex, tetraConnectivity, scalar1_f64, scalar2_f64);
auto scatterPlot4 = executeFilter(ds);
auto density4 = scatterPlot.GetField("density")
.GetData()
.AsArrayHandle<vtkm::cont::ArrayHandle<vtkm::Float32>>()
.ReadPortal();
testDensity(density4, 4, 0.888889f);
}
void TestContinuousScatterPlot()
{
// Projection forms 4 triangles in the data domain
TestSingleTetraProjectionQuadConvex();
TestSingleTetraProjectionQuadSelfIntersect();
TestSingleTetraProjectionQuadInverseOrder();
TestSingleTetraProjectionQuadSelfIntersectSecond();
// 3 triangles in the data domain
TestSingleTetra_projection_triangle_point0Inside();
TestSingleTetra_projection_triangle_point1Inside();
TestSingleTetra_projection_triangle_point2Inside();
TestSingleTetra_projection_triangle_point3Inside();
// // Edge cases
TestNullSpatialVolume();
TestNullDataSurface();
TestMultipleTetra();
TestNonTetra();
TestNonPointFields();
TestDataTypes();
}
}
int UnitTestContinuousScatterPlot(int argc, char* argv[])
{
return vtkm::cont::testing::Testing::Run(TestContinuousScatterPlot, argc, argv);
}

1
vtkm/filter/density_estimate/vtkm.module
View File

@ -6,6 +6,7 @@ DEPENDS
vtkm_filter_core
PRIVATE_DEPENDS
vtkm_worklet
vtkm_filter_geometry_refinement
OPTIONAL_DEPENDS
MPI::MPI_CXX
TEST_DEPENDS

1
vtkm/filter/density_estimate/worklet/CMakeLists.txt
View File

@ -9,6 +9,7 @@
##============================================================================
set(headers
ContinuousScatterPlot.h
FieldEntropy.h
FieldHistogram.h
NDimsEntropy.h

473
vtkm/filter/density_estimate/worklet/ContinuousScatterPlot.h Normal file
View File

@ -0,0 +1,473 @@
//============================================================================
// 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_worklet_ContinuousScatterPlot_h
#define vtk_m_worklet_ContinuousScatterPlot_h
#include <vtkm/VectorAnalysis.h>
#include <vtkm/cont/Algorithm.h>
#include <vtkm/cont/ArrayHandleGroupVec.h>
#include <vtkm/cont/ArrayHandleGroupVecVariable.h>
#include <vtkm/cont/ConvertNumComponentsToOffsets.h>
#include <vtkm/cont/Invoker.h>
#include <vtkm/worklet/ScatterCounting.h>
#include <vtkm/worklet/WorkletMapTopology.h>
namespace vtkm
{
namespace worklet
{
class ContinuousScatterPlot
{
public:
template <typename FieldType>
struct ClassifyTetra : vtkm::worklet::WorkletVisitCellsWithPoints
{
using ControlSignature = void(CellSetIn,
FieldInPoint scalars,
FieldOutCell pointsOrder,
FieldOutCell numberOfPoints,
FieldOutCell numberOfTris);
using ExecutionSignature = void(_2 scalar, _3 pointsOrder, _4 numberOfPoints, _5 numberOfTris);
using InputDomain = _1;
VTKM_EXEC FieldType ZCrossProduct(const vtkm::Pair<FieldType, FieldType>& point1,
const vtkm::Pair<FieldType, FieldType>& point2,
const vtkm::Pair<FieldType, FieldType>& point3) const
{
return vtkm::DifferenceOfProducts(point2.first - point1.first,
point3.second - point1.second,
point2.second - point1.second,
point3.first - point1.first);
}
VTKM_EXEC bool DifferentSign(const FieldType& value1, const FieldType& value2) const
{
return (vtkm::IsNegative(value1) != vtkm::IsNegative(value2));
}
VTKM_EXEC bool AllSameSign(const FieldType& value1,
const FieldType& value2,
const FieldType& value3) const
{
return (!(DifferentSign(value1, value2) || DifferentSign(value2, value3)));
}
template <typename ScalarField, typename PointsOutOrder>
VTKM_EXEC void operator()(const ScalarField& scalar,
PointsOutOrder& pointsOrder,
vtkm::IdComponent& numberOfPoints,
vtkm::IdComponent& numberOfTris) const
{
// To classify our tetras following their projection in the 2D scalar domain,
// We consider them as quads, with their coordinates being their respective scalar values.
// To identify the projection, we want to know if the polygon formed by the 4 points of the quad is convex.
// For this, we compute the Z component of the cross product of the vectors of the polygon's edges.
vtkm::Vec<FieldType, 4> scalarCrossProduct{ ZCrossProduct(scalar[0], scalar[1], scalar[2]),
ZCrossProduct(scalar[1], scalar[2], scalar[3]),
ZCrossProduct(scalar[2], scalar[3], scalar[0]),
ZCrossProduct(scalar[3], scalar[0], scalar[1]) };
// If every cross product of consecutive edges of the quad is the same sign, it means that it is convex
// In the case 2 of them are negative and 2 positive, the quad is self-intersecting.
// If one or 3 of them are negative, we have found a non-convex quad, projecting 3 triangles.
if ((DifferentSign(scalarCrossProduct[0], scalarCrossProduct[1]) !=
DifferentSign(scalarCrossProduct[2], scalarCrossProduct[3])))
{
numberOfPoints = 4;
numberOfTris = 3;
// Here, one of the 4 points is in the triangle formed by the 3 other.
// We can guess which one it is by analyzing which cross product has a different sign than the other 3.
// Assign this specific point's id to element 3 of our order array.
if (AllSameSign(scalarCrossProduct[1], scalarCrossProduct[2], scalarCrossProduct[3]))
{
pointsOrder[0] = 0;
pointsOrder[1] = 2;
pointsOrder[2] = 3;
pointsOrder[3] = 1;
}
else if (AllSameSign(scalarCrossProduct[0], scalarCrossProduct[2], scalarCrossProduct[3]))
{
pointsOrder[0] = 0;
pointsOrder[1] = 1;
pointsOrder[2] = 3;
pointsOrder[3] = 2;
}
else if (AllSameSign(scalarCrossProduct[0], scalarCrossProduct[1], scalarCrossProduct[3]))
{
pointsOrder[0] = 0;
pointsOrder[1] = 1;
pointsOrder[2] = 2;
pointsOrder[3] = 3;
}
else
{
pointsOrder[0] = 1;
pointsOrder[1] = 2;
pointsOrder[2] = 3;
pointsOrder[3] = 0;
}
}
else
{
// This tetra projection yields 4 triangles,
// And forms a convex quad in the data plane
numberOfPoints = 5;
numberOfTris = 4;
// Find an order of points which forms a non self-intersecting quad,
// Still using the Z cross-product signs
if (DifferentSign(scalarCrossProduct[0], scalarCrossProduct[1]))
{
// Polygon Self-intersects on the second diagonal
pointsOrder[0] = 0;
pointsOrder[1] = 2;
pointsOrder[2] = 3;
pointsOrder[3] = 1;
}
else if (DifferentSign(scalarCrossProduct[1], scalarCrossProduct[2]))
{
// Self-intersect on the first diagonal
pointsOrder[0] = 0;
pointsOrder[1] = 2;
pointsOrder[2] = 1;
pointsOrder[3] = 3;
}
else
{
// Convex
pointsOrder[0] = 0;
pointsOrder[1] = 1;
pointsOrder[2] = 2;
pointsOrder[3] = 3;
}
}
}
};
template <typename FieldType>
struct VolumeMeasure : vtkm::worklet::WorkletVisitCellsWithPoints
{
using ControlSignature = void(CellSetIn,
FieldInPoint scalars,
FieldInPoint coords,
FieldInCell numberOfTris,
FieldInCell ord,
FieldOutCell newCoords,
FieldOutCell density);
using ExecutionSignature =
void(_2 scalars, _3 coords, _4 numberOfTris, _5 ord, _6 newCoords, _7 density);
using InputDomain = _1;
using Vec3t = vtkm::Vec<FieldType, 3>;
template <typename VecType>
VTKM_EXEC Vec3t DiagonalIntersection(const VecType& point1,
const VecType& point2,
const VecType& point3,
const VecType& point4) const
{
FieldType denominator = vtkm::DifferenceOfProducts(
point1[0] - point2[0], point3[1] - point4[1], point1[1] - point2[1], point3[0] - point4[0]);
// In case the points are aligned, return arbitrarily the first point as intersection.
// These points represent the diagonals of a convex polygon,
// so they are either intersecting or lying on the same line
// The surface area of the quad in the data domain will be null in this case anyways.
if (denominator == 0)
{
return point1;
}
// Otherwise, compute the intersection point of the quad's diagonals, defined by 4 points.
// This vector is the point we get when equating the line equations (point1,point2) and (point3,point4)
return Vec3t{ vtkm::DifferenceOfProducts(
point3[0] - point4[0],
vtkm::DifferenceOfProducts(point1[0], point2[1], point1[1], point2[0]),
point1[0] - point2[0],
vtkm::DifferenceOfProducts(point3[0], point4[1], point3[1], point4[0])) /
denominator,
vtkm::DifferenceOfProducts(
point3[1] - point4[1],
vtkm::DifferenceOfProducts(point1[0], point2[1], point1[1], point2[0]),
point1[1] - point2[1],
vtkm::DifferenceOfProducts(point3[0], point4[1], point3[1], point4[0])) /
denominator,
0.0f };
}
VTKM_EXEC FieldType TriangleArea(Vec3t point1, Vec3t point2, Vec3t point3) const
{
return 0.5f * vtkm::Magnitude(vtkm::Cross(point2 - point1, point3 - point1));
}
VTKM_EXEC FieldType Distance(Vec3t point1, Vec3t point2) const
{
return vtkm::Magnitude(point1 - point2);
}
template <typename ScalarField,
typename CoordsType,
typename NewCoordsType,
typename PointsOutOrder,
typename DensityField>
VTKM_EXEC void operator()(const ScalarField& scalar,
const CoordsType& coords,
const vtkm::IdComponent& numberOfTris,
PointsOutOrder& ord,
NewCoordsType& newCoords,
DensityField& density) const
{
// Write points coordinates in the 2D plane based on provided scalar values
for (int i = 0; i < 4; i++)
{
newCoords[i] = Vec3t{ scalar[i].first, scalar[i].second, 0.0f };
}
if (numberOfTris == 4)
{
// If the projection has 4 triangles, the fifth imaginary point is defined as the intersection of the quad's diagonals.
newCoords[4] = DiagonalIntersection(
newCoords[ord[0]], newCoords[ord[2]], newCoords[ord[1]], newCoords[ord[3]]);
}
// Compute densities
// The density on the borders of the data domain is always null.
// For each tetra projection the only density > 0 is associated either to the point
// located inside the triangle formed by the others (for 3 triangles projection),
// Or to the imaginary point calculated above.
density[ord[0]] = 0.0f;
density[ord[1]] = 0.0f;
density[ord[2]] = 0.0f;
density[ord[3]] = 0.0f;
// Pre-compute some vectors to reuse later
vtkm::Vec<Vec3t, 3> spatialVector = { coords[1] - coords[0],
coords[2] - coords[0],
coords[3] - coords[0] };
vtkm::Vec<Vec3t, 3> spatialCrossProducts = { vtkm::Cross(spatialVector[1], spatialVector[0]),
vtkm::Cross(spatialVector[0], spatialVector[2]),
vtkm::Cross(spatialVector[2],
spatialVector[1]) };
vtkm::Vec<Vec3t, 2> scalarArray = { Vec3t{ scalar[1].first - scalar[0].first,
scalar[2].first - scalar[0].first,
scalar[3].first - scalar[0].first },
Vec3t{ Vec3t{ scalar[1].second - scalar[0].second,
scalar[2].second - scalar[0].second,
scalar[3].second - scalar[0].second } } };
// We need to calculate the determinant in the spatial domain, using the triple product formula
FieldType determinant = vtkm::Dot(spatialVector[2], spatialCrossProducts[0]);
// Calculate the spatial gradient for both scalar fields in the tetrahedron
vtkm::Vec<Vec3t, 2> scalar_gradient;
// The determinant is null, therefore the volume is null
if (determinant == 0.0f)
{
scalar_gradient = { Vec3t{ 0.0f, 0.0f, 0.0f }, Vec3t{ 0.0f, 0.0f, 0.0f } };
}
else
{
// This gradient formulation is derived from the calculation of the inverse Jacobian matrix,
// dividing the adjugate matrix of the transformation by the determinant.
// Each column of the matrix is then multiplied by the scalar difference
// with the scalar value for point with index 0 and summed to get the gradient, for each scalar field.
FieldType inv_determinant = 1.0f / determinant;
scalar_gradient = vtkm::Vec<Vec3t, 2>{ inv_determinant *
(spatialCrossProducts[0] * scalarArray[0][2] +
spatialCrossProducts[1] * scalarArray[0][1] +
spatialCrossProducts[2] * scalarArray[0][0]),
inv_determinant *
(spatialCrossProducts[0] * scalarArray[1][2] +
spatialCrossProducts[1] * scalarArray[1][1] +
spatialCrossProducts[2] * scalarArray[1][0]) };
}
// We get the volume measure, defined as the magnitude of the cross product of the gradient
// See formula (10) in the "continuous scatterplots" paper for the demonstration
FieldType volume = vtkm::Magnitude(vtkm::Cross(scalar_gradient[0], scalar_gradient[1]));
if (numberOfTris == 3)
{
// Calculate the area of the triangle on the backface of the projected tetra
FieldType fullArea =
this->TriangleArea(newCoords[ord[0]], newCoords[ord[1]], newCoords[ord[2]]);
if (volume == 0.0f || fullArea == 0.0f)
{
// For a tetrahedra of null volume, the density is infinite
density[ord[3]] = vtkm::Infinity<FieldType>();
return;
}
// The density for the central point is the distance to the backface divided by the volume
// We interpolate the position of point [3] using the scalar values of the other points
Vec3t contribs{
this->TriangleArea(newCoords[ord[1]], newCoords[ord[2]], newCoords[ord[3]]) / fullArea,
this->TriangleArea(newCoords[ord[0]], newCoords[ord[2]], newCoords[ord[3]]) / fullArea,
this->TriangleArea(newCoords[ord[0]], newCoords[ord[1]], newCoords[ord[3]]) / fullArea
};
Vec3t backfaceProjectionInterpolated = contribs[0] * coords[ord[0]] +
contribs[1] * coords[ord[1]] + contribs[2] * coords[ord[2]];
density[ord[3]] = this->Distance(coords[ord[3]], backfaceProjectionInterpolated) / volume;
}
else
{
// 4 triangles projection
FieldType distanceToIntersection1 = this->Distance(newCoords[4], newCoords[ord[0]]);
FieldType diagonalLength1 = this->Distance(newCoords[ord[2]], newCoords[ord[0]]);
FieldType distanceToIntersection2 = this->Distance(newCoords[4], newCoords[ord[1]]);
FieldType diagonalLength2 = this->Distance(newCoords[ord[1]], newCoords[ord[3]]);
// Spatial volume is null, or data domain surface is null
if (volume == 0.0f || diagonalLength1 == 0.0f || diagonalLength2 == 0.0f)
{
density[4] = vtkm::Infinity<FieldType>();
return;
}
// Interpolate the intersection of diagonals to get its scalar values
FieldType interpolateRatio1 = distanceToIntersection1 / diagonalLength1;
FieldType interpolateRatio2 = distanceToIntersection2 / diagonalLength2;
Vec3t interpolatedPos1 =
coords[ord[0]] + (coords[ord[2]] - coords[ord[0]]) * interpolateRatio1;
Vec3t interpolatedPos2 =
coords[ord[1]] + (coords[ord[3]] - coords[ord[1]]) * interpolateRatio2;
// Eventually, the density is calculated by dividing the scalar mass density by the volume of the cell
density[4] = this->Distance(interpolatedPos1, interpolatedPos2) / volume;
}
}
};
struct ComputeTriangles : vtkm::worklet::WorkletVisitCellsWithPoints
{
using ControlSignature = void(CellSetIn,
FieldInCell pointsOrder,
FieldInCell numberOfTris,
FieldInCell offsets,
FieldOutCell connectivity);
using ExecutionSignature =
void(_2 pointsOrder, _3 numberOfTris, _4 offsets, _5 connectivity, VisitIndex, InputIndex);
using InputDomain = _1;
using ScatterType = vtkm::worklet::ScatterCounting;
template <typename OrderType, typename ConnectivityType>
VTKM_EXEC void operator()(const OrderType& order,
const vtkm::IdComponent& numberOfTris,
const vtkm::Id& offsets,
ConnectivityType& connectivity,
const vtkm::Id& visitIndex,
const vtkm::Id& cellId) const
{
vtkm::Id secondPoint, thirdPoint;
if (numberOfTris == 3)
{
secondPoint = order[(visitIndex + 1) % 3]; // 1,2,0
// The one point in the triangle formed by the 3 other is the common vertex to all 3 visible triangles,
// And has index 3 in the order array
thirdPoint = order[3];
}
else
{
secondPoint = order[(visitIndex + 1) % 4]; // 1,2,3,0
thirdPoint =
4; // Imaginary point at the intersection of diagonals, connected to every triangle
}
connectivity[0] = cellId + offsets + order[static_cast<vtkm::IdComponent>(visitIndex)];
connectivity[1] = cellId + offsets + secondPoint;
connectivity[2] = cellId + offsets + thirdPoint;
}
};
template <typename CoordsComType,
typename CoordsComTypeOut,
typename CoordsInStorageType,
typename OutputCellSetType,
typename CoordsOutStorageType,
typename FieldType>
void Run(const vtkm::cont::CellSetSingleType<>& inputCellSet,
const vtkm::cont::ArrayHandle<vtkm::Vec<CoordsComType, 3>, CoordsInStorageType>& coords,
vtkm::cont::ArrayHandle<vtkm::Vec<CoordsComTypeOut, 3>, CoordsOutStorageType>& newCoords,
vtkm::cont::ArrayHandle<FieldType>& density,
const vtkm::cont::ArrayHandle<FieldType>& field1,
const vtkm::cont::ArrayHandle<FieldType>& field2,
OutputCellSetType& outputCellset)
{
vtkm::cont::Invoker invoke;
// Use zip to pass both scalar fields to worklets as a single argument
auto scalars = vtkm::cont::make_ArrayHandleZip(field1, field2);
// We want to project every tetrahedron in the 2-dimensional data domain using its scalar values,
// following the tetra projection algorithm
// (see "A polygonal approximation to direct scalar volume rendering" by Shirley and Tuchman)
// Minus degenerate cases, this projection makes 3 or 4 triangles in the 2D plane
// This first worklet generates the number of points and triangles needed to project a tetrahedron,
// And the order in which to take them to build the cells
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::IdComponent, 4>> pointsOrder;
vtkm::cont::ArrayHandle<vtkm::IdComponent> numberOfPoints;
vtkm::cont::ArrayHandle<vtkm::IdComponent> numberOfTris;
invoke(
ClassifyTetra<FieldType>{}, inputCellSet, scalars, pointsOrder, numberOfPoints, numberOfTris);
vtkm::Id totalPoints;
vtkm::cont::ArrayHandle<vtkm::Id> offsets =
vtkm::cont::ConvertNumComponentsToOffsets(numberOfPoints, totalPoints);
// Now, compute the tetra's coordinates in the data plane,
// and the density of each projected point
vtkm::cont::ArrayHandle<FieldType> volumeMeasure;
newCoords.Allocate(totalPoints);
density.Allocate(totalPoints);
invoke(VolumeMeasure<FieldType>{},
inputCellSet,
scalars,
coords,
numberOfTris,
pointsOrder,
vtkm::cont::make_ArrayHandleGroupVecVariable(newCoords, offsets),
vtkm::cont::make_ArrayHandleGroupVecVariable(density, offsets));
// Finally, write triangle connectivity in the data domain
vtkm::worklet::ScatterCounting scatter(numberOfTris, true);
vtkm::cont::ArrayHandle<vtkm::Id> offsets_connectivity = scatter.GetInputToOutputMap();
vtkm::cont::ArrayHandle<vtkm::Id> outConnectivity;
invoke(ComputeTriangles{},
scatter,
inputCellSet,
pointsOrder,
numberOfTris,
offsets_connectivity,
vtkm::cont::make_ArrayHandleGroupVec<3>(outConnectivity));
// Create the new dataset
outputCellset.Fill(totalPoints, vtkm::CellShapeTagTriangle::Id, 3, outConnectivity);
}
};
} // namespace worklet
} // namespace vtkm
#endif // vtk_m_worklet_ContinuousScatterPlot_h

11
vtkm/filter/flow/FilterParticleAdvection.h
View File

@ -91,6 +91,16 @@ public:
VTKM_CONT
void SetUseThreadedAlgorithm(bool val) { this->UseThreadedAlgorithm = val; }
VTKM_CONT
void SetUseAsynchronousCommunication() { this->UseAsynchronousCommunication = true; }
VTKM_CONT
bool GetUseAsynchronousCommunication() { return this->UseAsynchronousCommunication; }
VTKM_CONT
void SetUseSynchronousCommunication() { this->UseAsynchronousCommunication = false; }
VTKM_CONT
bool GetUseSynchronousCommunication() { return !this->GetUseAsynchronousCommunication(); }
protected:
VTKM_CONT virtual void ValidateOptions() const;
@ -103,6 +113,7 @@ protected:
vtkm::filter::flow::IntegrationSolverType SolverType =
vtkm::filter::flow::IntegrationSolverType::RK4_TYPE;
vtkm::FloatDefault StepSize = 0;
bool UseAsynchronousCommunication = true;
bool UseThreadedAlgorithm = false;
vtkm::filter::flow::VectorFieldType VecFieldType =
vtkm::filter::flow::VectorFieldType::VELOCITY_FIELD_TYPE;

8
vtkm/filter/flow/FilterParticleAdvectionSteadyState.cxx
View File

@ -92,8 +92,12 @@ VTKM_CONT vtkm::cont::PartitionedDataSet FilterParticleAdvectionSteadyState::DoE
auto dsi = CreateDataSetIntegrators(
input, variant, boundsMap, this->SolverType, this->VecFieldType, this->GetResultType());
vtkm::filter::flow::internal::ParticleAdvector<DSIType> pav(
boundsMap, dsi, this->UseThreadedAlgorithm, this->GetResultType());
vtkm::filter::flow::internal::ParticleAdvector<DSIType> pav(boundsMap,
dsi,
this->UseThreadedAlgorithm,
this->UseAsynchronousCommunication,
this->GetResultType());
return pav.Execute(this->NumberOfSteps, this->StepSize, this->Seeds);
}

7
vtkm/filter/flow/FilterParticleAdvectionUnsteadyState.cxx
View File

@ -75,8 +75,11 @@ VTKM_CONT vtkm::cont::PartitionedDataSet FilterParticleAdvectionUnsteadyState::D
this->VecFieldType,
this->GetResultType());
vtkm::filter::flow::internal::ParticleAdvector<DSIType> pav(
boundsMap, dsi, this->UseThreadedAlgorithm, this->GetResultType());
vtkm::filter::flow::internal::ParticleAdvector<DSIType> pav(boundsMap,
dsi,
this->UseThreadedAlgorithm,
this->UseAsynchronousCommunication,
this->GetResultType());
return pav.Execute(this->NumberOfSteps, this->StepSize, this->Seeds);
}

52
vtkm/filter/flow/internal/AdvectAlgorithm.h
View File

@ -29,11 +29,14 @@ template <typename DSIType, template <typename> class ResultType, typename Parti
class AdvectAlgorithm
{
public:
AdvectAlgorithm(const vtkm::filter::flow::internal::BoundsMap& bm, std::vector<DSIType>& blocks)
AdvectAlgorithm(const vtkm::filter::flow::internal::BoundsMap& bm,
std::vector<DSIType>& blocks,
bool useAsyncComm)
: Blocks(blocks)
, BoundsMap(bm)
, NumRanks(this->Comm.size())
, Rank(this->Comm.rank())
, UseAsynchronousCommunication(useAsyncComm)
{
}
@ -96,7 +99,7 @@ public:
virtual void Go()
{
vtkm::filter::flow::internal::ParticleMessenger<ParticleType> messenger(
this->Comm, this->BoundsMap, 1, 128);
this->Comm, this->UseAsynchronousCommunication, this->BoundsMap, 1, 128);
vtkm::Id nLocal = static_cast<vtkm::Id>(this->Active.size() + this->Inactive.size());
this->ComputeTotalNumParticles(nLocal);
@ -205,6 +208,11 @@ public:
std::vector<ParticleType> incoming;
std::unordered_map<vtkm::Id, std::vector<vtkm::Id>> incomingBlockIDs;
numTermMessages = 0;
bool block = false;
#ifdef VTKM_ENABLE_MPI
block = this->GetBlockAndWait(messenger.UsingSyncCommunication(), numLocalTerminations);
#endif
messenger.Exchange(outgoing,
outgoingRanks,
this->ParticleBlockIDsMap,
@ -212,7 +220,7 @@ public:
incoming,
incomingBlockIDs,
numTermMessages,
this->GetBlockAndWait(numLocalTerminations));
block);
//Cleanup what was sent.
for (const auto& p : outgoing)
@ -320,20 +328,29 @@ public:
return numTerm;
}
virtual bool GetBlockAndWait(const vtkm::Id& numLocalTerm)
virtual bool GetBlockAndWait(const bool& syncComm, const vtkm::Id& numLocalTerm)
{
//There are only two cases where blocking would deadlock.
//1. There are active particles.
//2. numLocalTerm + this->TotalNumberOfTerminatedParticles == this->TotalNumberOfParticles
//So, if neither are true, we can safely block and wait for communication to come in.
bool haveNoWork = this->Active.empty() && this->Inactive.empty();
if (this->Active.empty() && this->Inactive.empty() &&
(numLocalTerm + this->TotalNumTerminatedParticles < this->TotalNumParticles))
//Using syncronous communication we should only block and wait if we have no particles
if (syncComm)
{
return true;
return haveNoWork;
}
else
{
//Otherwise, for asyncronous communication, there are only two cases where blocking would deadlock.
//1. There are active particles.
//2. numLocalTerm + this->TotalNumberOfTerminatedParticles == this->TotalNumberOfParticles
//So, if neither are true, we can safely block and wait for communication to come in.
return false;
if (haveNoWork &&
(numLocalTerm + this->TotalNumTerminatedParticles < this->TotalNumParticles))
return true;
return false;
}
}
//Member data
@ -350,16 +367,7 @@ public:
vtkm::FloatDefault StepSize;
vtkm::Id TotalNumParticles = 0;
vtkm::Id TotalNumTerminatedParticles = 0;
//Active particles:
// {blockID : {particles}
// std::unordered_map<vtkm::Id, std::vector<ParticleType>> Active2;
// {ParticleID : {blockIds}
// std::unordered_map<vtkm::Id, std::vector<vtkm::Id>> ActiveBlockIds;
// std::vector<ParticleType> Inactive2;
// {ParticleID : {blockIds}
// std::unordered_map<vtkm::Id, std::vector<vtkm::Id>> InactiveBlockIds;
bool UseAsynchronousCommunication = true;
};
}

25
vtkm/filter/flow/internal/AdvectAlgorithmThreaded.h
View File

@ -33,8 +33,9 @@ class AdvectAlgorithmThreaded : public AdvectAlgorithm<DSIType, ResultType, Part
{
public:
AdvectAlgorithmThreaded(const vtkm::filter::flow::internal::BoundsMap& bm,
std::vector<DSIType>& blocks)
: AdvectAlgorithm<DSIType, ResultType, ParticleType>(bm, blocks)
std::vector<DSIType>& blocks,
bool useAsyncComm)
: AdvectAlgorithm<DSIType, ResultType, ParticleType>(bm, blocks, useAsyncComm)
, Done(false)
, WorkerActivate(false)
{
@ -133,8 +134,16 @@ protected:
void Manage()
{
if (!this->UseAsynchronousCommunication)
{
VTKM_LOG_S(vtkm::cont::LogLevel::Info,
"Synchronous communication not supported for AdvectAlgorithmThreaded. Forcing "
"asynchronous communication.");
}
bool useAsync = true;
vtkm::filter::flow::internal::ParticleMessenger<ParticleType> messenger(
this->Comm, this->BoundsMap, 1, 128);
this->Comm, useAsync, this->BoundsMap, 1, 128);
while (this->TotalNumTerminatedParticles < this->TotalNumParticles)
{
@ -163,13 +172,15 @@ protected:
this->SetDone();
}
bool GetBlockAndWait(const vtkm::Id& numLocalTerm) override
bool GetBlockAndWait(const bool& syncComm, const vtkm::Id& numLocalTerm) override
{
std::lock_guard<std::mutex> lock(this->Mutex);
if (this->Done)
return true;
return (
this->AdvectAlgorithm<DSIType, ResultType, ParticleType>::GetBlockAndWait(numLocalTerm) &&
!this->WorkerActivate && this->WorkerResults.empty());
return (this->AdvectAlgorithm<DSIType, ResultType, ParticleType>::GetBlockAndWait(
syncComm, numLocalTerm) &&
!this->WorkerActivate && this->WorkerResults.empty());
}
void GetWorkerResults(std::unordered_map<vtkm::Id, std::vector<DSIHelperInfoType>>& results)

136
vtkm/filter/flow/internal/Messenger.cxx
View File

@ -8,6 +8,9 @@
// PURPOSE. See the above copyright notice for more information.
//============================================================================
int R = 1;
#include <vtkm/Math.h>
#include <vtkm/cont/ErrorFilterExecution.h>
#include <vtkm/filter/flow/internal/Messenger.h>
@ -30,13 +33,14 @@ namespace internal
VTKM_CONT
#ifdef VTKM_ENABLE_MPI
Messenger::Messenger(vtkmdiy::mpi::communicator& comm)
Messenger::Messenger(vtkmdiy::mpi::communicator& comm, bool useAsyncComm)
: MPIComm(vtkmdiy::mpi::mpi_cast(comm.handle()))
, MsgID(0)
, NumRanks(comm.size())
, Rank(comm.rank())
, UseAsynchronousCommunication(useAsyncComm)
#else
Messenger::Messenger(vtkmdiy::mpi::communicator& vtkmNotUsed(comm))
Messenger::Messenger(vtkmdiy::mpi::communicator& vtkmNotUsed(comm), bool vtkmNotUsed(useAsyncComm))
#endif
{
}
@ -244,7 +248,7 @@ void Messenger::PrepareForSend(int tag,
}
}
void Messenger::SendData(int dst, int tag, const vtkmdiy::MemoryBuffer& buff)
void Messenger::SendDataAsync(int dst, int tag, const vtkmdiy::MemoryBuffer& buff)
{
std::vector<char*> bufferList;
@ -266,25 +270,25 @@ void Messenger::SendData(int dst, int tag, const vtkmdiy::MemoryBuffer& buff)
}
}
bool Messenger::RecvData(int tag, std::vector<vtkmdiy::MemoryBuffer>& buffers, bool blockAndWait)
void Messenger::SendDataSync(int dst, int tag, vtkmdiy::MemoryBuffer& buff)
{
std::set<int> setTag;
setTag.insert(tag);
std::vector<std::pair<int, vtkmdiy::MemoryBuffer>> b;
buffers.resize(0);
if (this->RecvData(setTag, b, blockAndWait))
auto entry = std::make_pair(dst, std::move(buff));
auto it = this->SyncSendBuffers.find(tag);
if (it == this->SyncSendBuffers.end())
{
buffers.resize(b.size());
for (std::size_t i = 0; i < b.size(); i++)
buffers[i] = std::move(b[i].second);
return true;
std::vector<std::pair<int, vtkmdiy::MemoryBuffer>> vec;
vec.push_back(std::move(entry));
this->SyncSendBuffers.insert(std::make_pair(tag, std::move(vec)));
}
return false;
else
it->second.emplace_back(std::move(entry));
it = this->SyncSendBuffers.find(tag);
}
bool Messenger::RecvData(const std::set<int>& tags,
std::vector<std::pair<int, vtkmdiy::MemoryBuffer>>& buffers,
bool blockAndWait)
bool Messenger::RecvDataAsync(const std::set<int>& tags,
std::vector<std::pair<int, vtkmdiy::MemoryBuffer>>& buffers,
bool blockAndWait)
{
buffers.resize(0);
@ -316,6 +320,104 @@ bool Messenger::RecvData(const std::set<int>& tags,
return !buffers.empty();
}
bool Messenger::RecvDataSync(const std::set<int>& tags,
std::vector<std::pair<int, vtkmdiy::MemoryBuffer>>& buffers,
bool blockAndWait)
{
buffers.resize(0);
if (!blockAndWait)
return false;
//Exchange data
for (auto tag : tags)
{
//Determine the number of messages being sent to each rank and the maximum size.
std::vector<int> maxBuffSize(this->NumRanks, 0), numMessages(this->NumRanks, 0);
const auto& it = this->SyncSendBuffers.find(tag);
if (it != this->SyncSendBuffers.end())
{
for (const auto& i : it->second)
{
int buffSz = i.second.size();
maxBuffSize[i.first] =
vtkm::Max(maxBuffSize[i.first], buffSz); //static_cast<int>(i.second.size()));
numMessages[i.first]++;
}
}
int err = MPI_Allreduce(
MPI_IN_PLACE, maxBuffSize.data(), this->NumRanks, MPI_INT, MPI_MAX, this->MPIComm);
if (err != MPI_SUCCESS)
throw vtkm::cont::ErrorFilterExecution("Error in MPI_Isend inside Messenger::RecvDataSync");
err = MPI_Allreduce(
MPI_IN_PLACE, numMessages.data(), this->NumRanks, MPI_INT, MPI_SUM, this->MPIComm);
if (err != MPI_SUCCESS)
throw vtkm::cont::ErrorFilterExecution("Error in MPI_Isend inside Messenger::RecvDataSync");
MPI_Status status;
std::vector<char> recvBuff;
for (int r = 0; r < this->NumRanks; r++)
{
int numMsgs = numMessages[r];
if (numMsgs == 0)
continue;
//Rank r needs some stuff..
if (this->Rank == r)
{
int maxSz = maxBuffSize[r];
recvBuff.resize(maxSz);
for (int n = 0; n < numMsgs; n++)
{
err =
MPI_Recv(recvBuff.data(), maxSz, MPI_BYTE, MPI_ANY_SOURCE, 0, this->MPIComm, &status);
if (err != MPI_SUCCESS)
throw vtkm::cont::ErrorFilterExecution(
"Error in MPI_Recv inside Messenger::RecvDataSync");
std::pair<int, vtkmdiy::MemoryBuffer> entry;
entry.first = tag;
entry.second.buffer = recvBuff;
buffers.emplace_back(std::move(entry));
}
}
else
{
if (it != this->SyncSendBuffers.end())
{
//it = <tag, <dst,buffer>>
for (const auto& dst_buff : it->second)
{
if (dst_buff.first == r)
{
err = MPI_Send(dst_buff.second.buffer.data(),
dst_buff.second.size(),
MPI_BYTE,
r,
0,
this->MPIComm);
if (err != MPI_SUCCESS)
throw vtkm::cont::ErrorFilterExecution(
"Error in MPI_Send inside Messenger::RecvDataSync");
}
}
}
}
}
//Clean up the send buffers.
if (it != this->SyncSendBuffers.end())
this->SyncSendBuffers.erase(it);
}
MPI_Barrier(this->MPIComm);
return !buffers.empty();
}
void Messenger::ProcessReceivedBuffers(std::vector<char*>& incomingBuffers,
std::vector<std::pair<int, vtkmdiy::MemoryBuffer>>& buffers)
{

50
vtkm/filter/flow/internal/Messenger.h
View File

@ -36,7 +36,7 @@ namespace internal
class VTKM_FILTER_FLOW_EXPORT Messenger
{
public:
VTKM_CONT Messenger(vtkmdiy::mpi::communicator& comm);
VTKM_CONT Messenger(vtkmdiy::mpi::communicator& comm, bool useAsyncComm);
VTKM_CONT virtual ~Messenger()
{
#ifdef VTKM_ENABLE_MPI
@ -50,18 +50,41 @@ public:
#ifdef VTKM_ENABLE_MPI
VTKM_CONT void RegisterTag(int tag, std::size_t numRecvs, std::size_t size);
bool UsingSyncCommunication() const { return !this->UsingAsyncCommunication(); }
bool UsingAsyncCommunication() const { return this->UseAsynchronousCommunication; }
protected:
static std::size_t CalcMessageBufferSize(std::size_t msgSz);
void InitializeBuffers();
void CheckPendingSendRequests();
void CleanupRequests(int tag = TAG_ANY);
void SendData(int dst, int tag, const vtkmdiy::MemoryBuffer& buff);
void SendData(int dst, int tag, vtkmdiy::MemoryBuffer& buff)
{
if (this->UseAsynchronousCommunication)
this->SendDataAsync(dst, tag, buff);
else
this->SendDataSync(dst, tag, buff);
}
bool RecvData(const std::set<int>& tags,
std::vector<std::pair<int, vtkmdiy::MemoryBuffer>>& buffers,
bool blockAndWait = false);
bool blockAndWait = false)
{
if (this->UseAsynchronousCommunication)
return this->RecvDataAsync(tags, buffers, blockAndWait);
else
return this->RecvDataSync(tags, buffers, blockAndWait);
}
private:
void SendDataAsync(int dst, int tag, const vtkmdiy::MemoryBuffer& buff);
void SendDataSync(int dst, int tag, vtkmdiy::MemoryBuffer& buff);
bool RecvDataAsync(const std::set<int>& tags,
std::vector<std::pair<int, vtkmdiy::MemoryBuffer>>& buffers,
bool blockAndWait);
bool RecvDataSync(const std::set<int>& tags,
std::vector<std::pair<int, vtkmdiy::MemoryBuffer>>& buffers,
bool blockAndWait);
void PostRecv(int tag);
void PostRecv(int tag, std::size_t sz, int src = -1);
@ -73,8 +96,6 @@ private:
std::size_t id, numPackets, packet, packetSz, dataSz;
} Header;
bool RecvData(int tag, std::vector<vtkmdiy::MemoryBuffer>& buffers, bool blockAndWait = false);
void PrepareForSend(int tag, const vtkmdiy::MemoryBuffer& buff, std::vector<char*>& buffList);
vtkm::Id GetMsgID() { return this->MsgID++; }
static bool PacketCompare(const char* a, const char* b);
@ -86,6 +107,8 @@ private:
using RankIdPair = std::pair<int, int>;
//Member data
// <tag, {dst, buffer}>
std::map<int, std::vector<std::pair<int, vtkmdiy::MemoryBuffer>>> SyncSendBuffers;
std::map<int, std::pair<std::size_t, std::size_t>> MessageTagInfo;
MPI_Comm MPIComm;
std::size_t MsgID;
@ -95,11 +118,13 @@ private:
std::map<RankIdPair, std::list<char*>> RecvPackets;
std::map<RequestTagPair, char*> SendBuffers;
static constexpr int TAG_ANY = -1;
bool UseAsynchronousCommunication = true;
void CheckRequests(const std::map<RequestTagPair, char*>& buffer,
const std::set<int>& tags,
bool BlockAndWait,
std::vector<RequestTagPair>& reqTags);
#else
protected:
static constexpr int NumRanks = 1;
@ -107,6 +132,21 @@ protected:
#endif
};
template <typename T>
std::ostream& operator<<(std::ostream& os, const std::vector<T>& v)
{
os << "[";
for (std::size_t i = 0; i < v.size(); ++i)
{
os << v[i];
if (i != v.size() - 1)
os << ", ";
}
os << "]";
return os;
}
}
}
}

5
vtkm/filter/flow/internal/ParticleAdvector.h
View File

@ -32,10 +32,12 @@ public:
ParticleAdvector(const vtkm::filter::flow::internal::BoundsMap& bm,
const std::vector<DSIType>& blocks,
const bool& useThreaded,
const bool& useAsyncComm,
const vtkm::filter::flow::FlowResultType& parType)
: Blocks(blocks)
, BoundsMap(bm)
, ResultType(parType)
, UseAsynchronousCommunication(useAsyncComm)
, UseThreadedAlgorithm(useThreaded)
{
}
@ -61,7 +63,7 @@ private:
vtkm::FloatDefault stepSize,
const vtkm::cont::ArrayHandle<ParticleType>& seeds)
{
AlgorithmType algo(this->BoundsMap, this->Blocks);
AlgorithmType algo(this->BoundsMap, this->Blocks, this->UseAsynchronousCommunication);
algo.Execute(numSteps, stepSize, seeds);
return algo.GetOutput();
}
@ -113,6 +115,7 @@ private:
std::vector<DSIType> Blocks;
vtkm::filter::flow::internal::BoundsMap BoundsMap;
FlowResultType ResultType;
bool UseAsynchronousCommunication = true;
bool UseThreadedAlgorithm;
};

4
vtkm/filter/flow/internal/ParticleMessenger.h
View File

@ -44,6 +44,7 @@ class VTKM_FILTER_FLOW_EXPORT ParticleMessenger : public vtkm::filter::flow::int
public:
VTKM_CONT ParticleMessenger(vtkmdiy::mpi::communicator& comm,
bool useAsyncComm,
const vtkm::filter::flow::internal::BoundsMap& bm,
int msgSz = 1,
int numParticles = 128,
@ -113,11 +114,12 @@ VTKM_CONT
template <typename ParticleType>
ParticleMessenger<ParticleType>::ParticleMessenger(
vtkmdiy::mpi::communicator& comm,
bool useAsyncComm,
const vtkm::filter::flow::internal::BoundsMap& boundsMap,
int msgSz,
int numParticles,
int numBlockIds)
: Messenger(comm)
: Messenger(comm, useAsyncComm)
#ifdef VTKM_ENABLE_MPI
, BoundsMap(boundsMap)
#endif

8
vtkm/filter/flow/testing/UnitTestParticleMessengerMPI.cxx
View File

@ -27,11 +27,12 @@ class TestMessenger : public vtkm::filter::flow::internal::ParticleMessenger<vtk
{
public:
TestMessenger(vtkmdiy::mpi::communicator& comm,
bool useAsyncComm,
const vtkm::filter::flow::internal::BoundsMap& bm,
int msgSz = 1,
int numParticles = 1,
int numBlockIds = 1)
: ParticleMessenger(comm, bm, msgSz, numParticles, numBlockIds)
: ParticleMessenger(comm, useAsyncComm, bm, msgSz, numParticles, numBlockIds)
{
}
@ -127,7 +128,8 @@ void TestParticleMessenger()
int maxMsgSz = 100;
int maxNumParticles = 128;
int maxNumBlockIds = 5 * comm.size();
TestMessenger messenger(comm, boundsMap, maxMsgSz / 2, maxNumParticles / 2, maxNumBlockIds / 2);
TestMessenger messenger(
comm, true, boundsMap, maxMsgSz / 2, maxNumParticles / 2, maxNumBlockIds / 2);
//create some data.
std::vector<std::vector<vtkm::Particle>> particles(comm.size());
@ -254,7 +256,7 @@ void TestBufferSizes()
for (const auto& numP : numPs)
for (const auto& nBids : numBids)
{
TestMessenger messenger(comm, boundsMap);
TestMessenger messenger(comm, true, boundsMap);
std::size_t pSize, mSize;
messenger.GetBufferSizes(numP, nBids, mSz, pSize, mSize);

114
vtkm/filter/flow/testing/UnitTestStreamlineFilterMPI.cxx
View File

@ -106,6 +106,7 @@ void SetFilter(FilterType& filter,
const std::string& fieldName,
vtkm::cont::ArrayHandle<vtkm::Particle> seedArray,
bool useThreaded,
bool useAsyncComm,
bool useBlockIds,
const std::vector<vtkm::Id>& blockIds)
{
@ -114,30 +115,42 @@ void SetFilter(FilterType& filter,
filter.SetSeeds(seedArray);
filter.SetActiveField(fieldName);
filter.SetUseThreadedAlgorithm(useThreaded);
// if (useAsyncComm)
// filter.SetUseAsynchronousCommunication();
// else
// filter.SetUseSynchronousCommunication();
if (useBlockIds)
filter.SetBlockIDs(blockIds);
}
void TestAMRStreamline(FilterType fType, bool useThreaded)
void TestAMRStreamline(FilterType fType, bool useThreaded, bool useAsyncComm)
{
switch (fType)
{
case PARTICLE_ADVECTION:
std::cout << "Particle advection";
break;
case STREAMLINE:
std::cout << "Streamline";
break;
case PATHLINE:
std::cout << "Pathline";
break;
}
if (useThreaded)
std::cout << " - using threaded";
std::cout << " - on an AMR data set" << std::endl;
auto comm = vtkm::cont::EnvironmentTracker::GetCommunicator();
if (comm.rank() == 0)
{
switch (fType)
{
case PARTICLE_ADVECTION:
std::cout << "Particle advection";
break;
case STREAMLINE:
std::cout << "Streamline";
break;
case PATHLINE:
std::cout << "Pathline";
break;
}
if (useThreaded)
std::cout << " - using threaded";
if (useAsyncComm)
std::cout << " - usingAsyncComm";
else
std::cout << " - usingSyncComm";
std::cout << " - on an AMR data set" << std::endl;
}
if (comm.size() < 2)
return;
@ -211,13 +224,22 @@ void TestAMRStreamline(FilterType fType, bool useThreaded)
if (fType == STREAMLINE)
{
vtkm::filter::flow::Streamline streamline;
SetFilter(streamline, stepSize, numSteps, fieldName, seedArray, useThreaded, false, {});
SetFilter(streamline,
stepSize,
numSteps,
fieldName,
seedArray,
useThreaded,
useAsyncComm,
false,
{});
out = streamline.Execute(pds);
}
else if (fType == PATHLINE)
{
vtkm::filter::flow::Pathline pathline;
SetFilter(pathline, stepSize, numSteps, fieldName, seedArray, useThreaded, false, {});
SetFilter(
pathline, stepSize, numSteps, fieldName, seedArray, useThreaded, useAsyncComm, false, {});
//Create timestep 2
auto pds2 = vtkm::cont::PartitionedDataSet(pds);
pathline.SetPreviousTime(0);
@ -386,6 +408,7 @@ void TestPartitionedDataSet(vtkm::Id nPerRank,
bool useGhost,
FilterType fType,
bool useThreaded,
bool useAsyncComm,
bool useBlockIds,
bool duplicateBlocks)
{
@ -409,6 +432,11 @@ void TestPartitionedDataSet(vtkm::Id nPerRank,
std::cout << " - using ghost cells";
if (useThreaded)
std::cout << " - using threaded";
if (useAsyncComm)
std::cout << " - usingAsyncComm";
else
std::cout << " - usingSyncComm";
if (useBlockIds)
std::cout << " - using block IDs";
if (duplicateBlocks)
@ -461,8 +489,15 @@ void TestPartitionedDataSet(vtkm::Id nPerRank,
if (fType == STREAMLINE)
{
vtkm::filter::flow::Streamline streamline;
SetFilter(
streamline, stepSize, numSteps, fieldName, seedArray, useThreaded, useBlockIds, blockIds);
SetFilter(streamline,
stepSize,
numSteps,
fieldName,
seedArray,
useThreaded,
useAsyncComm,
useBlockIds,
blockIds);
auto out = streamline.Execute(pds);
vtkm::Id numOutputs = out.GetNumberOfPartitions();
@ -486,6 +521,7 @@ void TestPartitionedDataSet(vtkm::Id nPerRank,
fieldName,
seedArray,
useThreaded,
useAsyncComm,
useBlockIds,
blockIds);
@ -514,8 +550,15 @@ void TestPartitionedDataSet(vtkm::Id nPerRank,
AddVectorFields(pds2, fieldName, vecX);
vtkm::filter::flow::Pathline pathline;
SetFilter(
pathline, stepSize, numSteps, fieldName, seedArray, useThreaded, useBlockIds, blockIds);
SetFilter(pathline,
stepSize,
numSteps,
fieldName,
seedArray,
useThreaded,
useAsyncComm,
useBlockIds,
blockIds);
pathline.SetPreviousTime(time0);
pathline.SetNextTime(time1);
@ -548,21 +591,26 @@ void TestStreamlineFiltersMPI()
for (auto useGhost : flags)
for (auto fType : filterTypes)
for (auto useThreaded : flags)
for (auto useBlockIds : flags)
{
//Run blockIds with and without block duplication.
if (useBlockIds && comm.size() > 1)
for (auto useAsyncComm : flags)
for (auto useBlockIds : flags)
{
TestPartitionedDataSet(n, useGhost, fType, useThreaded, useBlockIds, false);
TestPartitionedDataSet(n, useGhost, fType, useThreaded, useBlockIds, true);
//Run blockIds with and without block duplication.
if (useBlockIds && comm.size() > 1)
{
TestPartitionedDataSet(
n, useGhost, fType, useThreaded, useAsyncComm, useBlockIds, false);
TestPartitionedDataSet(
n, useGhost, fType, useThreaded, useAsyncComm, useBlockIds, true);
}
else
TestPartitionedDataSet(
n, useGhost, fType, useThreaded, useAsyncComm, useBlockIds, false);
}
else
TestPartitionedDataSet(n, useGhost, fType, useThreaded, useBlockIds, false);
}
for (auto fType : filterTypes)
for (auto useThreaded : flags)
TestAMRStreamline(fType, useThreaded);
for (auto useAsyncComm : flags)
TestAMRStreamline(fType, useThreaded, useAsyncComm);
}
}

5
vtkm/worklet/ScatterPermutation.h
View File

@ -59,7 +59,10 @@ public:
VTKM_CONT OutputToInputMapType GetOutputToInputMap() const { return this->Permutation; }
VTKM_CONT
VisitArrayType GetVisitArray(vtkm::Id inputRange) const { return VisitArrayType(0, inputRange); }
VisitArrayType GetVisitArray(vtkm::Id inputRange) const
{
return VisitArrayType(0, this->GetOutputRange(inputRange));
}
VTKM_CONT
VisitArrayType GetVisitArray(vtkm::Id3 inputRange) const

2
vtkm/worklet/internal/DispatcherBase.h
View File

@ -759,7 +759,7 @@ private:
// Get the arrays used for masking output elements.
typename MaskType::ThreadToOutputMapType threadToOutputMap =
this->Mask.GetThreadToOutputMap(inputRange);
this->Mask.GetThreadToOutputMap(outputRange);
// Replace the parameters in the invocation with the execution object and
// pass to next step of Invoke. Also add the scatter information.