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

Browse Source
This commit is contained in:
Dave Pugmire 2023-04-18 15:41:31 -04:00
commit d79b2cb5be
34 changed files with 2295 additions and 312 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

13
vtkm/filter/flow/FilterParticleAdvection.cxx
View File

@ -10,9 +10,12 @@
#include <vtkm/Particle.h>
#include <vtkm/cont/EnvironmentTracker.h>
#include <vtkm/cont/ErrorFilterExecution.h>
#include <vtkm/filter/flow/FilterParticleAdvection.h>
#include <vtkm/thirdparty/diy/diy.h>
namespace vtkm
{
namespace filter
@ -33,7 +36,15 @@ VTKM_CONT void FilterParticleAdvection::ValidateOptions() const
{
if (this->GetUseCoordinateSystemAsField())
throw vtkm::cont::ErrorFilterExecution("Coordinate system as field not supported");
if (this->Seeds.GetNumberOfValues() == 0)
vtkm::Id numSeeds = this->Seeds.GetNumberOfValues();
#ifdef VTKM_ENABLE_MPI
vtkmdiy::mpi::communicator comm = vtkm::cont::EnvironmentTracker::GetCommunicator();
vtkm::Id totalNumSeeds = 0;
vtkmdiy::mpi::all_reduce(comm, numSeeds, totalNumSeeds, std::plus<vtkm::Id>{});
numSeeds = totalNumSeeds;
#endif
if (numSeeds == 0)
throw vtkm::cont::ErrorFilterExecution("No seeds provided.");
if (!this->Seeds.IsBaseComponentType<vtkm::Particle>() &&
!this->Seeds.IsBaseComponentType<vtkm::ChargedParticle>())

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

@ -53,6 +53,12 @@ public:
this->Seeds = vtkm::cont::make_ArrayHandle(seeds, copyFlag);
}
VTKM_CONT void SetBlockIDs(const std::vector<vtkm::Id>& blockIds)
{
this->BlockIds = blockIds;
this->BlockIdsSet = true;
}
VTKM_CONT
void SetSolverRK4() { this->SolverType = vtkm::filter::flow::IntegrationSolverType::RK4_TYPE; }
@ -85,16 +91,29 @@ 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;
VTKM_CONT virtual vtkm::filter::flow::FlowResultType GetResultType() const = 0;
bool BlockIdsSet = false;
std::vector<vtkm::Id> BlockIds;
vtkm::Id NumberOfSteps = 0;
vtkm::cont::UnknownArrayHandle Seeds;
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;

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

@ -83,12 +83,21 @@ VTKM_CONT vtkm::cont::PartitionedDataSet FilterParticleAdvectionSteadyState::DoE
variant.Emplace<DSIType::ElectroMagneticFieldNameType>(electric, magnetic);
}
vtkm::filter::flow::internal::BoundsMap boundsMap(input);
vtkm::filter::flow::internal::BoundsMap boundsMap;
if (this->BlockIdsSet)
boundsMap = vtkm::filter::flow::internal::BoundsMap(input, this->BlockIds);
else
boundsMap = vtkm::filter::flow::internal::BoundsMap(input);
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);
}

145
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)
{
}
@ -77,13 +80,17 @@ public:
const ParticleType p = portal.Get(i);
std::vector<vtkm::Id> ids = this->BoundsMap.FindBlocks(p.GetPosition());
if (!ids.empty() && this->BoundsMap.FindRank(ids[0]) == this->Rank)
//Note: For duplicate blocks, this will give the seeds to the rank that are first in the list.
if (!ids.empty())
{
particles.emplace_back(p);
blockIDs.emplace_back(ids);
auto ranks = this->BoundsMap.FindRank(ids[0]);
if (!ranks.empty() && this->Rank == ranks[0])
{
particles.emplace_back(p);
blockIDs.emplace_back(ids);
}
}
}
this->SetSeedArray(particles, blockIDs);
}
@ -91,7 +98,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);
@ -119,7 +126,6 @@ public:
}
}
virtual void ClearParticles()
{
this->Active.clear();
@ -160,11 +166,13 @@ public:
bit++;
}
//Update Active
this->Active.insert(this->Active.end(), particles.begin(), particles.end());
}
virtual bool GetActiveParticles(std::vector<ParticleType>& particles, vtkm::Id& blockId)
{
//DRP: particles = std::move(this->Active2[blockId]);
particles.clear();
blockId = -1;
if (this->Active.empty())
@ -187,23 +195,97 @@ public:
return !particles.empty();
}
virtual void Communicate(vtkm::filter::flow::internal::ParticleMessenger<ParticleType>& messenger,
vtkm::Id numLocalTerminations,
vtkm::Id& numTermMessages)
void Communicate(vtkm::filter::flow::internal::ParticleMessenger<ParticleType>& messenger,
vtkm::Id numLocalTerminations,
vtkm::Id& numTermMessages)
{
std::vector<ParticleType> outgoing;
std::vector<vtkm::Id> outgoingRanks;
this->GetOutgoingParticles(outgoing, outgoingRanks);
std::vector<ParticleType> incoming;
std::unordered_map<vtkm::Id, std::vector<vtkm::Id>> incomingIDs;
std::unordered_map<vtkm::Id, std::vector<vtkm::Id>> incomingBlockIDs;
numTermMessages = 0;
messenger.Exchange(this->Inactive,
bool block = false;
#ifdef VTKM_ENABLE_MPI
block = this->GetBlockAndWait(messenger.UsingSyncCommunication(), numLocalTerminations);
#endif
messenger.Exchange(outgoing,
outgoingRanks,
this->ParticleBlockIDsMap,
numLocalTerminations,
incoming,
incomingIDs,
incomingBlockIDs,
numTermMessages,
this->GetBlockAndWait(numLocalTerminations));
block);
//Cleanup what was sent.
for (const auto& p : outgoing)
this->ParticleBlockIDsMap.erase(p.GetID());
this->UpdateActive(incoming, incomingBlockIDs);
}
void GetOutgoingParticles(std::vector<ParticleType>& outgoing,
std::vector<vtkm::Id>& outgoingRanks)
{
outgoing.clear();
outgoingRanks.clear();
outgoing.reserve(this->Inactive.size());
outgoingRanks.reserve(this->Inactive.size());
std::vector<ParticleType> particlesStaying;
std::unordered_map<vtkm::Id, std::vector<vtkm::Id>> particlesStayingBlockIDs;
//Send out Everything.
for (const auto& p : this->Inactive)
{
const auto& bid = this->ParticleBlockIDsMap[p.GetID()];
VTKM_ASSERT(!bid.empty());
auto ranks = this->BoundsMap.FindRank(bid[0]);
VTKM_ASSERT(!ranks.empty());
if (ranks.size() == 1)
{
if (ranks[0] == this->Rank)
{
particlesStaying.emplace_back(p);
particlesStayingBlockIDs[p.GetID()] = this->ParticleBlockIDsMap[p.GetID()];
}
else
{
outgoing.emplace_back(p);
outgoingRanks.emplace_back(ranks[0]);
}
}
else
{
//Decide where it should go...
//Random selection:
vtkm::Id outRank = std::rand() % ranks.size();
if (outRank == this->Rank)
{
particlesStayingBlockIDs[p.GetID()] = this->ParticleBlockIDsMap[p.GetID()];
particlesStaying.emplace_back(p);
}
else
{
outgoing.emplace_back(p);
outgoingRanks.emplace_back(outRank);
}
}
}
this->Inactive.clear();
this->UpdateActive(incoming, incomingIDs);
VTKM_ASSERT(outgoing.size() == outgoingRanks.size());
VTKM_ASSERT(particlesStaying.size() == particlesStayingBlockIDs.size());
if (!particlesStaying.empty())
this->UpdateActive(particlesStaying, particlesStayingBlockIDs);
}
virtual void UpdateActive(const std::vector<ParticleType>& particles,
@ -231,8 +313,8 @@ public:
vtkm::Id UpdateResult(const DSIHelperInfo<ParticleType>& stuff)
{
this->UpdateActive(stuff.A, stuff.IdMapA);
this->UpdateInactive(stuff.I, stuff.IdMapI);
this->UpdateActive(stuff.InBounds.Particles, stuff.InBounds.BlockIDs);
this->UpdateInactive(stuff.OutOfBounds.Particles, stuff.OutOfBounds.BlockIDs);
vtkm::Id numTerm = static_cast<vtkm::Id>(stuff.TermID.size());
//Update terminated particles.
@ -245,20 +327,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
@ -269,11 +360,13 @@ public:
std::vector<ParticleType> Inactive;
vtkm::Id NumberOfSteps;
vtkm::Id NumRanks;
//{particleId : {block IDs}}
std::unordered_map<vtkm::Id, std::vector<vtkm::Id>> ParticleBlockIDsMap;
vtkm::Id Rank;
vtkm::FloatDefault StepSize;
vtkm::Id TotalNumParticles = 0;
vtkm::Id TotalNumTerminatedParticles = 0;
bool UseAsynchronousCommunication = true;
};
}

28
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,11 +134,20 @@ 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)
{
std::unordered_map<vtkm::Id, std::vector<DSIHelperInfoType>> workerResults;
this->GetWorkerResults(workerResults);
@ -145,7 +155,9 @@ protected:
for (auto& it : workerResults)
{
for (auto& r : it.second)
{
numTerm += this->UpdateResult(r.Get<DSIHelperInfo<ParticleType>>());
}
}
vtkm::Id numTermMessages = 0;
@ -160,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)

181
vtkm/filter/flow/internal/BoundsMap.h
View File

@ -15,6 +15,7 @@
#include <vtkm/cont/AssignerPartitionedDataSet.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/EnvironmentTracker.h>
#include <vtkm/cont/ErrorFilterExecution.h>
#include <vtkm/cont/Field.h>
#include <vtkm/cont/PartitionedDataSet.h>
@ -25,6 +26,9 @@
#include <vtkm/thirdparty/diy/mpi-cast.h>
#endif
#include <algorithm>
#include <iostream>
#include <set>
#include <vector>
namespace vtkm
@ -39,31 +43,22 @@ namespace internal
class VTKM_ALWAYS_EXPORT BoundsMap
{
public:
BoundsMap()
: LocalNumBlocks(0)
, TotalNumBlocks(0)
BoundsMap() {}
BoundsMap(const vtkm::cont::DataSet& dataSet) { this->Init({ dataSet }); }
BoundsMap(const vtkm::cont::DataSet& dataSet, const vtkm::Id& blockId)
{
this->Init({ dataSet }, { blockId });
}
BoundsMap(const vtkm::cont::DataSet& dataSet)
: LocalNumBlocks(1)
, TotalNumBlocks(0)
{
this->Init({ dataSet });
}
BoundsMap(const std::vector<vtkm::cont::DataSet>& dataSets) { this->Init(dataSets); }
BoundsMap(const std::vector<vtkm::cont::DataSet>& dataSets)
: LocalNumBlocks(static_cast<vtkm::Id>(dataSets.size()))
, TotalNumBlocks(0)
{
this->Init(dataSets);
}
BoundsMap(const vtkm::cont::PartitionedDataSet& pds) { this->Init(pds.GetPartitions()); }
BoundsMap(const vtkm::cont::PartitionedDataSet& pds)
: LocalNumBlocks(pds.GetNumberOfPartitions())
, TotalNumBlocks(0)
BoundsMap(const vtkm::cont::PartitionedDataSet& pds, const std::vector<vtkm::Id>& blockIds)
{
this->Init(pds.GetPartitions());
this->Init(pds.GetPartitions(), blockIds);
}
vtkm::Id GetLocalBlockId(vtkm::Id idx) const
@ -72,11 +67,12 @@ public:
return this->LocalIDs[static_cast<std::size_t>(idx)];
}
int FindRank(vtkm::Id blockId) const
std::vector<int> FindRank(vtkm::Id blockId) const
{
auto it = this->BlockToRankMap.find(blockId);
if (it == this->BlockToRankMap.end())
return -1;
return {};
return it->second;
}
@ -110,8 +106,94 @@ public:
vtkm::Id GetLocalNumBlocks() const { return this->LocalNumBlocks; }
private:
void Init(const std::vector<vtkm::cont::DataSet>& dataSets, const std::vector<vtkm::Id>& blockIds)
{
if (dataSets.size() != blockIds.size())
throw vtkm::cont::ErrorFilterExecution("Number of datasets and block ids must match");
this->LocalIDs = blockIds;
this->LocalNumBlocks = dataSets.size();
vtkmdiy::mpi::communicator comm = vtkm::cont::EnvironmentTracker::GetCommunicator();
//1. Get the min/max blockId
vtkm::Id locMinId = 0, locMaxId = 0;
if (!this->LocalIDs.empty())
{
locMinId = *std::min_element(this->LocalIDs.begin(), this->LocalIDs.end());
locMaxId = *std::max_element(this->LocalIDs.begin(), this->LocalIDs.end());
}
vtkm::Id globalMinId = 0, globalMaxId = 0;
vtkmdiy::mpi::all_reduce(comm, locMinId, globalMinId, vtkmdiy::mpi::minimum<vtkm::Id>{});
vtkmdiy::mpi::all_reduce(comm, locMaxId, globalMaxId, vtkmdiy::mpi::maximum<vtkm::Id>{});
if (globalMinId != 0 || (globalMaxId - globalMinId) < 1)
throw vtkm::cont::ErrorFilterExecution("Invalid block ids");
//2. Find out how many blocks everyone has.
std::vector<vtkm::Id> locBlockCounts(comm.size(), 0), globalBlockCounts(comm.size(), 0);
locBlockCounts[comm.rank()] = this->LocalIDs.size();
vtkmdiy::mpi::all_reduce(comm, locBlockCounts, globalBlockCounts, std::plus<vtkm::Id>{});
//note: there might be duplicates...
vtkm::Id globalNumBlocks =
std::accumulate(globalBlockCounts.begin(), globalBlockCounts.end(), 0);
//3. given the counts per rank, calc offset for this rank.
vtkm::Id offset = 0;
for (int i = 0; i < comm.rank(); i++)
offset += globalBlockCounts[i];
//4. calc the blocks on each rank.
std::vector<vtkm::Id> localBlockIds(globalNumBlocks, 0);
vtkm::Id idx = 0;
for (const auto& bid : this->LocalIDs)
localBlockIds[offset + idx++] = bid;
//use an MPI_Alltoallv instead.
std::vector<vtkm::Id> globalBlockIds(globalNumBlocks, 0);
vtkmdiy::mpi::all_reduce(comm, localBlockIds, globalBlockIds, std::plus<vtkm::Id>{});
//5. create a rank -> blockId map.
// rankToBlockIds[rank] = {this->LocalIDs on rank}.
std::vector<std::vector<vtkm::Id>> rankToBlockIds(comm.size());
offset = 0;
for (int rank = 0; rank < comm.size(); rank++)
{
vtkm::Id numBIds = globalBlockCounts[rank];
rankToBlockIds[rank].resize(numBIds);
for (vtkm::Id i = 0; i < numBIds; i++)
rankToBlockIds[rank][i] = globalBlockIds[offset + i];
offset += numBIds;
}
//6. there might be duplicates, so count number of UNIQUE blocks.
std::set<vtkm::Id> globalUniqueBlockIds;
globalUniqueBlockIds.insert(globalBlockIds.begin(), globalBlockIds.end());
this->TotalNumBlocks = globalUniqueBlockIds.size();
//Build a vector of : blockIdsToRank[blockId] = {ranks that have this blockId}
std::vector<std::vector<vtkm::Id>> blockIdsToRank(this->TotalNumBlocks);
for (int rank = 0; rank < comm.size(); rank++)
{
for (const auto& bid : rankToBlockIds[rank])
{
blockIdsToRank[bid].push_back(rank);
this->BlockToRankMap[bid].push_back(rank);
}
}
this->Build(dataSets);
}
void Init(const std::vector<vtkm::cont::DataSet>& dataSets)
{
this->LocalNumBlocks = dataSets.size();
vtkm::cont::AssignerPartitionedDataSet assigner(this->LocalNumBlocks);
this->TotalNumBlocks = assigner.nblocks();
std::vector<int> ids;
@ -122,7 +204,7 @@ private:
this->LocalIDs.emplace_back(static_cast<vtkm::Id>(i));
for (vtkm::Id id = 0; id < this->TotalNumBlocks; id++)
this->BlockToRankMap[id] = assigner.rank(static_cast<int>(id));
this->BlockToRankMap[id] = { assigner.rank(static_cast<int>(id)) };
this->Build(dataSets);
}
@ -131,48 +213,61 @@ private:
std::vector<vtkm::Float64> vals(static_cast<std::size_t>(this->TotalNumBlocks * 6), 0);
std::vector<vtkm::Float64> vals2(vals.size());
std::vector<vtkm::Float64> localMins((this->TotalNumBlocks * 3),
std::numeric_limits<vtkm::Float64>::max());
std::vector<vtkm::Float64> localMaxs((this->TotalNumBlocks * 3),
-std::numeric_limits<vtkm::Float64>::max());
for (std::size_t i = 0; i < this->LocalIDs.size(); i++)
{
const vtkm::cont::DataSet& ds = dataSets[static_cast<std::size_t>(i)];
vtkm::Bounds bounds = ds.GetCoordinateSystem().GetBounds();
std::size_t idx = static_cast<std::size_t>(this->LocalIDs[i] * 6);
vals[idx++] = bounds.X.Min;
vals[idx++] = bounds.X.Max;
vals[idx++] = bounds.Y.Min;
vals[idx++] = bounds.Y.Max;
vals[idx++] = bounds.Z.Min;
vals[idx++] = bounds.Z.Max;
vtkm::Id localID = this->LocalIDs[i];
localMins[localID * 3 + 0] = bounds.X.Min;
localMins[localID * 3 + 1] = bounds.Y.Min;
localMins[localID * 3 + 2] = bounds.Z.Min;
localMaxs[localID * 3 + 0] = bounds.X.Max;
localMaxs[localID * 3 + 1] = bounds.Y.Max;
localMaxs[localID * 3 + 2] = bounds.Z.Max;
}
std::vector<vtkm::Float64> globalMins, globalMaxs;
#ifdef VTKM_ENABLE_MPI
vtkmdiy::mpi::communicator Comm = vtkm::cont::EnvironmentTracker::GetCommunicator();
MPI_Comm mpiComm = vtkmdiy::mpi::mpi_cast(Comm.handle());
MPI_Allreduce(vals.data(), vals2.data(), vals.size(), MPI_DOUBLE, MPI_SUM, mpiComm);
globalMins.resize(this->TotalNumBlocks * 3);
globalMaxs.resize(this->TotalNumBlocks * 3);
vtkmdiy::mpi::communicator comm = vtkm::cont::EnvironmentTracker::GetCommunicator();
vtkmdiy::mpi::all_reduce(comm, localMins, globalMins, vtkmdiy::mpi::minimum<vtkm::Float64>{});
vtkmdiy::mpi::all_reduce(comm, localMaxs, globalMaxs, vtkmdiy::mpi::maximum<vtkm::Float64>{});
#else
vals2 = vals;
globalMins = localMins;
globalMaxs = localMaxs;
#endif
this->BlockBounds.resize(static_cast<std::size_t>(this->TotalNumBlocks));
this->GlobalBounds = vtkm::Bounds();
std::size_t idx = 0;
for (auto& block : this->BlockBounds)
{
block = vtkm::Bounds(vals2[idx + 0],
vals2[idx + 1],
vals2[idx + 2],
vals2[idx + 3],
vals2[idx + 4],
vals2[idx + 5]);
block = vtkm::Bounds(globalMins[idx + 0],
globalMaxs[idx + 0],
globalMins[idx + 1],
globalMaxs[idx + 1],
globalMins[idx + 2],
globalMaxs[idx + 2]);
this->GlobalBounds.Include(block);
idx += 6;
idx += 3;
}
}
vtkm::Id LocalNumBlocks;
vtkm::Id LocalNumBlocks = 0;
std::vector<vtkm::Id> LocalIDs;
std::map<vtkm::Id, vtkm::Int32> BlockToRankMap;
vtkm::Id TotalNumBlocks;
std::map<vtkm::Id, std::vector<vtkm::Int32>> BlockToRankMap;
vtkm::Id TotalNumBlocks = 0;
std::vector<vtkm::Bounds> BlockBounds;
vtkm::Bounds GlobalBounds;
};

110
vtkm/filter/flow/internal/DataSetIntegrator.h
View File

@ -42,16 +42,63 @@ public:
const std::unordered_map<vtkm::Id, std::vector<vtkm::Id>>& particleBlockIDsMap)
: BoundsMap(boundsMap)
, ParticleBlockIDsMap(particleBlockIDsMap)
, V(v)
, Particles(v)
{
}
struct ParticleBlockIds
{
void Clear()
{
this->Particles.clear();
this->BlockIDs.clear();
}
void Add(const ParticleType& p, const std::vector<vtkm::Id>& bids)
{
this->Particles.emplace_back(p);
this->BlockIDs[p.GetID()] = std::move(bids);
}
std::vector<ParticleType> Particles;
std::unordered_map<vtkm::Id, std::vector<vtkm::Id>> BlockIDs;
};
void Clear()
{
this->InBounds.Clear();
this->OutOfBounds.Clear();
this->TermIdx.clear();
this->TermID.clear();
}
void Validate(vtkm::Id num)
{
#ifndef NDEBUG
//Make sure we didn't miss anything. Every particle goes into a single bucket.
VTKM_ASSERT(static_cast<std::size_t>(num) ==
(this->InBounds.Particles.size() + this->OutOfBounds.Particles.size() +
this->TermIdx.size()));
VTKM_ASSERT(this->InBounds.Particles.size() == this->InBounds.BlockIDs.size());
VTKM_ASSERT(this->OutOfBounds.Particles.size() == this->OutOfBounds.BlockIDs.size());
VTKM_ASSERT(this->TermIdx.size() == this->TermID.size());
#endif
}
void AddTerminated(vtkm::Id idx, vtkm::Id pID)
{
this->TermIdx.emplace_back(idx);
this->TermID.emplace_back(pID);
}
const vtkm::filter::flow::internal::BoundsMap BoundsMap;
const std::unordered_map<vtkm::Id, std::vector<vtkm::Id>> ParticleBlockIDsMap;
std::vector<ParticleType> A, I, V;
std::unordered_map<vtkm::Id, std::vector<vtkm::Id>> IdMapA, IdMapI;
std::vector<vtkm::Id> TermIdx, TermID;
ParticleBlockIds InBounds;
ParticleBlockIds OutOfBounds;
std::vector<ParticleType> Particles;
std::vector<vtkm::Id> TermID;
std::vector<vtkm::Id> TermIdx;
};
using DSIHelperInfoType =
@ -147,12 +194,13 @@ VTKM_CONT inline void DataSetIntegrator<Derived>::ClassifyParticles(
const vtkm::cont::ArrayHandle<ParticleType>& particles,
DSIHelperInfo<ParticleType>& dsiInfo) const
{
dsiInfo.A.clear();
dsiInfo.I.clear();
dsiInfo.TermID.clear();
dsiInfo.TermIdx.clear();
dsiInfo.IdMapI.clear();
dsiInfo.IdMapA.clear();
/*
each particle: --> T,I,A
if T: update TermIdx, TermID
if A: update IdMapA;
if I: update IdMapI;
*/
dsiInfo.Clear();
auto portal = particles.WritePortal();
vtkm::Id n = portal.GetNumberOfValues();
@ -161,13 +209,15 @@ VTKM_CONT inline void DataSetIntegrator<Derived>::ClassifyParticles(
{
auto p = portal.Get(i);
//Terminated.
if (p.GetStatus().CheckTerminate())
{
dsiInfo.TermIdx.emplace_back(i);
dsiInfo.TermID.emplace_back(p.GetID());
dsiInfo.AddTerminated(i, p.GetID());
}
else
{
//Didn't terminate.
//Get the blockIDs.
const auto& it = dsiInfo.ParticleBlockIDsMap.find(p.GetID());
VTKM_ASSERT(it != dsiInfo.ParticleBlockIDsMap.end());
auto currBIDs = it->second;
@ -175,9 +225,17 @@ VTKM_CONT inline void DataSetIntegrator<Derived>::ClassifyParticles(
std::vector<vtkm::Id> newIDs;
if (p.GetStatus().CheckSpatialBounds() && !p.GetStatus().CheckTookAnySteps())
{
//particle is OUTSIDE but didn't take any steps.
//this means that the particle wasn't in the block.
//assign newIDs to currBIDs[1:]
newIDs.assign(std::next(currBIDs.begin(), 1), currBIDs.end());
}
else
{
//Otherwise, get new blocks from the current position.
newIDs = dsiInfo.BoundsMap.FindBlocks(p.GetPosition(), currBIDs);
}
//reset the particle status.
p.GetStatus() = vtkm::ParticleStatus();
@ -185,19 +243,19 @@ VTKM_CONT inline void DataSetIntegrator<Derived>::ClassifyParticles(
if (newIDs.empty()) //No blocks, we're done.
{
p.GetStatus().SetTerminate();
dsiInfo.TermIdx.emplace_back(i);
dsiInfo.TermID.emplace_back(p.GetID());
dsiInfo.AddTerminated(i, p.GetID());
}
else
{
//If we have more than blockId, we want to minimize communication
//If we have more than one blockId, we want to minimize communication
//and put any blocks owned by this rank first.
if (newIDs.size() > 1)
{
for (auto idit = newIDs.begin(); idit != newIDs.end(); idit++)
{
vtkm::Id bid = *idit;
if (dsiInfo.BoundsMap.FindRank(bid) == this->Rank)
auto ranks = dsiInfo.BoundsMap.FindRank(bid);
if (std::find(ranks.begin(), ranks.end(), this->Rank) != ranks.end())
{
newIDs.erase(idit);
newIDs.insert(newIDs.begin(), bid);
@ -206,26 +264,14 @@ VTKM_CONT inline void DataSetIntegrator<Derived>::ClassifyParticles(
}
}
int dstRank = dsiInfo.BoundsMap.FindRank(newIDs[0]);
if (dstRank == this->Rank)
{
dsiInfo.A.emplace_back(p);
dsiInfo.IdMapA[p.GetID()] = newIDs;
}
else
{
dsiInfo.I.emplace_back(p);
dsiInfo.IdMapI[p.GetID()] = newIDs;
}
dsiInfo.OutOfBounds.Add(p, newIDs);
}
portal.Set(i, p);
}
portal.Set(i, p);
}
//Make sure we didn't miss anything. Every particle goes into a single bucket.
VTKM_ASSERT(static_cast<std::size_t>(n) ==
(dsiInfo.A.size() + dsiInfo.I.size() + dsiInfo.TermIdx.size()));
VTKM_ASSERT(dsiInfo.TermIdx.size() == dsiInfo.TermID.size());
//Make sure everything is copacetic.
dsiInfo.Validate(n);
}
template <typename Derived>

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

@ -199,7 +199,7 @@ VTKM_CONT inline void DataSetIntegratorSteadyState::DoAdvect(DSIHelperInfo<vtkm:
using ArrayType = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
auto copyFlag = (this->CopySeedArray ? vtkm::CopyFlag::On : vtkm::CopyFlag::Off);
auto seedArray = vtkm::cont::make_ArrayHandle(b.V, copyFlag);
auto seedArray = vtkm::cont::make_ArrayHandle(b.Particles, copyFlag);
if (this->VecFieldType == VectorFieldType::VELOCITY_FIELD_TYPE)
{
@ -238,7 +238,7 @@ VTKM_CONT inline void DataSetIntegratorSteadyState::DoAdvect(
using ArrayType = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
auto copyFlag = (this->CopySeedArray ? vtkm::CopyFlag::On : vtkm::CopyFlag::Off);
auto seedArray = vtkm::cont::make_ArrayHandle(b.V, copyFlag);
auto seedArray = vtkm::cont::make_ArrayHandle(b.Particles, copyFlag);
if (this->VecFieldType == VectorFieldType::ELECTRO_MAGNETIC_FIELD_TYPE)
{

2
vtkm/filter/flow/internal/DataSetIntegratorUnsteadyState.h
View File

@ -200,7 +200,7 @@ VTKM_CONT inline void DataSetIntegratorUnsteadyState::DoAdvect(DSIHelperInfo<vtk
using ArrayType = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
auto copyFlag = (this->CopySeedArray ? vtkm::CopyFlag::On : vtkm::CopyFlag::Off);
auto seedArray = vtkm::cont::make_ArrayHandle(b.V, copyFlag);
auto seedArray = vtkm::cont::make_ArrayHandle(b.Particles, copyFlag);
using AHType = internal::AdvectHelper<internal::UnsteadyStateGridEvalType, vtkm::Particle>;

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;
};

20
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,
@ -51,6 +52,7 @@ public:
VTKM_CONT ~ParticleMessenger() {}
VTKM_CONT void Exchange(const std::vector<ParticleType>& outData,
const std::vector<vtkm::Id>& outRanks,
const std::unordered_map<vtkm::Id, std::vector<vtkm::Id>>& outBlockIDsMap,
vtkm::Id numLocalTerm,
std::vector<ParticleType>& inData,
@ -96,6 +98,7 @@ protected:
VTKM_CONT void SerialExchange(
const std::vector<ParticleType>& outData,
const std::vector<vtkm::Id>& outRanks,
const std::unordered_map<vtkm::Id, std::vector<vtkm::Id>>& outBlockIDsMap,
vtkm::Id numLocalTerm,
std::vector<ParticleType>& inData,
@ -111,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
@ -166,6 +170,7 @@ VTKM_CONT
template <typename ParticleType>
void ParticleMessenger<ParticleType>::SerialExchange(
const std::vector<ParticleType>& outData,
const std::vector<vtkm::Id>& vtkmNotUsed(outRanks),
const std::unordered_map<vtkm::Id, std::vector<vtkm::Id>>& outBlockIDsMap,
vtkm::Id vtkmNotUsed(numLocalTerm),
std::vector<ParticleType>& inData,
@ -184,6 +189,7 @@ VTKM_CONT
template <typename ParticleType>
void ParticleMessenger<ParticleType>::Exchange(
const std::vector<ParticleType>& outData,
const std::vector<vtkm::Id>& outRanks,
const std::unordered_map<vtkm::Id, std::vector<vtkm::Id>>& outBlockIDsMap,
vtkm::Id numLocalTerm,
std::vector<ParticleType>& inData,
@ -191,23 +197,25 @@ void ParticleMessenger<ParticleType>::Exchange(
vtkm::Id& numTerminateMessages,
bool blockAndWait)
{
VTKM_ASSERT(outData.size() == outRanks.size());
numTerminateMessages = 0;
inDataBlockIDsMap.clear();
if (this->GetNumRanks() == 1)
return this->SerialExchange(
outData, outBlockIDsMap, numLocalTerm, inData, inDataBlockIDsMap, blockAndWait);
outData, outRanks, outBlockIDsMap, numLocalTerm, inData, inDataBlockIDsMap, blockAndWait);
#ifdef VTKM_ENABLE_MPI
//dstRank, vector of (particles,blockIDs)
std::unordered_map<int, std::vector<ParticleCommType>> sendData;
for (const auto& p : outData)
std::size_t numP = outData.size();
for (std::size_t i = 0; i < numP; i++)
{
const auto& bids = outBlockIDsMap.find(p.GetID())->second;
int dstRank = this->BoundsMap.FindRank(bids[0]);
sendData[dstRank].emplace_back(std::make_pair(p, bids));
const auto& bids = outBlockIDsMap.find(outData[i].GetID())->second;
sendData[outRanks[i]].emplace_back(std::make_pair(outData[i], bids));
}
//Do all the sends first.

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);

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

@ -45,40 +45,112 @@ void AddVectorFields(vtkm::cont::PartitionedDataSet& pds,
ds.AddPointField(fieldName, CreateConstantVectorField(ds.GetNumberOfPoints(), vec));
}
std::vector<vtkm::cont::PartitionedDataSet> CreateAllDataSetBounds(vtkm::Id nPerRank, bool useGhost)
{
auto comm = vtkm::cont::EnvironmentTracker::GetCommunicator();
vtkm::Id totNumBlocks = nPerRank * comm.size();
vtkm::Id numDims = 5;
vtkm::FloatDefault x0 = 0;
vtkm::FloatDefault x1 = x0 + static_cast<vtkm::FloatDefault>(numDims - 1);
vtkm::FloatDefault dx = x1 - x0;
vtkm::FloatDefault y0 = 0, y1 = numDims - 1, z0 = 0, z1 = numDims - 1;
if (useGhost)
{
numDims = numDims + 2; //add 1 extra on each side
x0 = x0 - 1;
x1 = x1 + 1;
dx = x1 - x0 - 2;
y0 = y0 - 1;
y1 = y1 + 1;
z0 = z0 - 1;
z1 = z1 + 1;
}
//Create ALL of the blocks.
std::vector<vtkm::Bounds> bounds;
for (vtkm::Id i = 0; i < totNumBlocks; i++)
{
bounds.push_back(vtkm::Bounds(x0, x1, y0, y1, z0, z1));
x0 += dx;
x1 += dx;
}
const vtkm::Id3 dims(numDims, numDims, numDims);
auto allPDS = vtkm::worklet::testing::CreateAllDataSets(bounds, dims, useGhost);
return allPDS;
}
std::vector<vtkm::Range> ExtractMaxXRanges(const vtkm::cont::PartitionedDataSet& pds, bool useGhost)
{
std::vector<vtkm::Range> xMaxRanges;
for (const auto& ds : pds.GetPartitions())
{
auto bounds = ds.GetCoordinateSystem().GetBounds();
auto xMax = bounds.X.Max;
if (useGhost)
xMax = xMax - 1;
xMaxRanges.push_back(vtkm::Range(xMax, xMax + static_cast<vtkm::FloatDefault>(.5)));
}
return xMaxRanges;
}
template <typename FilterType>
void SetFilter(FilterType& filter,
vtkm::FloatDefault stepSize,
vtkm::Id numSteps,
const std::string& fieldName,
vtkm::cont::ArrayHandle<vtkm::Particle> seedArray,
bool useThreaded)
bool useThreaded,
bool useAsyncComm,
bool useBlockIds,
const std::vector<vtkm::Id>& blockIds)
{
filter.SetStepSize(stepSize);
filter.SetNumberOfSteps(numSteps);
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;
@ -152,13 +224,22 @@ void TestAMRStreamline(FilterType fType, bool useThreaded)
if (fType == STREAMLINE)
{
vtkm::filter::flow::Streamline streamline;
SetFilter(streamline, stepSize, numSteps, fieldName, seedArray, useThreaded);
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);
SetFilter(
pathline, stepSize, numSteps, fieldName, seedArray, useThreaded, useAsyncComm, false, {});
//Create timestep 2
auto pds2 = vtkm::cont::PartitionedDataSet(pds);
pathline.SetPreviousTime(0);
@ -277,17 +358,21 @@ void TestAMRStreamline(FilterType fType, bool useThreaded)
void ValidateOutput(const vtkm::cont::DataSet& out,
vtkm::Id numSeeds,
vtkm::Range xMaxRange,
FilterType fType)
const vtkm::Range& xMaxRange,
FilterType fType,
bool checkEndPoint,
bool blockDuplication)
{
//Validate the result is correct.
VTKM_TEST_ASSERT(out.GetNumberOfCoordinateSystems() == 1,
"Wrong number of coordinate systems in the output dataset");
vtkm::cont::UnknownCellSet dcells = out.GetCellSet();
out.PrintSummary(std::cout);
std::cout << " nSeeds= " << numSeeds << std::endl;
VTKM_TEST_ASSERT(dcells.GetNumberOfCells() == numSeeds, "Wrong number of cells");
vtkm::Id numCells = out.GetNumberOfCells();
if (!blockDuplication)
VTKM_TEST_ASSERT(numCells == numSeeds, "Wrong number of cells");
auto coords = out.GetCoordinateSystem().GetDataAsMultiplexer();
auto ptPortal = coords.ReadPortal();
@ -296,89 +381,94 @@ void ValidateOutput(const vtkm::cont::DataSet& out,
vtkm::cont::CellSetExplicit<> explicitCells;
VTKM_TEST_ASSERT(dcells.IsType<vtkm::cont::CellSetExplicit<>>(), "Wrong cell type.");
explicitCells = dcells.AsCellSet<vtkm::cont::CellSetExplicit<>>();
for (vtkm::Id j = 0; j < numSeeds; j++)
for (vtkm::Id j = 0; j < numCells; j++)
{
vtkm::cont::ArrayHandle<vtkm::Id> indices;
explicitCells.GetIndices(j, indices);
vtkm::Id nPts = indices.GetNumberOfValues();
auto iPortal = indices.ReadPortal();
vtkm::Vec3f lastPt = ptPortal.Get(iPortal.Get(nPts - 1));
VTKM_TEST_ASSERT(xMaxRange.Contains(lastPt[0]), "Wrong end point for seed");
if (checkEndPoint)
VTKM_TEST_ASSERT(xMaxRange.Contains(lastPt[0]), "Wrong end point for seed");
}
}
else if (fType == PARTICLE_ADVECTION)
{
VTKM_TEST_ASSERT(out.GetNumberOfPoints() == numSeeds, "Wrong number of coordinates");
for (vtkm::Id i = 0; i < numSeeds; i++)
VTKM_TEST_ASSERT(xMaxRange.Contains(ptPortal.Get(i)[0]), "Wrong end point for seed");
if (!blockDuplication)
VTKM_TEST_ASSERT(out.GetNumberOfPoints() == numSeeds, "Wrong number of coordinates");
if (checkEndPoint)
{
for (vtkm::Id i = 0; i < numCells; i++)
VTKM_TEST_ASSERT(xMaxRange.Contains(ptPortal.Get(i)[0]), "Wrong end point for seed");
}
}
}
void TestPartitionedDataSet(vtkm::Id nPerRank, bool useGhost, FilterType fType, bool useThreaded)
void TestPartitionedDataSet(vtkm::Id nPerRank,
bool useGhost,
FilterType fType,
bool useThreaded,
bool useAsyncComm,
bool useBlockIds,
bool duplicateBlocks)
{
switch (fType)
{
case PARTICLE_ADVECTION:
std::cout << "Particle advection";
break;
case STREAMLINE:
std::cout << "Streamline";
break;
case PATHLINE:
std::cout << "Pathline";
break;
}
if (useGhost)
std::cout << " - using ghost cells";
if (useThreaded)
std::cout << " - using threaded";
std::cout << " - on a partitioned data set" << std::endl;
auto comm = vtkm::cont::EnvironmentTracker::GetCommunicator();
vtkm::Id numDims = 5;
vtkm::FloatDefault x0 = static_cast<vtkm::FloatDefault>((numDims - 1) * (nPerRank * comm.rank()));
vtkm::FloatDefault x1 = x0 + static_cast<vtkm::FloatDefault>(numDims - 1);
vtkm::FloatDefault dx = x1 - x0;
vtkm::FloatDefault y0 = 0, y1 = numDims - 1, z0 = 0, z1 = numDims - 1;
if (useGhost)
if (comm.rank() == 0)
{
numDims = numDims + 2; //add 1 extra on each side
x0 = x0 - 1;
x1 = x1 + 1;
dx = x1 - x0 - 2;
y0 = y0 - 1;
y1 = y1 + 1;
z0 = z0 - 1;
z1 = z1 + 1;
}
vtkm::FloatDefault time0 = 0;
vtkm::FloatDefault time1 = (numDims - 1) * (nPerRank * comm.size());
std::vector<vtkm::Bounds> bounds;
for (vtkm::Id i = 0; i < nPerRank; i++)
{
bounds.push_back(vtkm::Bounds(x0, x1, y0, y1, z0, z1));
x0 += dx;
x1 += dx;
}
std::vector<vtkm::Range> xMaxRanges;
for (const auto& b : bounds)
{
auto xMax = b.X.Max;
switch (fType)
{
case PARTICLE_ADVECTION:
std::cout << "Particle advection";
break;
case STREAMLINE:
std::cout << "Streamline";
break;
case PATHLINE:
std::cout << "Pathline";
break;
}
std::cout << " blocksPerRank= " << nPerRank;
if (useGhost)
xMax = xMax - 1;
xMaxRanges.push_back(vtkm::Range(xMax, xMax + static_cast<vtkm::FloatDefault>(.5)));
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)
std::cout << " - with duplicate blocks";
std::cout << " - on a partitioned data set" << std::endl;
}
std::vector<vtkm::Id> blockIds;
//Uniform assignment.
for (vtkm::Id i = 0; i < nPerRank; i++)
blockIds.push_back(comm.rank() * nPerRank + i);
//For block duplication, give everyone the 2nd to last block.
//We want to keep the last block on the last rank for validation.
if (duplicateBlocks && blockIds.size() > 1)
{
vtkm::Id totNumBlocks = comm.size() * nPerRank;
vtkm::Id dupBlock = totNumBlocks - 2;
for (int r = 0; r < comm.size(); r++)
{
if (std::find(blockIds.begin(), blockIds.end(), dupBlock) == blockIds.end())
blockIds.push_back(dupBlock);
}
}
std::vector<vtkm::cont::PartitionedDataSet> allPDS, allPDS2;
const vtkm::Id3 dims(numDims, numDims, numDims);
allPDS = vtkm::worklet::testing::CreateAllDataSets(bounds, dims, useGhost);
allPDS2 = vtkm::worklet::testing::CreateAllDataSets(bounds, dims, useGhost);
allPDS = CreateAllDataSetBounds(nPerRank, useGhost);
allPDS2 = CreateAllDataSetBounds(nPerRank, useGhost);
auto xMaxRanges = ExtractMaxXRanges(allPDS[0], useGhost);
vtkm::FloatDefault time0 = 0;
vtkm::FloatDefault time1 = xMaxRanges[xMaxRanges.size() - 1].Max;
vtkm::Vec3f vecX(1, 0, 0);
std::string fieldName = "vec";
@ -386,7 +476,9 @@ void TestPartitionedDataSet(vtkm::Id nPerRank, bool useGhost, FilterType fType,
vtkm::Id numSteps = 100000;
for (std::size_t n = 0; n < allPDS.size(); n++)
{
auto pds = allPDS[n];
vtkm::cont::PartitionedDataSet pds;
for (const auto& bid : blockIds)
pds.AppendPartition(allPDS[n].GetPartition(bid));
AddVectorFields(pds, fieldName, vecX);
vtkm::cont::ArrayHandle<vtkm::Particle> seedArray;
@ -397,62 +489,129 @@ void TestPartitionedDataSet(vtkm::Id nPerRank, bool useGhost, FilterType fType,
if (fType == STREAMLINE)
{
vtkm::filter::flow::Streamline streamline;
SetFilter(streamline, stepSize, numSteps, fieldName, seedArray, useThreaded);
SetFilter(streamline,
stepSize,
numSteps,
fieldName,
seedArray,
useThreaded,
useAsyncComm,
useBlockIds,
blockIds);
auto out = streamline.Execute(pds);
for (vtkm::Id i = 0; i < nPerRank; i++)
ValidateOutput(out.GetPartition(i), numSeeds, xMaxRanges[i], fType);
vtkm::Id numOutputs = out.GetNumberOfPartitions();
bool checkEnds = numOutputs == static_cast<vtkm::Id>(blockIds.size());
for (vtkm::Id i = 0; i < numOutputs; i++)
{
ValidateOutput(out.GetPartition(i),
numSeeds,
xMaxRanges[blockIds[i]],
fType,
checkEnds,
duplicateBlocks);
}
}
else if (fType == PARTICLE_ADVECTION)
{
vtkm::filter::flow::ParticleAdvection particleAdvection;
SetFilter(particleAdvection, stepSize, numSteps, fieldName, seedArray, useThreaded);
SetFilter(particleAdvection,
stepSize,
numSteps,
fieldName,
seedArray,
useThreaded,
useAsyncComm,
useBlockIds,
blockIds);
auto out = particleAdvection.Execute(pds);
//Particles end up in last rank.
if (comm.rank() == comm.size() - 1)
{
bool checkEnds = out.GetNumberOfPartitions() == static_cast<vtkm::Id>(blockIds.size());
VTKM_TEST_ASSERT(out.GetNumberOfPartitions() == 1, "Wrong number of partitions in output");
ValidateOutput(out.GetPartition(0), numSeeds, xMaxRanges[xMaxRanges.size() - 1], fType);
ValidateOutput(out.GetPartition(0),
numSeeds,
xMaxRanges[xMaxRanges.size() - 1],
fType,
checkEnds,
duplicateBlocks);
}
else
VTKM_TEST_ASSERT(out.GetNumberOfPartitions() == 0, "Wrong number of partitions in output");
}
else if (fType == PATHLINE)
{
auto pds2 = allPDS2[n];
vtkm::cont::PartitionedDataSet pds2;
for (const auto& bid : blockIds)
pds2.AppendPartition(allPDS2[n].GetPartition(bid));
AddVectorFields(pds2, fieldName, vecX);
vtkm::filter::flow::Pathline pathline;
SetFilter(pathline, stepSize, numSteps, fieldName, seedArray, useThreaded);
SetFilter(pathline,
stepSize,
numSteps,
fieldName,
seedArray,
useThreaded,
useAsyncComm,
useBlockIds,
blockIds);
pathline.SetPreviousTime(time0);
pathline.SetNextTime(time1);
pathline.SetNextDataSet(pds2);
auto out = pathline.Execute(pds);
for (vtkm::Id i = 0; i < nPerRank; i++)
ValidateOutput(out.GetPartition(i), numSeeds, xMaxRanges[i], fType);
vtkm::Id numOutputs = out.GetNumberOfPartitions();
bool checkEnds = numOutputs == static_cast<vtkm::Id>(blockIds.size());
for (vtkm::Id i = 0; i < numOutputs; i++)
ValidateOutput(out.GetPartition(i),
numSeeds,
xMaxRanges[blockIds[i]],
fType,
checkEnds,
duplicateBlocks);
}
}
}
void TestStreamlineFiltersMPI()
{
auto comm = vtkm::cont::EnvironmentTracker::GetCommunicator();
if (comm.rank() == 0)
std::cout << std::endl << "*** TestStreamlineFiltersMPI" << std::endl;
std::vector<bool> flags = { true, false };
std::vector<FilterType> filterTypes = { PARTICLE_ADVECTION, STREAMLINE, PATHLINE };
for (int n = 1; n < 3; n++)
{
for (auto useGhost : flags)
for (auto fType : filterTypes)
for (auto useThreaded : flags)
TestPartitionedDataSet(n, useGhost, fType, useThreaded);
}
for (auto useAsyncComm : flags)
for (auto useBlockIds : flags)
{
//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);
}
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.