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

This commit is contained in:
Dave Pugmire 2022-08-05 16:05:37 -04:00
commit f6326f19c0
136 changed files with 4496 additions and 3594 deletions

@ -164,8 +164,9 @@ stages:
interruptible: true
before_script:
- *install_cmake
- .gitlab/ci/config/sccache.sh
- "cmake -VV -P .gitlab/ci/config/ninja.cmake"
- export PATH=$PWD/.gitlab:$PATH
- .gitlab/ci/config/sccache.sh
- SCCACHE_IDLE_TIMEOUT=0 sccache --start-server
- sccache --show-stats
- .gitlab/ci/config/google_benchmarks.sh
@ -185,6 +186,8 @@ stages:
interruptible: true
before_script:
- *install_cmake
- "cmake -VV -P .gitlab/ci/config/ninja.cmake"
- export PATH=$PWD/.gitlab:$PATH
script:
- "ctest $CTEST_TIMEOUT -VV -S .gitlab/ci/ctest_test.cmake"
extends:

@ -0,0 +1,61 @@
##============================================================================
## 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.
##============================================================================
cmake_minimum_required(VERSION 3.0 FATAL_ERROR)
set(version 4.6.1)
set(arch x86_64)
if(CMAKE_HOST_SYSTEM_NAME STREQUAL "Linux")
set(sha256sum da1e1781bc1c4b019216fa16391af3e1daaee7e7f49a8ec9b0cdc8a1d05c50e2)
set(base_url https://github.com/ccache/ccache/releases/download)
set(platform linux)
set(extension tar.xz)
elseif(CMAKE_HOST_SYSTEM_NAME STREQUAL "Darwin")
set(sha256sum 3e36ba8c80fbf7f2b95fe0227b9dd1ca6143d721aab052caf0d5729769138059)
set(full_url https://gitlab.kitware.com/utils/ci-utilities/-/package_files/534/download)
set(filename ccache)
set(extension tar.gz)
elseif(CMAKE_HOST_SYSTEM_NAME STREQUAL "Windows")
set(sha256sum a6c6311973aa3d2aae22424895f2f968e5d661be003b25f1bd854a5c0cd57563)
set(base_url https://github.com/ccache/ccache/releases/download)
set(platform windows)
set(extension zip)
else()
message(FATAL_ERROR "Unrecognized platform ${CMAKE_HOST_SYSTEM_NAME}")
endif()
if(NOT DEFINED filename)
set(filename "ccache-${version}-${platform}-${arch}")
endif()
set(tarball "${filename}.${extension}")
if(NOT DEFINED full_url)
set(full_url "${base_url}/v${version}/${tarball}")
endif()
file(DOWNLOAD
"${full_url}" .gitlab/${tarball}
EXPECTED_HASH SHA256=${sha256sum}
SHOW_PROGRESS
)
execute_process(
COMMAND ${CMAKE_COMMAND} -E tar xf ${tarball}
WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/.gitlab
RESULT_VARIABLE extract_results
)
if(extract_results)
message(FATAL_ERROR "Extracting `${tarball}` failed: ${extract_results}.")
endif()
file(RENAME .gitlab/${filename} .gitlab/ccache)

@ -0,0 +1,44 @@
##============================================================================
## 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.
##============================================================================
cmake_minimum_required(VERSION 3.0 FATAL_ERROR)
set(version 1.11.0)
if(CMAKE_HOST_SYSTEM_NAME STREQUAL "Linux")
set(sha256sum 9726e730d5b8599f82654dc80265e64a10a8a817552c34153361ed0c017f9f02)
set(platform linux)
elseif(CMAKE_HOST_SYSTEM_NAME STREQUAL "Darwin")
set(sha256sum 21915277db59756bfc61f6f281c1f5e3897760b63776fd3d360f77dd7364137f)
set(platform mac)
elseif(CMAKE_HOST_SYSTEM_NAME STREQUAL "Windows")
set(sha256sum d0ee3da143211aa447e750085876c9b9d7bcdd637ab5b2c5b41349c617f22f3b)
set(platform win)
else()
message(FATAL_ERROR "Unrecognized platform ${CMAKE_HOST_SYSTEM_NAME}")
endif()
set(tarball "ninja-${platform}.zip")
file(DOWNLOAD
"https://github.com/ninja-build/ninja/releases/download/v${version}/${tarball}" .gitlab/${tarball}
EXPECTED_HASH SHA256=${sha256sum}
SHOW_PROGRESS
)
execute_process(
COMMAND ${CMAKE_COMMAND} -E tar xf ${tarball}
WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/.gitlab
RESULT_VARIABLE extract_results
)
if(extract_results)
message(FATAL_ERROR "Extracting `${tarball}` failed: ${extract_results}.")
endif()

@ -42,7 +42,13 @@ test:macos_xcode13:
before_script:
- .gitlab/ci/config/cmake.sh
- export PATH=$PWD/.gitlab/cmake/bin:$PATH
- "cmake -VV -P .gitlab/ci/config/ccache.cmake"
- export PATH=$PWD/.gitlab/ccache:$PATH
- "cmake -VV -P .gitlab/ci/config/ninja.cmake"
- export PATH=$PWD/.gitlab:$PATH
- "cmake --version"
- "ccache --version"
- "ninja --version"
- "cmake -V -P .gitlab/ci/config/fetch_vtkm_tags.cmake"
- "cmake -V -P .gitlab/ci/config/gitlab_ci_setup.cmake"
- "ctest -VV -S .gitlab/ci/ctest_configure.cmake"
@ -59,8 +65,11 @@ test:macos_xcode13:
interruptible: true
before_script:
- .gitlab/ci/config/cmake.sh
- export PATH=$PWD/.gitlab/cmake/bin:$PATH
- export PATH=.gitlab/cmake/bin:$PATH
- "cmake -VV -P .gitlab/ci/config/ninja.cmake"
- export PATH=$PWD/.gitlab:$PATH
- cmake --version
- ninja --version
script:
- "ctest $CTEST_TIMEOUT -VV -S .gitlab/ci/ctest_test.cmake"
extends:

@ -29,8 +29,12 @@
- Invoke-Expression -Command .gitlab/ci/config/cmake.ps1
- Invoke-Expression -Command .gitlab/ci/config/vcvarsall.ps1
- $pwdpath = $pwd.Path
- Set-Item -Force -Path "env:PATH" -Value "$pwdpath\.gitlab;$pwdpath\.gitlab\cmake\bin;$env:PATH"
- Set-Item -Force -Path "env:PATH" -Value "$pwdpath\.gitlab\cmake\bin;$env:PATH"
- "cmake --version"
- "cmake -V -P .gitlab/ci/config/ccache.cmake"
- Set-Item -Force -Path "env:PATH" -Value "$pwdpath\.gitlab\ccache;$env:PATH"
- "cmake -V -P .gitlab/ci/config/ninja.cmake"
- Set-Item -Force -Path "env:PATH" -Value "$pwdpath\.gitlab;$env:PATH"
- "cmake -V -P .gitlab/ci/config/gitlab_ci_setup.cmake"
- "ctest -VV -S .gitlab/ci/ctest_configure.cmake"
script:
@ -67,7 +71,10 @@
- Invoke-Expression -Command .gitlab/ci/config/cmake.ps1
- Invoke-Expression -Command .gitlab/ci/config/vcvarsall.ps1
- $pwdpath = $pwd.Path
- Set-Item -Force -Path "env:PATH" -Value "$pwdpath\.gitlab;$pwdpath\.gitlab\cmake\bin;$env:PATH"
- Set-Item -Force -Path "env:PATH" -Value "$pwdpath\.gitlab\cmake\bin;$env:PATH"
- "cmake --version"
- "cmake -V -P .gitlab/ci/config/ninja.cmake"
- Set-Item -Force -Path "env:PATH" -Value "$pwdpath\.gitlab;$env:PATH"
script:
- "ctest -VV -S .gitlab/ci/ctest_test.cmake"

@ -70,32 +70,29 @@
# Use TBBConfig.cmake if possible.
# Disabling this as it running the TBBConfig.cmake on dragnipur is
# causing a CMake error. I don't know if this is an install problem
# or an issue with version 2018.0.
# set(_tbb_find_quiet)
# if (TBB_FIND_QUIETLY)
# set(_tbb_find_quiet QUIET)
# endif ()
# set(_tbb_find_components)
# set(_tbb_find_optional_components)
# foreach (_tbb_find_component IN LISTS TBB_FIND_COMPONENTS)
# if (TBB_FIND_REQUIRED_${_tbb_find_component})
# list(APPEND _tbb_find_components "${_tbb_find_component}")
# else ()
# list(APPEND _tbb_find_optional_components "${_tbb_find_component}")
# endif ()
# endforeach ()
# unset(_tbb_find_component)
# find_package(TBB CONFIG ${_tbb_find_quiet}
# COMPONENTS ${_tbb_find_components}
# OPTIONAL_COMPONENTS ${_tbb_find_optional_components})
# unset(_tbb_find_quiet)
# unset(_tbb_find_components)
# unset(_tbb_find_optional_components)
# if (TBB_FOUND)
# return ()
# endif ()
set(_tbb_find_quiet)
if (TBB_FIND_QUIETLY)
set(_tbb_find_quiet QUIET)
endif ()
set(_tbb_find_components)
set(_tbb_find_optional_components)
foreach (_tbb_find_component IN LISTS TBB_FIND_COMPONENTS)
if (TBB_FIND_REQUIRED_${_tbb_find_component})
list(APPEND _tbb_find_components "${_tbb_find_component}")
else ()
list(APPEND _tbb_find_optional_components "${_tbb_find_component}")
endif ()
endforeach ()
unset(_tbb_find_component)
find_package(TBB CONFIG ${_tbb_find_quiet}
COMPONENTS ${_tbb_find_components}
OPTIONAL_COMPONENTS ${_tbb_find_optional_components})
unset(_tbb_find_quiet)
unset(_tbb_find_components)
unset(_tbb_find_optional_components)
if (TBB_FOUND)
return ()
endif ()
#====================================================
# Fix the library path in case it is a linker script

@ -42,7 +42,11 @@ function(vtkm_extract_real_library library real_library)
endfunction()
if(VTKm_ENABLE_TBB AND NOT TARGET vtkm::tbb)
find_package(TBB REQUIRED)
# Skip find_package(TBB) if we already have it
if (NOT TARGET TBB::tbb)
find_package(TBB REQUIRED)
endif()
add_library(vtkmTBB INTERFACE)
add_library(vtkm::tbb ALIAS vtkmTBB)
target_link_libraries(vtkmTBB INTERFACE TBB::tbb)

@ -14,7 +14,10 @@
#include <sstream>
#include <vtkm/ImplicitFunction.h>
#include <vtkm/Particle.h>
#include <vtkm/cont/Algorithm.h>
#include <vtkm/cont/ArrayCopy.h>
#include <vtkm/cont/BoundsCompute.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/FieldRangeCompute.h>
@ -25,9 +28,9 @@
#include <vtkm/cont/internal/OptionParser.h>
#include <vtkm/filter/Streamline.h>
#include <vtkm/filter/contour/Contour.h>
#include <vtkm/filter/contour/Slice.h>
#include <vtkm/filter/flow/Streamline.h>
#include <vtkm/filter/geometry_refinement/Tetrahedralize.h>
#include <vtkm/filter/geometry_refinement/Tube.h>
#include <vtkm/filter/vector_analysis/Gradient.h>
@ -423,7 +426,7 @@ void AddField(vtkm::cont::PartitionedDataSet& input,
}
template <typename DataSetType>
DataSetType RunStreamlinesHelper(vtkm::filter::Streamline& filter, const DataSetType& input)
DataSetType RunStreamlinesHelper(vtkm::filter::flow::Streamline& filter, const DataSetType& input)
{
auto dataSetBounds = vtkm::cont::BoundsCompute(input);
vtkm::cont::ArrayHandle<vtkm::Particle> seedArray;
@ -457,7 +460,7 @@ void BenchStreamlines(::benchmark::State& state)
BuildInputDataSet(cycle, isStructured, isMultiBlock, DataSetDim);
inputGenTimer.Stop();
vtkm::filter::Streamline streamline;
vtkm::filter::flow::Streamline streamline;
streamline.SetStepSize(0.1f);
streamline.SetNumberOfSteps(1000);
streamline.SetActiveField(PointVectorsName);

@ -10,17 +10,14 @@
#include "Benchmarker.h"
#include <vtkm/Particle.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/DataSetBuilderUniform.h>
#include <vtkm/cont/ErrorInternal.h>
#include <vtkm/cont/Logging.h>
#include <vtkm/cont/RuntimeDeviceTracker.h>
#include <vtkm/cont/Timer.h>
#include <vtkm/cont/internal/OptionParser.h>
#include <vtkm/filter/ParticleAdvection.h>
#include <vtkm/worklet/particleadvection/EulerIntegrator.h>
#include <vtkm/worklet/particleadvection/RK4Integrator.h>
#include <vtkm/filter/flow/ParticleAdvection.h>
namespace
{
@ -50,7 +47,7 @@ void BenchParticleAdvection(::benchmark::State& state)
vtkm::Particle(vtkm::Vec3f(.2f, 2.0f, .2f), 1),
vtkm::Particle(vtkm::Vec3f(.2f, 3.0f, .2f), 2) });
vtkm::filter::ParticleAdvection particleAdvection;
vtkm::filter::flow::ParticleAdvection particleAdvection;
particleAdvection.SetStepSize(vtkm::FloatDefault(1) / state.range(0));
particleAdvection.SetNumberOfSteps(static_cast<vtkm::Id>(state.range(0)));

@ -548,14 +548,8 @@ int main(int argc, char* argv[])
static_cast<vtkm::Id>(dims[1]),
static_cast<vtkm::Id>(0));
vtkm::cont::ArrayHandle<vtkm::Id3> localBlockIndices;
vtkm::cont::ArrayHandle<vtkm::Id3> localBlockOrigins;
vtkm::cont::ArrayHandle<vtkm::Id3> localBlockSizes;
localBlockIndices.Allocate(blocksPerRank);
localBlockOrigins.Allocate(blocksPerRank);
localBlockSizes.Allocate(blocksPerRank);
auto localBlockIndicesPortal = localBlockIndices.WritePortal();
auto localBlockOriginsPortal = localBlockOrigins.WritePortal();
auto localBlockSizesPortal = localBlockSizes.WritePortal();
{
vtkm::Id lastDimSize =
@ -604,30 +598,29 @@ int main(int argc, char* argv[])
vtkm::Vec<ValueType, 2> origin(0, blockIndex * blockSize);
vtkm::Vec<ValueType, 2> spacing(1, 1);
ds = dsb.Create(vdims, origin, spacing);
vtkm::cont::CellSetStructured<2> cs;
cs.SetPointDimensions(vdims);
cs.SetGlobalPointDimensions(vtkm::Id2{ globalSize[0], globalSize[1] });
cs.SetGlobalPointIndexStart(vtkm::Id2{ 0, blockIndex * blockSize });
ds.SetCellSet(cs);
localBlockIndicesPortal.Set(localBlockIndex, vtkm::Id3(blockIndex, 0, 0));
localBlockOriginsPortal.Set(localBlockIndex,
vtkm::Id3((blockStart / blockSliceSize), 0, 0));
localBlockSizesPortal.Set(localBlockIndex,
vtkm::Id3(currBlockSize, static_cast<vtkm::Id>(dims[0]), 0));
localBlockIndicesPortal.Set(localBlockIndex, vtkm::Id3(0, blockIndex, 0));
}
// 3D data
else
{
vtkm::Id3 vdims;
vdims[0] = static_cast<vtkm::Id>(dims[1]);
vdims[1] = static_cast<vtkm::Id>(dims[0]);
vdims[0] = static_cast<vtkm::Id>(dims[0]);
vdims[1] = static_cast<vtkm::Id>(dims[1]);
vdims[2] = static_cast<vtkm::Id>(currBlockSize);
vtkm::Vec<ValueType, 3> origin(0, 0, (blockIndex * blockSize));
vtkm::Vec<ValueType, 3> origin(0, 0, blockIndex * blockSize);
vtkm::Vec<ValueType, 3> spacing(1, 1, 1);
ds = dsb.Create(vdims, origin, spacing);
localBlockIndicesPortal.Set(localBlockIndex, vtkm::Id3(0, 0, blockIndex));
localBlockOriginsPortal.Set(localBlockIndex,
vtkm::Id3(0, 0, (blockStart / blockSliceSize)));
localBlockSizesPortal.Set(
localBlockIndex,
vtkm::Id3(static_cast<vtkm::Id>(dims[0]), static_cast<vtkm::Id>(dims[1]), currBlockSize));
vtkm::cont::CellSetStructured<3> cs;
cs.SetPointDimensions(vdims);
cs.SetGlobalPointDimensions(globalSize);
cs.SetGlobalPointIndexStart(vtkm::Id3{ 0, 0, blockIndex * blockSize });
ds.SetCellSet(cs);
}
std::vector<vtkm::Float32> subValues((values.begin() + blockStart),
@ -648,8 +641,7 @@ int main(int argc, char* argv[])
computeRegularStructure);
#ifdef WITH_MPI
filter.SetSpatialDecomposition(
blocksPerDim, globalSize, localBlockIndices, localBlockOrigins, localBlockSizes);
filter.SetBlockIndices(blocksPerDim, localBlockIndices);
#endif
filter.SetActiveField("values");

@ -401,14 +401,8 @@ int main(int argc, char* argv[])
vtkm::Id3 globalSize;
vtkm::Id3 blocksPerDim;
vtkm::cont::ArrayHandle<vtkm::Id3> localBlockIndices;
vtkm::cont::ArrayHandle<vtkm::Id3> localBlockOrigins;
vtkm::cont::ArrayHandle<vtkm::Id3> localBlockSizes;
localBlockIndices.Allocate(blocksPerRank);
localBlockOrigins.Allocate(blocksPerRank);
localBlockSizes.Allocate(blocksPerRank);
auto localBlockIndicesPortal = localBlockIndices.WritePortal();
auto localBlockOriginsPortal = localBlockOrigins.WritePortal();
auto localBlockSizesPortal = localBlockSizes.WritePortal();
// Read the pre-split data files
if (preSplitFiles)
@ -595,6 +589,11 @@ int main(int argc, char* argv[])
static_cast<ValueType>(offset[1]) };
const vtkm::Vec<ValueType, 2> v_spacing{ 1, 1 };
ds = dsb.Create(v_dims, v_origin, v_spacing);
vtkm::cont::CellSetStructured<2> cs;
cs.SetPointDimensions(v_dims);
cs.SetGlobalPointDimensions(vtkm::Id2{ globalSize[0], globalSize[1] });
cs.SetGlobalPointIndexStart(vtkm::Id2{ offset[0], offset[1] });
ds.SetCellSet(cs);
}
else
{
@ -607,6 +606,11 @@ int main(int argc, char* argv[])
static_cast<ValueType>(offset[2]) };
vtkm::Vec<ValueType, 3> v_spacing(1, 1, 1);
ds = dsb.Create(v_dims, v_origin, v_spacing);
vtkm::cont::CellSetStructured<3> cs;
cs.SetPointDimensions(v_dims);
cs.SetGlobalPointDimensions(globalSize);
cs.SetGlobalPointIndexStart(vtkm::Id3{ offset[0], offset[1], offset[2] });
ds.SetCellSet(cs);
}
ds.AddPointField("values", values);
// and add to partition
@ -617,14 +621,6 @@ int main(int argc, char* argv[])
vtkm::Id3{ static_cast<vtkm::Id>(blockIndex[0]),
static_cast<vtkm::Id>(blockIndex[1]),
static_cast<vtkm::Id>(nDims == 3 ? blockIndex[2] : 0) });
localBlockOriginsPortal.Set(blockNo,
vtkm::Id3{ static_cast<vtkm::Id>(offset[0]),
static_cast<vtkm::Id>(offset[1]),
static_cast<vtkm::Id>(nDims == 3 ? offset[2] : 0) });
localBlockSizesPortal.Set(blockNo,
vtkm::Id3{ static_cast<vtkm::Id>(dims[0]),
static_cast<vtkm::Id>(dims[1]),
static_cast<vtkm::Id>(nDims == 3 ? dims[2] : 0) });
if (blockNo == 0)
{
@ -753,7 +749,6 @@ int main(int argc, char* argv[])
: vtkm::Id3(static_cast<vtkm::Id>(dims[0]),
static_cast<vtkm::Id>(dims[1]),
static_cast<vtkm::Id>(1));
std::cout << blocksPerDim << " " << globalSize << std::endl;
{
vtkm::Id lastDimSize =
(nDims == 2) ? static_cast<vtkm::Id>(dims[1]) : static_cast<vtkm::Id>(dims[2]);
@ -801,12 +796,12 @@ int main(int argc, char* argv[])
vtkm::Vec<ValueType, 2> origin(0, blockIndex * blockSize);
vtkm::Vec<ValueType, 2> spacing(1, 1);
ds = dsb.Create(vdims, origin, spacing);
vtkm::cont::CellSetStructured<2> cs;
cs.SetPointDimensions(vdims);
cs.SetGlobalPointDimensions(vtkm::Id2{ globalSize[0], globalSize[1] });
cs.SetGlobalPointIndexStart(vtkm::Id2{ 0, (blockStart / blockSliceSize) });
ds.SetCellSet(cs);
localBlockIndicesPortal.Set(localBlockIndex, vtkm::Id3(0, blockIndex, 0));
localBlockOriginsPortal.Set(localBlockIndex,
vtkm::Id3(0, (blockStart / blockSliceSize), 0));
localBlockSizesPortal.Set(localBlockIndex,
vtkm::Id3(static_cast<vtkm::Id>(dims[0]), currBlockSize, 0));
}
// 3D data
else
@ -818,14 +813,12 @@ int main(int argc, char* argv[])
vtkm::Vec<ValueType, 3> origin(0, 0, (blockIndex * blockSize));
vtkm::Vec<ValueType, 3> spacing(1, 1, 1);
ds = dsb.Create(vdims, origin, spacing);
vtkm::cont::CellSetStructured<3> cs;
cs.SetPointDimensions(vdims);
cs.SetGlobalPointDimensions(globalSize);
cs.SetGlobalPointIndexStart(vtkm::Id3(0, 0, blockStart / blockSliceSize));
ds.SetCellSet(cs);
localBlockIndicesPortal.Set(localBlockIndex, vtkm::Id3(0, 0, blockIndex));
localBlockOriginsPortal.Set(localBlockIndex,
vtkm::Id3(0, 0, (blockStart / blockSliceSize)));
localBlockSizesPortal.Set(localBlockIndex,
vtkm::Id3(static_cast<vtkm::Id>(dims[0]),
static_cast<vtkm::Id>(dims[1]),
currBlockSize));
}
std::vector<vtkm::Float32> subValues((values.begin() + blockStart),
@ -863,13 +856,9 @@ int main(int argc, char* argv[])
prevTime = currTime;
// Convert the mesh of values into contour tree, pairs of vertex ids
vtkm::filter::ContourTreeUniformDistributed filter(blocksPerDim,
globalSize,
localBlockIndices,
localBlockOrigins,
localBlockSizes,
timingsLogLevel,
treeLogLevel);
vtkm::filter::scalar_topology::ContourTreeUniformDistributed filter(timingsLogLevel,
treeLogLevel);
filter.SetBlockIndices(blocksPerDim, localBlockIndices);
filter.SetUseBoundaryExtremaOnly(useBoundaryExtremaOnly);
filter.SetUseMarchingCubes(useMarchingCubes);
filter.SetAugmentHierarchicalTree(augmentHierarchicalTree);
@ -895,8 +884,7 @@ int main(int argc, char* argv[])
{
if (computeHierarchicalVolumetricBranchDecomposition)
{
vtkm::filter::scalar_topology::DistributedBranchDecompositionFilter bd_filter(
blocksPerDim, globalSize, localBlockIndices, localBlockOrigins, localBlockSizes);
vtkm::filter::scalar_topology::DistributedBranchDecompositionFilter bd_filter;
auto bd_result = bd_filter.Execute(result);
for (vtkm::Id ds_no = 0; ds_no < result.GetNumberOfPartitions(); ++ds_no)

@ -14,10 +14,4 @@ project(ParticleAdvection CXX)
find_package(VTKm REQUIRED QUIET)
add_executable(Particle_Advection ParticleAdvection.cxx)
target_link_libraries(Particle_Advection PRIVATE vtkm_filter vtkm_io)
vtkm_add_target_information(Particle_Advection
DROP_UNUSED_SYMBOLS MODIFY_CUDA_FLAGS
DEVICE_SOURCES ParticleAdvection.cxx)
if(TARGET vtkm::tbb)
target_compile_definitions(Particle_Advection PRIVATE BUILDING_TBB_VERSION)
endif()
target_link_libraries(Particle_Advection PRIVATE vtkm_filter_flow vtkm_io)

@ -8,9 +8,10 @@
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#include <vtkm/Particle.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/Initialize.h>
#include <vtkm/filter/Streamline.h>
#include <vtkm/filter/flow/Streamline.h>
#include <vtkm/io/VTKDataSetReader.h>
#include <vtkm/io/VTKDataSetWriter.h>
@ -74,7 +75,7 @@ int main(int argc, char** argv)
auto seedArray = vtkm::cont::make_ArrayHandle(seeds, vtkm::CopyFlag::Off);
//compute streamlines
vtkm::filter::Streamline streamline;
vtkm::filter::flow::Streamline streamline;
streamline.SetStepSize(stepSize);
streamline.SetNumberOfSteps(numSteps);

@ -16,12 +16,5 @@ find_package(VTKm REQUIRED QUIET)
if (VTKm_ENABLE_MPI)
add_executable(StreamlineMPI StreamlineMPI.cxx)
target_compile_definitions(StreamlineMPI PRIVATE "MPI_ENABLED")
target_link_libraries(StreamlineMPI PRIVATE vtkm_filter vtkm_io MPI::MPI_CXX)
vtkm_add_target_information(StreamlineMPI
DROP_UNUSED_SYMBOLS MODIFY_CUDA_FLAGS
DEVICE_SOURCES StreamlineMPI.cxx)
target_link_libraries(StreamlineMPI PRIVATE vtkm_filter_flow vtkm_io MPI::MPI_CXX)
endif()
#if(TARGET vtkm::tbb)
# target_compile_definitions(streamline_mpi PRIVATE BUILDING_TBB_VERSION)
#endif()

@ -8,13 +8,14 @@
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#include <vtkm/Particle.h>
#include <vtkm/cont/AssignerPartitionedDataSet.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/EnvironmentTracker.h>
#include <vtkm/cont/Field.h>
#include <vtkm/cont/Initialize.h>
#include <vtkm/cont/PartitionedDataSet.h>
#include <vtkm/filter/Streamline.h>
#include <vtkm/filter/flow/ParticleAdvection.h>
#include <vtkm/io/VTKDataSetReader.h>
#include <vtkm/io/VTKDataSetWriter.h>
#include <vtkm/io/reader/VTKDataSetReader.h>
@ -23,12 +24,6 @@
#include <vtkm/thirdparty/diy/diy.h>
#include <vtkm/thirdparty/diy/mpi-cast.h>
#include <vtkm/filter/ParticleAdvection.h>
#include <vtkm/filter/particleadvection/BoundsMap.h>
#include <vtkm/filter/particleadvection/ParticleMessenger.h>
void LoadData(std::string& fname, std::vector<vtkm::cont::DataSet>& dataSets, int rank, int nRanks)
{
std::string buff;
@ -99,7 +94,7 @@ int main(int argc, char** argv)
std::vector<vtkm::cont::DataSet> dataSets;
LoadData(dataFile, dataSets, rank, size);
vtkm::filter::ParticleAdvection pa;
vtkm::filter::flow::ParticleAdvection pa;
vtkm::cont::ArrayHandle<vtkm::Particle> seedArray;
seedArray = vtkm::cont::make_ArrayHandle({ vtkm::Particle(vtkm::Vec3f(.1f, .1f, .9f), 0),

@ -12,11 +12,14 @@
#include <string>
#include <vector>
#include <vtkm/Particle.h>
//#include <vtkm/cont/Algorithm.h>
#include <vtkm/cont/ArrayCopy.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/Initialize.h>
#include <vtkm/cont/Timer.h>
//#include <vtkm/cont/Timer.h>
#include <vtkm/filter/Pathline.h>
#include <vtkm/filter/flow/Pathline.h>
#include <vtkm/io/VTKDataSetReader.h>
#include <vtkm/io/VTKDataSetWriter.h>
@ -91,7 +94,7 @@ int main(int argc, char** argv)
// Instantiate the filter by providing necessary parameters.
// Necessary parameters are :
vtkm::filter::Pathline pathlineFilter;
vtkm::filter::flow::Pathline pathlineFilter;
pathlineFilter.SetActiveField(fieldName);
// 1. The current and next time slice. The current time slice is passed
// through the parameter to the Execute method.

@ -142,6 +142,17 @@ public:
vtkm::Id NumSteps = 0;
vtkm::ParticleStatus Status;
vtkm::FloatDefault Time = 0;
static size_t Sizeof()
{
constexpr std::size_t sz = sizeof(vtkm::Vec3f) // Pos
+ sizeof(vtkm::Id) // ID
+ sizeof(vtkm::Id) // NumSteps
+ sizeof(vtkm::UInt8) // Status
+ sizeof(vtkm::FloatDefault); // Time
return sz;
}
};
class ChargedParticle
@ -228,6 +239,14 @@ public:
return this->Pos + translation;
}
inline VTKM_CONT friend std::ostream& operator<<(std::ostream& out,
const vtkm::ChargedParticle& p)
{
out << "v(" << p.Time << ") = " << p.Pos << ", ID: " << p.ID << ", NumSteps: " << p.NumSteps
<< ", Status: " << p.Status;
return out;
}
vtkm::Vec3f Pos;
vtkm::Id ID = -1;
vtkm::Id NumSteps = 0;
@ -243,6 +262,23 @@ private:
static_cast<vtkm::FloatDefault>(2.99792458e8);
friend struct mangled_diy_namespace::Serialization<vtkm::ChargedParticle>;
public:
static size_t Sizeof()
{
constexpr std::size_t sz = sizeof(vtkm::Vec3f) // Pos
+ sizeof(vtkm::Id) // ID
+ sizeof(vtkm::Id) // NumSteps
+ sizeof(vtkm::UInt8) // Status
+ sizeof(vtkm::FloatDefault) // Time
+ sizeof(vtkm::FloatDefault) //Mass
+ sizeof(vtkm::FloatDefault) //Charge
+ sizeof(vtkm::FloatDefault) //Weighting
+ sizeof(vtkm::Vec3f) //Momentum
+ sizeof(vtkm::FloatDefault); //Speed_of_light
return sz;
}
};
} //namespace vtkm

@ -51,6 +51,11 @@ public:
this->Structure.SetPointDimensions(dimensions);
}
void SetGlobalPointDimensions(SchedulingRangeType dimensions)
{
this->Structure.SetGlobalPointDimensions(dimensions);
}
void SetGlobalPointIndexStart(SchedulingRangeType start)
{
this->Structure.SetGlobalPointIndexStart(start);
@ -58,8 +63,18 @@ public:
SchedulingRangeType GetPointDimensions() const { return this->Structure.GetPointDimensions(); }
SchedulingRangeType GetGlobalPointDimensions() const
{
return this->Structure.GetGlobalPointDimensions();
}
SchedulingRangeType GetCellDimensions() const { return this->Structure.GetCellDimensions(); }
SchedulingRangeType GetGlobalCellDimensions() const
{
return this->Structure.GetGlobalCellDimensions();
}
SchedulingRangeType GetGlobalPointIndexStart() const
{
return this->Structure.GetGlobalPointIndexStart();
@ -225,17 +240,20 @@ public:
static VTKM_CONT void save(BinaryBuffer& bb, const Type& cs)
{
vtkmdiy::save(bb, cs.GetPointDimensions());
vtkmdiy::save(bb, cs.GetGlobalPointDimensions());
vtkmdiy::save(bb, cs.GetGlobalPointIndexStart());
}
static VTKM_CONT void load(BinaryBuffer& bb, Type& cs)
{
typename Type::SchedulingRangeType dims, start;
typename Type::SchedulingRangeType dims, gdims, start;
vtkmdiy::load(bb, dims);
vtkmdiy::load(bb, gdims);
vtkmdiy::load(bb, start);
cs = Type{};
cs.SetPointDimensions(dims);
cs.SetGlobalPointDimensions(gdims);
cs.SetGlobalPointIndexStart(start);
}
};

@ -14,6 +14,23 @@ namespace vtkm
{
namespace cont
{
VTKM_CONT std::string& GlobalGhostCellFieldName() noexcept
{
static std::string GhostCellName("vtkGhostCells");
return GhostCellName;
}
VTKM_CONT const std::string& GetGlobalGhostCellFieldName() noexcept
{
return GlobalGhostCellFieldName();
}
VTKM_CONT void SetGlobalGhostCellFieldName(const std::string& name) noexcept
{
GlobalGhostCellFieldName() = name;
}
void DataSet::Clear()
{
this->CoordSystems.clear();
@ -39,6 +56,7 @@ void DataSet::CopyStructure(const vtkm::cont::DataSet& source)
{
this->CoordSystems = source.CoordSystems;
this->CellSet = source.CellSet;
this->GhostCellName = source.GhostCellName;
}
const vtkm::cont::Field& DataSet::GetField(vtkm::Id index) const
@ -59,16 +77,35 @@ vtkm::cont::Field& DataSet::GetField(vtkm::Id index)
vtkm::Id DataSet::GetFieldIndex(const std::string& name, vtkm::cont::Field::Association assoc) const
{
bool found;
vtkm::Id index = this->FindFieldIndex(name, assoc, found);
if (found)
const auto it = this->Fields.find({ name, assoc });
if (it != this->Fields.end())
{
return index;
return static_cast<vtkm::Id>(std::distance(this->Fields.begin(), it));
}
else
return -1;
}
const vtkm::cont::Field& DataSet::GetField(const std::string& name,
vtkm::cont::Field::Association assoc) const
{
auto idx = this->GetFieldIndex(name, assoc);
if (idx == -1)
{
throw vtkm::cont::ErrorBadValue("No field with requested name: " + name);
}
return this->GetField(idx);
}
vtkm::cont::Field& DataSet::GetField(const std::string& name, vtkm::cont::Field::Association assoc)
{
auto idx = this->GetFieldIndex(name, assoc);
if (idx == -1)
{
throw vtkm::cont::ErrorBadValue("No field with requested name: " + name);
}
return this->GetField(idx);
}
const vtkm::cont::CoordinateSystem& DataSet::GetCoordinateSystem(vtkm::Id index) const
@ -150,21 +187,6 @@ void DataSet::PrintSummary(std::ostream& out) const
out.flush();
}
vtkm::Id DataSet::FindFieldIndex(const std::string& name,
vtkm::cont::Field::Association association,
bool& found) const
{
const auto it = this->Fields.find({ name, association });
if (it != this->Fields.end())
{
found = true;
return static_cast<vtkm::Id>(std::distance(this->Fields.begin(), it));
}
found = false;
return -1;
}
VTKM_CONT void DataSet::AddField(const Field& field)
{
// map::insert will not replace the duplicated element

@ -24,9 +24,25 @@ namespace vtkm
namespace cont
{
VTKM_CONT_EXPORT VTKM_CONT std::string& GlobalGhostCellFieldName() noexcept;
VTKM_CONT_EXPORT VTKM_CONT const std::string& GetGlobalGhostCellFieldName() noexcept;
VTKM_CONT_EXPORT VTKM_CONT void SetGlobalGhostCellFieldName(const std::string& name) noexcept;
class VTKM_CONT_EXPORT DataSet
{
public:
DataSet() = default;
DataSet(vtkm::cont::DataSet&&) = default;
DataSet(const vtkm::cont::DataSet&) = default;
vtkm::cont::DataSet& operator=(vtkm::cont::DataSet&&) = default;
vtkm::cont::DataSet& operator=(const vtkm::cont::DataSet&) = default;
VTKM_CONT void Clear();
/// Get the number of cells contained in this DataSet
@ -50,25 +66,45 @@ public:
bool HasField(const std::string& name,
vtkm::cont::Field::Association assoc = vtkm::cont::Field::Association::Any) const
{
bool found = false;
this->FindFieldIndex(name, assoc, found);
return found;
return (this->GetFieldIndex(name, assoc) != -1);
}
VTKM_CONT
bool HasCellField(const std::string& name) const
{
bool found = false;
this->FindFieldIndex(name, vtkm::cont::Field::Association::Cells, found);
return found;
return (this->GetFieldIndex(name, vtkm::cont::Field::Association::Cells) != -1);
}
VTKM_CONT
bool HasGhostCellField() const
{
if (this->GhostCellName)
{
return this->HasCellField(*this->GhostCellName);
}
else
{
return false;
}
}
VTKM_CONT
const std::string& GetGhostCellFieldName() const
{
if (this->HasGhostCellField())
{
return *this->GhostCellName;
}
else
{
throw vtkm::cont::ErrorBadValue("No Ghost Cell Field Name");
}
}
VTKM_CONT
bool HasPointField(const std::string& name) const
{
bool found = false;
this->FindFieldIndex(name, vtkm::cont::Field::Association::Points, found);
return found;
return (this->GetFieldIndex(name, vtkm::cont::Field::Association::Points) != -1);
}
@ -85,18 +121,12 @@ public:
VTKM_CONT
const vtkm::cont::Field& GetField(
const std::string& name,
vtkm::cont::Field::Association assoc = vtkm::cont::Field::Association::Any) const
{
return this->GetField(this->GetFieldIndex(name, assoc));
}
vtkm::cont::Field::Association assoc = vtkm::cont::Field::Association::Any) const;
VTKM_CONT
vtkm::cont::Field& GetField(
const std::string& name,
vtkm::cont::Field::Association assoc = vtkm::cont::Field::Association::Any)
{
return this->GetField(this->GetFieldIndex(name, assoc));
}
vtkm::cont::Field::Association assoc = vtkm::cont::Field::Association::Any);
//@}
/// Returns the first cell field that matches the provided name.
@ -115,6 +145,23 @@ public:
}
//@}
/// Returns the first cell field that matches the provided name.
/// Will throw an exception if no match is found
//@{
VTKM_CONT
const vtkm::cont::Field& GetGhostCellField() const
{
if (this->HasGhostCellField())
{
return this->GetCellField(*(this->GhostCellName));
}
else
{
throw vtkm::cont::ErrorBadValue("No Ghost Cell Field");
}
}
//@}
/// Returns the first point field that matches the provided name.
/// Will throw an exception if no match is found
//@{
@ -165,6 +212,25 @@ public:
this->AddField(make_FieldCell(fieldName, field));
}
VTKM_CONT
void AddGhostCellField(const std::string& fieldName, const vtkm::cont::UnknownArrayHandle& field)
{
this->GhostCellName.reset(new std::string(fieldName));
this->AddField(make_FieldCell(fieldName, field));
}
VTKM_CONT
void AddGhostCellField(const vtkm::cont::UnknownArrayHandle& field)
{
this->AddGhostCellField(GetGlobalGhostCellFieldName(), field);
}
VTKM_CONT
void AddGhostCellField(const vtkm::cont::Field& field)
{
this->AddGhostCellField(field.GetName(), field.GetData());
}
template <typename T, typename Storage>
VTKM_CONT void AddCellField(const std::string& fieldName,
const vtkm::cont::ArrayHandle<T, Storage>& field)
@ -296,11 +362,7 @@ private:
std::vector<vtkm::cont::CoordinateSystem> CoordSystems;
std::map<FieldCompare::Key, vtkm::cont::Field, FieldCompare> Fields;
vtkm::cont::UnknownCellSet CellSet;
VTKM_CONT
vtkm::Id FindFieldIndex(const std::string& name,
vtkm::cont::Field::Association association,
bool& found) const;
std::shared_ptr<std::string> GhostCellName;
};
} // namespace cont

@ -15,6 +15,13 @@
#include <vtkm/cont/ErrorBadValue.h>
#include <vtkm/cont/Field.h>
#include <vtkm/cont/PartitionedDataSet.h>
#include <vtkm/internal/Configure.h>
// clang-format off
VTKM_THIRDPARTY_PRE_INCLUDE
#include <vtkm/thirdparty/diy/diy.h>
VTKM_THIRDPARTY_POST_INCLUDE
// clang-format on
namespace vtkm
{
@ -73,6 +80,19 @@ vtkm::Id PartitionedDataSet::GetNumberOfPartitions() const
return static_cast<vtkm::Id>(this->Partitions.size());
}
VTKM_CONT
vtkm::Id PartitionedDataSet::GetGlobalNumberOfPartitions() const
{
#ifdef VTKM_ENABLE_MPI
auto comm = vtkm::cont::EnvironmentTracker::GetCommunicator();
vtkm::Id globalSize = 0;
vtkmdiy::mpi::all_reduce(comm, GetNumberOfPartitions(), globalSize, std::plus<vtkm::Id>{});
return globalSize;
#else
return GetNumberOfPartitions();
#endif
}
VTKM_CONT
const vtkm::cont::DataSet& PartitionedDataSet::GetPartition(vtkm::Id blockId) const
{

@ -32,16 +32,16 @@ public:
using reference = typename StorageVec::reference;
using const_reference = typename StorageVec::const_reference;
/// create a new PartitionedDataSet containng a single DataSet @a ds
/// Create a new PartitionedDataSet containng a single DataSet @a ds.
VTKM_CONT
PartitionedDataSet(const vtkm::cont::DataSet& ds);
/// create a new PartitionedDataSet with the existing one @a src
/// Create a new PartitionedDataSet with the existing one @a src.
VTKM_CONT
PartitionedDataSet(const vtkm::cont::PartitionedDataSet& src);
/// create a new PartitionedDataSet with a DataSet vector @a partitions.
/// Create a new PartitionedDataSet with a DataSet vector @a partitions.
VTKM_CONT
explicit PartitionedDataSet(const std::vector<vtkm::cont::DataSet>& partitions);
/// create a new PartitionedDataSet with the capacity set to be @a size
/// Create a new PartitionedDataSet with the capacity set to be @a size.
VTKM_CONT
explicit PartitionedDataSet(vtkm::Id size);
@ -53,33 +53,41 @@ public:
VTKM_CONT
~PartitionedDataSet();
/// get the field @a field_name from partition @a partition_index
/// Get the field @a field_name from partition @a partition_index.
VTKM_CONT
vtkm::cont::Field GetField(const std::string& field_name, int partition_index) const;
/// Get number of DataSet objects stored in this PartitionedDataSet.
VTKM_CONT
vtkm::Id GetNumberOfPartitions() const;
/// Get number of partations across all MPI ranks.
/// @warning This method requires global communication (MPI_Allreduce) if MPI is enabled.
VTKM_CONT
vtkm::Id GetGlobalNumberOfPartitions() const;
/// Get the DataSet @a partId.
VTKM_CONT
const vtkm::cont::DataSet& GetPartition(vtkm::Id partId) const;
/// Get an STL vector of all DataSet objects stored in PartitionedDataSet.
VTKM_CONT
const std::vector<vtkm::cont::DataSet>& GetPartitions() const;
/// add DataSet @a ds to the end of the contained DataSet vector
/// Add DataSet @a ds to the end of the contained DataSet vector.
VTKM_CONT
void AppendPartition(const vtkm::cont::DataSet& ds);
/// add DataSet @a ds to position @a index of the contained DataSet vector
/// Add DataSet @a ds to position @a index of the contained DataSet vector.
VTKM_CONT
void InsertPartition(vtkm::Id index, const vtkm::cont::DataSet& ds);
/// replace the @a index positioned element of the contained DataSet vector
/// with @a ds
/// Replace the @a index positioned element of the contained DataSet vector
/// with @a ds.
VTKM_CONT
void ReplacePartition(vtkm::Id index, const vtkm::cont::DataSet& ds);
/// append the DataSet vector "partitions" to the end of the contained one
/// Append the DataSet vector @a partitions to the end of the contained one.
VTKM_CONT
void AppendPartitions(const std::vector<vtkm::cont::DataSet>& partitions);

@ -95,7 +95,17 @@ bool UnknownAHComponentInfo::operator==(const UnknownAHComponentInfo& rhs)
else
{
// Needs optimization based on platform. OSX cannot compare typeid across translation units?
return this->Type == rhs.Type;
bool typesEqual = (this->Type == rhs.Type);
// Could use optimization based on platform. OSX cannot compare typeid across translation
// units, so we have to also check the names. (Why doesn't the == operator above do that?)
// Are there other platforms that behave similarly?
if (!typesEqual)
{
typesEqual = (std::strcmp(this->Type.name(), rhs.Type.name()) == 0);
}
return typesEqual;
}
}
@ -151,7 +161,7 @@ VTKM_CONT bool UnknownArrayHandle::IsBaseComponentTypeImpl(
return false;
}
// Needs optimization based on platform. OSX cannot compare typeid across translation units?
// Note that detail::UnknownAHComponentInfo has a custom operator==
return this->Container->BaseComponentType == type;
}

@ -9,6 +9,7 @@
//============================================================================
#include <random>
#include <vtkm/cont/Algorithm.h>
#include <vtkm/cont/ParticleArrayCopy.h>
#include <vtkm/cont/testing/Testing.h>

@ -21,6 +21,7 @@ namespace test_variant
template <vtkm::IdComponent Index>
struct TypePlaceholder
{
vtkm::IdComponent Value = Index;
};
// A class that is trivially copiable but not totally trivial.
@ -313,6 +314,17 @@ struct TestFunctor
}
};
struct TestFunctorModify
{
template <vtkm::IdComponent Index>
void operator()(TypePlaceholder<Index>& object, vtkm::Id expectedValue)
{
VTKM_TEST_ASSERT(Index == expectedValue, "Index = ", Index, ", expected = ", expectedValue);
VTKM_TEST_ASSERT(object.Value == expectedValue);
++object.Value;
}
};
void TestGet()
{
std::cout << "Test Get" << std::endl;
@ -431,6 +443,11 @@ void TestCastAndCall()
VariantType variant26{ TypePlaceholder<26>{} };
result = variant26.CastAndCall(TestFunctor(), 26);
VTKM_TEST_ASSERT(test_equal(result, TestValue(26, vtkm::FloatDefault{})));
variant1.CastAndCall(TestFunctorModify{}, 1);
VTKM_TEST_ASSERT(variant1.Get<1>().Value == 2, "Variant object not updated.");
variant1.CastAndCall([](auto& object) { ++object.Value; });
}
void TestCopyInvalid()

@ -81,11 +81,11 @@ set(common_headers
FilterDataSetWithField.h
FilterField.h
Filter.h
FilterParticleAdvection.h
FilterTemporalParticleAdvection.h
FilterTraits.h
PolicyBase.h
PolicyDefault.h
Lagrangian.h
LagrangianStructures.h
)
vtkm_declare_headers(${common_headers})
@ -95,37 +95,12 @@ set(common_header_template_sources
FilterDataSetWithField.hxx
FilterField.hxx
Filter.hxx
FilterParticleAdvection.hxx
FilterTemporalParticleAdvection.hxx
Lagrangian.hxx
LagrangianStructures.hxx
)
vtkm_declare_headers(${common_header_template_sources})
set(extra_headers
Lagrangian.h
LagrangianStructures.h
ParticleAdvection.h
Pathline.h
PathParticle.h
Streamline.h
StreamSurface.h
)
set(extra_header_template_sources
Lagrangian.hxx
LagrangianStructures.hxx
ParticleAdvection.hxx
Pathline.hxx
PathParticle.hxx
Streamline.hxx
StreamSurface.hxx
)
set(extra_sources_device
particleadvection/Messenger.cxx
particleadvection/ParticleMessenger.cxx
)
set(core_headers
FieldSelection.h
NewFilter.h
@ -153,30 +128,16 @@ vtkm_library(
add_library(vtkm_filter INTERFACE)
vtkm_library(
NAME vtkm_filter_extra
TEMPLATE_SOURCES ${extra_header_template_sources}
HEADERS ${extra_headers}
DEVICE_SOURCES ${extra_sources_device}
USE_VTKM_JOB_POOL
)
set_target_properties(
vtkm_filter_core
vtkm_filter_extra
PROPERTIES
UNITY_BUILD ON
UNITY_BUILD_MODE BATCH
)
target_link_libraries(vtkm_filter_core PUBLIC vtkm_cont vtkm_worklet)
target_link_libraries(vtkm_filter_extra PUBLIC vtkm_cont vtkm_worklet)
if (VTKm_ENABLE_MPI)
target_link_libraries(vtkm_filter_extra PUBLIC MPI::MPI_CXX)
endif()
target_link_libraries(vtkm_filter PUBLIC INTERFACE
vtkm_filter_extra
vtkm_filter_core
)
@ -187,9 +148,9 @@ add_subdirectory(connected_components)
add_subdirectory(contour)
add_subdirectory(density_estimate)
add_subdirectory(entity_extraction)
add_subdirectory(flow)
add_subdirectory(image_processing)
add_subdirectory(internal)
add_subdirectory(particleadvection)
add_subdirectory(field_conversion)
add_subdirectory(field_transform)
add_subdirectory(geometry_refinement)

@ -1,92 +0,0 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
//
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#ifndef vtk_m_filter_FilterParticleAdvection_h
#define vtk_m_filter_FilterParticleAdvection_h
#include <vtkm/Particle.h>
#include <vtkm/filter/FilterDataSetWithField.h>
#include <vtkm/filter/particleadvection/BoundsMap.h>
#include <vtkm/filter/particleadvection/DataSetIntegrator.h>
namespace vtkm
{
namespace filter
{
/// \brief base class for advecting particles in a vector field.
/// Takes as input a vector field and seed locations and advects the seeds
/// through the flow field.
template <class Derived, typename ParticleType>
class FilterParticleAdvection : public vtkm::filter::FilterDataSetWithField<Derived>
{
public:
using SupportedTypes = vtkm::TypeListFieldVec3;
VTKM_CONT
FilterParticleAdvection();
VTKM_CONT
void SetStepSize(vtkm::FloatDefault s) { this->StepSize = s; }
VTKM_CONT
void SetNumberOfSteps(vtkm::Id n) { this->NumberOfSteps = n; }
VTKM_CONT
void SetSeeds(const std::vector<ParticleType>& seeds,
vtkm::CopyFlag copyFlag = vtkm::CopyFlag::On)
{
this->Seeds = vtkm::cont::make_ArrayHandle(seeds, copyFlag);
}
VTKM_CONT
void SetSeeds(vtkm::cont::ArrayHandle<ParticleType>& seeds) { this->Seeds = seeds; }
VTKM_CONT
bool GetUseThreadedAlgorithm() { return this->UseThreadedAlgorithm; }
VTKM_CONT
void SetUseThreadedAlgorithm(bool val) { this->UseThreadedAlgorithm = val; }
template <typename DerivedPolicy>
VTKM_CONT vtkm::cont::DataSet PrepareForExecution(const vtkm::cont::DataSet& input,
vtkm::filter::PolicyBase<DerivedPolicy> policy);
template <typename DerivedPolicy>
VTKM_CONT bool MapFieldOntoOutput(vtkm::cont::DataSet&,
const vtkm::cont::Field&,
vtkm::filter::PolicyBase<DerivedPolicy>)
{
return false;
}
protected:
VTKM_CONT virtual void ValidateOptions() const;
VTKM_CONT std::vector<vtkm::filter::particleadvection::DataSetIntegrator>
CreateDataSetIntegrators(const vtkm::cont::PartitionedDataSet& input,
const vtkm::filter::particleadvection::BoundsMap& boundsMap) const;
vtkm::Id NumberOfSteps;
vtkm::FloatDefault StepSize;
vtkm::cont::ArrayHandle<ParticleType> Seeds;
bool UseThreadedAlgorithm;
private:
};
}
} // namespace vtkm::filter
#ifndef vtk_m_filter_FilterParticleAdvection_hxx
#include <vtkm/filter/FilterParticleAdvection.hxx>
#endif
#endif // vtk_m_filter_FilterParticleAdvection_h

@ -1,83 +0,0 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
//
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#ifndef vtk_m_filter_FilterParticleAdvection_hxx
#define vtk_m_filter_FilterParticleAdvection_hxx
#include <vtkm/filter/FilterParticleAdvection.h>
namespace vtkm
{
namespace filter
{
//-----------------------------------------------------------------------------
template <typename Derived, typename ParticleType>
inline VTKM_CONT FilterParticleAdvection<Derived, ParticleType>::FilterParticleAdvection()
: vtkm::filter::FilterDataSetWithField<Derived>()
, NumberOfSteps(0)
, StepSize(0)
, UseThreadedAlgorithm(false)
{
}
template <typename Derived, typename ParticleType>
void FilterParticleAdvection<Derived, ParticleType>::ValidateOptions() const
{
if (this->GetUseCoordinateSystemAsField())
throw vtkm::cont::ErrorFilterExecution("Coordinate system as field not supported");
if (this->Seeds.GetNumberOfValues() == 0)
throw vtkm::cont::ErrorFilterExecution("No seeds provided.");
if (this->NumberOfSteps == 0)
throw vtkm::cont::ErrorFilterExecution("Number of steps not specified.");
if (this->StepSize == 0)
throw vtkm::cont::ErrorFilterExecution("Step size not specified.");
}
template <typename Derived, typename ParticleType>
std::vector<vtkm::filter::particleadvection::DataSetIntegrator>
FilterParticleAdvection<Derived, ParticleType>::CreateDataSetIntegrators(
const vtkm::cont::PartitionedDataSet& input,
const vtkm::filter::particleadvection::BoundsMap& boundsMap) const
{
using DSIType = vtkm::filter::particleadvection::DataSetIntegrator;
std::vector<DSIType> dsi;
if (boundsMap.GetTotalNumBlocks() == 0)
throw vtkm::cont::ErrorFilterExecution("No input datasets.");
std::string activeField = this->GetActiveFieldName();
for (vtkm::Id i = 0; i < input.GetNumberOfPartitions(); i++)
{
vtkm::Id blockId = boundsMap.GetLocalBlockId(i);
auto ds = input.GetPartition(i);
if (!ds.HasPointField(activeField) && !ds.HasCellField(activeField))
throw vtkm::cont::ErrorFilterExecution("Unsupported field assocation");
dsi.push_back(DSIType(ds, blockId, activeField));
}
return dsi;
}
//-----------------------------------------------------------------------------
template <typename Derived, typename ParticleType>
template <typename DerivedPolicy>
inline VTKM_CONT vtkm::cont::DataSet
FilterParticleAdvection<Derived, ParticleType>::PrepareForExecution(
const vtkm::cont::DataSet& input,
vtkm::filter::PolicyBase<DerivedPolicy> policy)
{
return (static_cast<Derived*>(this))->DoExecute(input, policy);
}
}
} // namespace vtkm::filter
#endif

@ -1,75 +0,0 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
//
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#ifndef vtk_m_filter_TemporalFilterParticleAdvection_h
#define vtk_m_filter_TemporalFilterParticleAdvection_h
#include <vtkm/Particle.h>
#include <vtkm/filter/FilterParticleAdvection.h>
#include <vtkm/filter/particleadvection/BoundsMap.h>
#include <vtkm/filter/particleadvection/DataSetIntegrator.h>
namespace vtkm
{
namespace filter
{
/// \brief base class for advecting particles in a vector field.
/// Takes as input a vector field and seed locations and advects the seeds
/// through the flow field.
template <class Derived, typename ParticleType>
class FilterTemporalParticleAdvection
: public vtkm::filter::FilterParticleAdvection<Derived, ParticleType>
{
public:
VTKM_CONT
FilterTemporalParticleAdvection();
VTKM_CONT
void SetPreviousTime(vtkm::FloatDefault t) { this->PreviousTime = t; }
VTKM_CONT
void SetNextTime(vtkm::FloatDefault t) { this->NextTime = t; }
VTKM_CONT
void SetNextDataSet(const vtkm::cont::DataSet& ds)
{
this->NextDataSet = vtkm::cont::PartitionedDataSet(ds);
}
VTKM_CONT
void SetNextDataSet(const vtkm::cont::PartitionedDataSet& pds) { this->NextDataSet = pds; }
template <typename DerivedPolicy>
VTKM_CONT vtkm::cont::DataSet PrepareForExecution(const vtkm::cont::DataSet& input,
vtkm::filter::PolicyBase<DerivedPolicy> policy);
protected:
VTKM_CONT void ValidateOptions(const vtkm::cont::PartitionedDataSet& input) const;
using vtkm::filter::FilterParticleAdvection<Derived, ParticleType>::ValidateOptions;
VTKM_CONT std::vector<vtkm::filter::particleadvection::TemporalDataSetIntegrator>
CreateDataSetIntegrators(const vtkm::cont::PartitionedDataSet& input,
const vtkm::filter::particleadvection::BoundsMap& boundsMap) const;
vtkm::FloatDefault PreviousTime;
vtkm::FloatDefault NextTime;
vtkm::cont::PartitionedDataSet NextDataSet;
private:
};
}
} // namespace vtkm::filter
#ifndef vtk_m_filter_FilterTemporalParticleAdvection_hxx
#include <vtkm/filter/FilterTemporalParticleAdvection.hxx>
#endif
#endif // vtk_m_filter_FilterTemporalParticleAdvection_h

@ -1,83 +0,0 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
//
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#ifndef vtk_m_filter_FilterTemporalParticleAdvection_hxx
#define vtk_m_filter_FilterTemporalParticleAdvection_hxx
#include <vtkm/filter/FilterTemporalParticleAdvection.h>
namespace vtkm
{
namespace filter
{
//-----------------------------------------------------------------------------
template <typename Derived, typename ParticleType>
inline VTKM_CONT
FilterTemporalParticleAdvection<Derived, ParticleType>::FilterTemporalParticleAdvection()
: vtkm::filter::FilterParticleAdvection<Derived, ParticleType>()
, PreviousTime(0)
, NextTime(0)
{
}
template <typename Derived, typename ParticleType>
void FilterTemporalParticleAdvection<Derived, ParticleType>::ValidateOptions(
const vtkm::cont::PartitionedDataSet& input) const
{
this->vtkm::filter::FilterParticleAdvection<Derived, ParticleType>::ValidateOptions();
if (this->NextDataSet.GetNumberOfPartitions() != input.GetNumberOfPartitions())
throw vtkm::cont::ErrorFilterExecution("Number of partitions do not match");
if (!(this->PreviousTime < this->NextTime))
throw vtkm::cont::ErrorFilterExecution("Previous time must be less than Next time.");
}
template <typename Derived, typename ParticleType>
std::vector<vtkm::filter::particleadvection::TemporalDataSetIntegrator>
FilterTemporalParticleAdvection<Derived, ParticleType>::CreateDataSetIntegrators(
const vtkm::cont::PartitionedDataSet& input,
const vtkm::filter::particleadvection::BoundsMap& boundsMap) const
{
using DSIType = vtkm::filter::particleadvection::TemporalDataSetIntegrator;
std::vector<DSIType> dsi;
std::string activeField = this->GetActiveFieldName();
if (boundsMap.GetTotalNumBlocks() == 0)
throw vtkm::cont::ErrorFilterExecution("No input datasets.");
for (vtkm::Id i = 0; i < input.GetNumberOfPartitions(); i++)
{
vtkm::Id blockId = boundsMap.GetLocalBlockId(i);
auto dsPrev = input.GetPartition(i);
auto dsNext = this->NextDataSet.GetPartition(i);
if (!(dsPrev.HasPointField(activeField) || dsPrev.HasCellField(activeField)) ||
!(dsNext.HasPointField(activeField) || dsNext.HasCellField(activeField)))
throw vtkm::cont::ErrorFilterExecution("Unsupported field assocation");
dsi.push_back(
DSIType(dsPrev, this->PreviousTime, dsNext, this->NextTime, blockId, activeField));
}
return dsi;
}
//-----------------------------------------------------------------------------
template <typename Derived, typename ParticleType>
template <typename DerivedPolicy>
inline VTKM_CONT vtkm::cont::DataSet
FilterTemporalParticleAdvection<Derived, ParticleType>::PrepareForExecution(
const vtkm::cont::DataSet& input,
vtkm::filter::PolicyBase<DerivedPolicy> policy)
{
return (static_cast<Derived*>(this))->DoExecute(input, policy);
}
}
} // namespace vtkm::filter
#endif

@ -21,13 +21,13 @@
#include <vtkm/cont/DataSetBuilderRectilinear.h>
#include <vtkm/cont/DeviceAdapter.h>
#include <vtkm/cont/ErrorFilterExecution.h>
#include <vtkm/worklet/ParticleAdvection.h>
#include <vtkm/filter/flow/worklet/Field.h>
#include <vtkm/filter/flow/worklet/GridEvaluators.h>
#include <vtkm/filter/flow/worklet/ParticleAdvection.h>
#include <vtkm/filter/flow/worklet/Particles.h>
#include <vtkm/filter/flow/worklet/RK4Integrator.h>
#include <vtkm/filter/flow/worklet/Stepper.h>
#include <vtkm/worklet/WorkletMapField.h>
#include <vtkm/worklet/particleadvection/Field.h>
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
#include <vtkm/worklet/particleadvection/Particles.h>
#include <vtkm/worklet/particleadvection/RK4Integrator.h>
#include <vtkm/worklet/particleadvection/Stepper.h>
#include <cstring>
#include <sstream>
@ -274,13 +274,13 @@ inline VTKM_CONT vtkm::cont::DataSet Lagrangian::DoExecute(
vtkm::Bounds bounds = input.GetCoordinateSystem().GetBounds();
using FieldHandle = vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>, StorageType>;
using FieldType = vtkm::worklet::particleadvection::VelocityField<FieldHandle>;
using GridEvalType = vtkm::worklet::particleadvection::GridEvaluator<FieldType>;
using RK4Type = vtkm::worklet::particleadvection::RK4Integrator<GridEvalType>;
using Stepper = vtkm::worklet::particleadvection::Stepper<RK4Type, GridEvalType>;
using FieldType = vtkm::worklet::flow::VelocityField<FieldHandle>;
using GridEvalType = vtkm::worklet::flow::GridEvaluator<FieldType>;
using RK4Type = vtkm::worklet::flow::RK4Integrator<GridEvalType>;
using Stepper = vtkm::worklet::flow::Stepper<RK4Type, GridEvalType>;
vtkm::worklet::ParticleAdvection particleadvection;
vtkm::worklet::ParticleAdvectionResult<vtkm::Particle> res;
vtkm::worklet::flow::ParticleAdvection particleadvection;
vtkm::worklet::flow::ParticleAdvectionResult<vtkm::Particle> res;
FieldType velocities(field, fieldMeta.GetAssociation());
GridEvalType gridEval(coords, cells, velocities);

@ -0,0 +1,9 @@
//============================================================================
// 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.
//============================================================================

@ -15,9 +15,9 @@
#include <vtkm/filter/FilterDataSetWithField.h>
#include <vtkm/worklet/ParticleAdvection.h>
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
#include <vtkm/worklet/particleadvection/Stepper.h>
#include <vtkm/filter/flow/worklet/GridEvaluators.h>
#include <vtkm/filter/flow/worklet/ParticleAdvection.h>
#include <vtkm/filter/flow/worklet/Stepper.h>
namespace vtkm
{

@ -14,11 +14,11 @@
#include <vtkm/cont/ArrayHandleIndex.h>
#include <vtkm/cont/ErrorFilterExecution.h>
#include <vtkm/cont/Invoker.h>
#include <vtkm/worklet/particleadvection/Field.h>
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
#include <vtkm/worklet/particleadvection/Particles.h>
#include <vtkm/worklet/particleadvection/RK4Integrator.h>
#include <vtkm/worklet/particleadvection/Stepper.h>
#include <vtkm/filter/flow/worklet/Field.h>
#include <vtkm/filter/flow/worklet/GridEvaluators.h>
#include <vtkm/filter/flow/worklet/Particles.h>
#include <vtkm/filter/flow/worklet/RK4Integrator.h>
#include <vtkm/filter/flow/worklet/Stepper.h>
#include <vtkm/worklet/LagrangianStructures.h>
@ -79,10 +79,10 @@ inline VTKM_CONT vtkm::cont::DataSet LagrangianStructures::DoExecute(
using Structured3DType = vtkm::cont::CellSetStructured<3>;
using FieldHandle = vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>, StorageType>;
using FieldType = vtkm::worklet::particleadvection::VelocityField<FieldHandle>;
using GridEvaluator = vtkm::worklet::particleadvection::GridEvaluator<FieldType>;
using IntegratorType = vtkm::worklet::particleadvection::RK4Integrator<GridEvaluator>;
using Stepper = vtkm::worklet::particleadvection::Stepper<IntegratorType, GridEvaluator>;
using FieldType = vtkm::worklet::flow::VelocityField<FieldHandle>;
using GridEvaluator = vtkm::worklet::flow::GridEvaluator<FieldType>;
using IntegratorType = vtkm::worklet::flow::RK4Integrator<GridEvaluator>;
using Stepper = vtkm::worklet::flow::Stepper<IntegratorType, GridEvaluator>;
vtkm::FloatDefault stepSize = this->GetStepSize();
vtkm::Id numberOfSteps = this->GetNumberOfSteps();
@ -135,8 +135,8 @@ inline VTKM_CONT vtkm::cont::DataSet LagrangianStructures::DoExecute(
FieldType velocities(field, fieldMeta.GetAssociation());
GridEvaluator evaluator(input.GetCoordinateSystem(), input.GetCellSet(), velocities);
Stepper integrator(evaluator, stepSize);
vtkm::worklet::ParticleAdvection particles;
vtkm::worklet::ParticleAdvectionResult<vtkm::Particle> advectionResult;
vtkm::worklet::flow::ParticleAdvection particles;
vtkm::worklet::flow::ParticleAdvectionResult<vtkm::Particle> advectionResult;
vtkm::cont::ArrayHandle<vtkm::Particle> advectionPoints;
invoke(detail::MakeParticles{}, lcsInputPoints, advectionPoints);
advectionResult = particles.Run(integrator, advectionPoints, numberOfSteps);

@ -1,47 +0,0 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
//
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#ifndef vtk_m_filter_ParticleAdvection_h
#define vtk_m_filter_ParticleAdvection_h
#include <vtkm/Particle.h>
#include <vtkm/filter/FilterParticleAdvection.h>
namespace vtkm
{
namespace filter
{
/// \brief Advect particles in a vector field.
/// Takes as input a vector field and seed locations and generates the
/// end points for each seed through the vector field.
template <typename ParticleType = vtkm::Particle>
class ParticleAdvectionBase
: public vtkm::filter::FilterParticleAdvection<ParticleAdvectionBase<ParticleType>, ParticleType>
{
public:
VTKM_CONT ParticleAdvectionBase();
template <typename DerivedPolicy>
vtkm::cont::PartitionedDataSet PrepareForExecution(
const vtkm::cont::PartitionedDataSet& input,
const vtkm::filter::PolicyBase<DerivedPolicy>& policy);
};
using ParticleAdvection = ParticleAdvectionBase<vtkm::Particle>;
}
} // namespace vtkm::filter
#ifndef vtk_m_filter_ParticleAdvection_hxx
#include <vtkm/filter/ParticleAdvection.hxx>
#endif
#endif // vtk_m_filter_ParticleAdvection_h

@ -1,57 +0,0 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
//
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#ifndef vtk_m_filter_ParticleAdvection_hxx
#define vtk_m_filter_ParticleAdvection_hxx
#include <vtkm/cont/ErrorFilterExecution.h>
#include <vtkm/filter/ParticleAdvection.h>
#include <vtkm/filter/particleadvection/BoundsMap.h>
#include <vtkm/filter/particleadvection/DataSetIntegrator.h>
#include <vtkm/filter/particleadvection/ParticleAdvectionAlgorithm.h>
namespace vtkm
{
namespace filter
{
//-----------------------------------------------------------------------------
template <typename ParticleType>
inline VTKM_CONT ParticleAdvectionBase<ParticleType>::ParticleAdvectionBase()
: vtkm::filter::FilterParticleAdvection<ParticleAdvectionBase<ParticleType>, ParticleType>()
{
}
//-----------------------------------------------------------------------------
template <typename ParticleType>
template <typename DerivedPolicy>
inline VTKM_CONT vtkm::cont::PartitionedDataSet
ParticleAdvectionBase<ParticleType>::PrepareForExecution(
const vtkm::cont::PartitionedDataSet& input,
const vtkm::filter::PolicyBase<DerivedPolicy>&)
{
using AlgorithmType = vtkm::filter::particleadvection::ParticleAdvectionAlgorithm;
using ThreadedAlgorithmType = vtkm::filter::particleadvection::ParticleAdvectionThreadedAlgorithm;
using DSIType = vtkm::filter::particleadvection::DataSetIntegrator;
this->ValidateOptions();
vtkm::filter::particleadvection::BoundsMap boundsMap(input);
auto dsi = this->CreateDataSetIntegrators(input, boundsMap);
if (this->GetUseThreadedAlgorithm())
return vtkm::filter::particleadvection::RunAlgo<DSIType, ThreadedAlgorithmType>(
boundsMap, dsi, this->NumberOfSteps, this->StepSize, this->Seeds);
else
return vtkm::filter::particleadvection::RunAlgo<DSIType, AlgorithmType>(
boundsMap, dsi, this->NumberOfSteps, this->StepSize, this->Seeds);
}
}
} // namespace vtkm::filter
#endif

@ -1,51 +0,0 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
//
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#ifndef vtk_m_filter_PathParticle_h
#define vtk_m_filter_PathParticle_h
#include <vtkm/filter/FilterTemporalParticleAdvection.h>
namespace vtkm
{
namespace filter
{
/// \brief Advect particles in a time varying vector field.
/// Takes as input a vector field and seed locations and generates the
/// paths taken by the seeds through the vector field.
template <typename ParticleType>
class PathParticleBase
: public vtkm::filter::FilterTemporalParticleAdvection<PathParticleBase<ParticleType>,
ParticleType>
{
public:
VTKM_CONT
PathParticleBase();
template <typename DerivedPolicy>
vtkm::cont::PartitionedDataSet PrepareForExecution(
const vtkm::cont::PartitionedDataSet& input,
const vtkm::filter::PolicyBase<DerivedPolicy>& policy);
protected:
private:
};
using PathParticle = PathParticleBase<vtkm::Particle>;
}
} // namespace vtkm::filter
#ifndef vtk_m_filter_PathParticle_hxx
#include <vtkm/filter/PathParticle.hxx>
#endif
#endif // vtk_m_filter_PathParticle_h

@ -1,58 +0,0 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
//
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#ifndef vtk_m_filter_PathParticle_hxx
#define vtk_m_filter_PathParticle_hxx
#include <vtkm/cont/ErrorFilterExecution.h>
#include <vtkm/filter/PathParticle.h>
#include <vtkm/filter/particleadvection/BoundsMap.h>
#include <vtkm/filter/particleadvection/DataSetIntegrator.h>
#include <vtkm/filter/particleadvection/ParticleAdvectionAlgorithm.h>
namespace vtkm
{
namespace filter
{
//-----------------------------------------------------------------------------
template <typename ParticleType>
inline VTKM_CONT PathParticleBase<ParticleType>::PathParticleBase()
: vtkm::filter::FilterTemporalParticleAdvection<PathParticleBase<ParticleType>, ParticleType>()
{
}
//-----------------------------------------------------------------------------
template <typename ParticleType>
template <typename DerivedPolicy>
inline VTKM_CONT vtkm::cont::PartitionedDataSet PathParticleBase<ParticleType>::PrepareForExecution(
const vtkm::cont::PartitionedDataSet& input,
const vtkm::filter::PolicyBase<DerivedPolicy>&)
{
using AlgorithmType = vtkm::filter::particleadvection::PathParticleAlgorithm;
using ThreadedAlgorithmType = vtkm::filter::particleadvection::PathParticleThreadedAlgorithm;
using TDSIType = vtkm::filter::particleadvection::TemporalDataSetIntegrator;
this->ValidateOptions(input);
vtkm::filter::particleadvection::BoundsMap boundsMap(input);
auto dsi = this->CreateDataSetIntegrators(input, boundsMap);
if (this->GetUseThreadedAlgorithm())
return vtkm::filter::particleadvection::RunAlgo<TDSIType, ThreadedAlgorithmType>(
boundsMap, dsi, this->NumberOfSteps, this->StepSize, this->Seeds);
else
return vtkm::filter::particleadvection::RunAlgo<TDSIType, AlgorithmType>(
boundsMap, dsi, this->NumberOfSteps, this->StepSize, this->Seeds);
}
}
} // namespace vtkm::filter
#endif

@ -1,50 +0,0 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
//
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#ifndef vtk_m_filter_Pathline_h
#define vtk_m_filter_Pathline_h
#include <vtkm/filter/FilterTemporalParticleAdvection.h>
namespace vtkm
{
namespace filter
{
/// \brief generate pathlines from a time sequence of vector fields.
/// Takes as input a vector field and seed locations and generates the
/// paths taken by the seeds through the vector field.
template <typename ParticleType>
class PathlineBase
: public vtkm::filter::FilterTemporalParticleAdvection<PathlineBase<ParticleType>, ParticleType>
{
public:
VTKM_CONT
PathlineBase();
template <typename DerivedPolicy>
vtkm::cont::PartitionedDataSet PrepareForExecution(
const vtkm::cont::PartitionedDataSet& input,
const vtkm::filter::PolicyBase<DerivedPolicy>& policy);
protected:
private:
};
using Pathline = PathlineBase<vtkm::Particle>;
}
} // namespace vtkm::filter
#ifndef vtk_m_filter_Pathline_hxx
#include <vtkm/filter/Pathline.hxx>
#endif
#endif // vtk_m_filter_Pathline_h

@ -1,58 +0,0 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
//
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#ifndef vtk_m_filter_Pathline_hxx
#define vtk_m_filter_Pathline_hxx
#include <vtkm/cont/ErrorFilterExecution.h>
#include <vtkm/filter/Pathline.h>
#include <vtkm/filter/particleadvection/BoundsMap.h>
#include <vtkm/filter/particleadvection/DataSetIntegrator.h>
#include <vtkm/filter/particleadvection/ParticleAdvectionAlgorithm.h>
namespace vtkm
{
namespace filter
{
//-----------------------------------------------------------------------------
template <typename ParticleType>
inline VTKM_CONT PathlineBase<ParticleType>::PathlineBase()
: vtkm::filter::FilterTemporalParticleAdvection<PathlineBase<ParticleType>, ParticleType>()
{
}
//-----------------------------------------------------------------------------
template <typename ParticleType>
template <typename DerivedPolicy>
inline VTKM_CONT vtkm::cont::PartitionedDataSet PathlineBase<ParticleType>::PrepareForExecution(
const vtkm::cont::PartitionedDataSet& input,
const vtkm::filter::PolicyBase<DerivedPolicy>&)
{
using AlgorithmType = vtkm::filter::particleadvection::PathlineAlgorithm;
using ThreadedAlgorithmType = vtkm::filter::particleadvection::PathlineThreadedAlgorithm;
using TDSIType = vtkm::filter::particleadvection::TemporalDataSetIntegrator;
this->ValidateOptions(input);
vtkm::filter::particleadvection::BoundsMap boundsMap(input);
auto dsi = this->CreateDataSetIntegrators(input, boundsMap);
if (this->GetUseThreadedAlgorithm())
return vtkm::filter::particleadvection::RunAlgo<TDSIType, ThreadedAlgorithmType>(
boundsMap, dsi, this->NumberOfSteps, this->StepSize, this->Seeds);
else
return vtkm::filter::particleadvection::RunAlgo<TDSIType, AlgorithmType>(
boundsMap, dsi, this->NumberOfSteps, this->StepSize, this->Seeds);
}
}
} // namespace vtkm::filter
#endif

@ -1,70 +0,0 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
//
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#ifndef vtk_m_filter_StreamSurface_h
#define vtk_m_filter_StreamSurface_h
#include <vtkm/filter/FilterDataSetWithField.h>
#include <vtkm/worklet/ParticleAdvection.h>
#include <vtkm/worklet/StreamSurface.h>
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
#include <vtkm/worklet/particleadvection/Stepper.h>
namespace vtkm
{
namespace filter
{
/// \brief generate streamlines from a vector field.
/// Takes as input a vector field and seed locations and generates the
/// paths taken by the seeds through the vector field.
class StreamSurface : public vtkm::filter::FilterDataSetWithField<StreamSurface>
{
public:
using SupportedTypes = vtkm::TypeListFieldVec3;
VTKM_CONT
StreamSurface();
VTKM_CONT
void SetStepSize(vtkm::FloatDefault s) { this->StepSize = s; }
VTKM_CONT
void SetNumberOfSteps(vtkm::Id n) { this->NumberOfSteps = n; }
VTKM_CONT
void SetSeeds(vtkm::cont::ArrayHandle<vtkm::Particle>& seeds) { this->Seeds = seeds; }
template <typename T, typename StorageType, typename DerivedPolicy>
VTKM_CONT vtkm::cont::DataSet DoExecute(
const vtkm::cont::DataSet& input,
const vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>, StorageType>& field,
const vtkm::filter::FieldMetadata& fieldMeta,
const vtkm::filter::PolicyBase<DerivedPolicy>& policy);
template <typename DerivedPolicy>
VTKM_CONT bool MapFieldOntoOutput(vtkm::cont::DataSet& result,
const vtkm::cont::Field& field,
vtkm::filter::PolicyBase<DerivedPolicy> policy);
private:
vtkm::Id NumberOfSteps;
vtkm::cont::ArrayHandle<vtkm::Particle> Seeds;
vtkm::FloatDefault StepSize;
vtkm::worklet::StreamSurface Worklet;
};
}
} // namespace vtkm::filter
#ifndef vtk_m_filter_StreamSurface_hxx
#include <vtkm/filter/StreamSurface.hxx>
#endif
#endif // vtk_m_filter_StreamSurface_h

@ -1,95 +0,0 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
//
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#ifndef vtk_m_filter_StreamSurface_hxx
#define vtk_m_filter_StreamSurface_hxx
#include <vtkm/filter/StreamSurface.h>
#include <vtkm/cont/ArrayCopy.h>
#include <vtkm/cont/ArrayHandleIndex.h>
#include <vtkm/cont/ErrorFilterExecution.h>
#include <vtkm/worklet/ParticleAdvection.h>
#include <vtkm/worklet/particleadvection/Field.h>
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
#include <vtkm/worklet/particleadvection/Particles.h>
#include <vtkm/worklet/particleadvection/RK4Integrator.h>
#include <vtkm/worklet/particleadvection/Stepper.h>
namespace vtkm
{
namespace filter
{
//-----------------------------------------------------------------------------
inline VTKM_CONT StreamSurface::StreamSurface()
: vtkm::filter::FilterDataSetWithField<StreamSurface>()
, Worklet()
{
}
//-----------------------------------------------------------------------------
template <typename T, typename StorageType, typename DerivedPolicy>
inline VTKM_CONT vtkm::cont::DataSet StreamSurface::DoExecute(
const vtkm::cont::DataSet& input,
const vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>, StorageType>& field,
const vtkm::filter::FieldMetadata& fieldMeta,
const vtkm::filter::PolicyBase<DerivedPolicy>&)
{
//Check for some basics.
if (this->Seeds.GetNumberOfValues() == 0)
throw vtkm::cont::ErrorFilterExecution("No seeds provided.");
const vtkm::cont::UnknownCellSet& cells = input.GetCellSet();
const vtkm::cont::CoordinateSystem& coords =
input.GetCoordinateSystem(this->GetActiveCoordinateSystemIndex());
using FieldHandle = vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>, StorageType>;
using FieldType = vtkm::worklet::particleadvection::VelocityField<FieldHandle>;
using GridEvalType = vtkm::worklet::particleadvection::GridEvaluator<FieldType>;
using RK4Type = vtkm::worklet::particleadvection::RK4Integrator<GridEvalType>;
using Stepper = vtkm::worklet::particleadvection::Stepper<RK4Type, GridEvalType>;
//compute streamlines
FieldType velocities(field, fieldMeta.GetAssociation());
GridEvalType eval(coords, cells, velocities);
Stepper rk4(eval, this->StepSize);
vtkm::worklet::Streamline streamline;
vtkm::cont::ArrayHandle<vtkm::Particle> seedArray;
vtkm::cont::ArrayCopy(this->Seeds, seedArray);
auto res = streamline.Run(rk4, seedArray, this->NumberOfSteps);
//compute surface from streamlines
vtkm::cont::ArrayHandle<vtkm::Vec3f> srfPoints;
vtkm::cont::CellSetSingleType<> srfCells;
vtkm::cont::CoordinateSystem slCoords("coordinates", res.Positions);
this->Worklet.Run(slCoords, res.PolyLines, srfPoints, srfCells);
vtkm::cont::DataSet outData;
vtkm::cont::CoordinateSystem outputCoords("coordinates", srfPoints);
outData.AddCoordinateSystem(outputCoords);
outData.SetCellSet(srfCells);
return outData;
}
//-----------------------------------------------------------------------------
template <typename DerivedPolicy>
inline VTKM_CONT bool StreamSurface::MapFieldOntoOutput(vtkm::cont::DataSet&,
const vtkm::cont::Field&,
vtkm::filter::PolicyBase<DerivedPolicy>)
{
return false;
}
}
} // namespace vtkm::filter
#endif

@ -1,48 +0,0 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
//
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#ifndef vtk_m_filter_Streamline_h
#define vtk_m_filter_Streamline_h
#include <vtkm/filter/FilterParticleAdvection.h>
namespace vtkm
{
namespace filter
{
/// \brief Generate streamlines from a vector field.
/// Takes as input a vector field and seed locations and generates the
/// paths taken by the seeds through the vector field.
template <typename ParticleType = vtkm::Particle>
class StreamlineBase
: public vtkm::filter::FilterParticleAdvection<StreamlineBase<ParticleType>, ParticleType>
{
public:
VTKM_CONT
StreamlineBase();
template <typename DerivedPolicy>
vtkm::cont::PartitionedDataSet PrepareForExecution(
const vtkm::cont::PartitionedDataSet& input,
const vtkm::filter::PolicyBase<DerivedPolicy>& policy);
private:
};
using Streamline = StreamlineBase<vtkm::Particle>;
}
} // namespace vtkm::filter
#ifndef vtk_m_filter_Streamline_hxx
#include <vtkm/filter/Streamline.hxx>
#endif
#endif // vtk_m_filter_Streamline_h

@ -1,56 +0,0 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
//
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#ifndef vtk_m_filter_Streamline_hxx
#define vtk_m_filter_Streamline_hxx
#include <vtkm/cont/ErrorFilterExecution.h>
#include <vtkm/filter/Streamline.h>
#include <vtkm/filter/particleadvection/BoundsMap.h>
#include <vtkm/filter/particleadvection/DataSetIntegrator.h>
#include <vtkm/filter/particleadvection/ParticleAdvectionAlgorithm.h>
namespace vtkm
{
namespace filter
{
//-----------------------------------------------------------------------------
template <typename ParticleType>
inline VTKM_CONT StreamlineBase<ParticleType>::StreamlineBase()
: vtkm::filter::FilterParticleAdvection<StreamlineBase<ParticleType>, ParticleType>()
{
}
//-----------------------------------------------------------------------------
template <typename ParticleType>
template <typename DerivedPolicy>
inline VTKM_CONT vtkm::cont::PartitionedDataSet StreamlineBase<ParticleType>::PrepareForExecution(
const vtkm::cont::PartitionedDataSet& input,
const vtkm::filter::PolicyBase<DerivedPolicy>&)
{
using AlgorithmType = vtkm::filter::particleadvection::StreamlineAlgorithm;
using ThreadedAlgorithmType = vtkm::filter::particleadvection::StreamlineThreadedAlgorithm;
using DSIType = vtkm::filter::particleadvection::DataSetIntegrator;
this->ValidateOptions();
vtkm::filter::particleadvection::BoundsMap boundsMap(input);
auto dsi = this->CreateDataSetIntegrators(input, boundsMap);
if (this->GetUseThreadedAlgorithm())
return vtkm::filter::particleadvection::RunAlgo<DSIType, ThreadedAlgorithmType>(
boundsMap, dsi, this->NumberOfSteps, this->StepSize, this->Seeds);
else
return vtkm::filter::particleadvection::RunAlgo<DSIType, AlgorithmType>(
boundsMap, dsi, this->NumberOfSteps, this->StepSize, this->Seeds);
}
}
} // namespace vtkm::filter
#endif

@ -308,8 +308,9 @@ namespace entity_extraction
//-----------------------------------------------------------------------------
VTKM_CONT GhostCellRemove::GhostCellRemove()
{
this->SetActiveField("vtkmGhostCells");
this->SetFieldsToPass("vtkmGhostCells", vtkm::filter::FieldSelection::Mode::Exclude);
this->SetActiveField(vtkm::cont::GetGlobalGhostCellFieldName());
this->SetFieldsToPass(vtkm::cont::GetGlobalGhostCellFieldName(),
vtkm::filter::FieldSelection::Mode::Exclude);
}
//-----------------------------------------------------------------------------

@ -94,7 +94,7 @@ vtkm::cont::DataSet MakeUniform(vtkm::Id numI,
ds = vtkm::cont::DataSetBuilderUniform::Create(vtkm::Id3(numI + 1, numJ + 1, numK + 1));
auto ghosts = StructuredGhostCellArray(numI, numJ, numK, numLayers, addMidGhost);
ds.AddCellField("vtkmGhostCells", ghosts);
ds.AddGhostCellField(ghosts);
return ds;
}
@ -128,7 +128,7 @@ vtkm::cont::DataSet MakeRectilinear(vtkm::Id numI,
auto ghosts = StructuredGhostCellArray(numI, numJ, numK, numLayers, addMidGhost);
ds.AddCellField("vtkmGhostCells", ghosts);
ds.AddGhostCellField(ghosts);
return ds;
}
@ -205,7 +205,7 @@ vtkm::cont::DataSet MakeExplicit(vtkm::Id numI, vtkm::Id numJ, vtkm::Id numK, in
auto ghosts = StructuredGhostCellArray(numI, numJ, numK, numLayers);
ds.AddCellField("vtkmGhostCells", ghosts);
ds.AddGhostCellField(ghosts);
return ds;
}

@ -0,0 +1,57 @@
##============================================================================
## Copyright (c) Kitware, Inc.
## All rights reserved.
## See LICENSE.txt for details.
##
## This software is distributed WITHOUT ANY WARRANTY; without even
## the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
## PURPOSE. See the above copyright notice for more information.
##============================================================================
set(flow_headers
FlowTypes.h
NewFilterParticleAdvection.h
NewFilterParticleAdvectionSteadyState.h
NewFilterParticleAdvectionUnsteadyState.h
ParticleAdvection.h
Pathline.h
PathParticle.h
Streamline.h
StreamSurface.h
)
set(flow_sources
internal/Messenger.cxx
NewFilterParticleAdvection.cxx
ParticleAdvection.cxx
Pathline.cxx
PathParticle.cxx
Streamline.cxx
)
set(flow_device_sources
NewFilterParticleAdvectionSteadyState.cxx
NewFilterParticleAdvectionUnsteadyState.cxx
StreamSurface.cxx
)
vtkm_library(
NAME vtkm_filter_flow
HEADERS ${flow_headers}
SOURCES ${flow_sources}
DEVICE_SOURCES ${flow_device_sources}
USE_VTKM_JOB_POOL
)
target_link_libraries(vtkm_filter_flow PUBLIC vtkm_worklet vtkm_filter_core)
target_link_libraries(vtkm_filter PUBLIC INTERFACE vtkm_filter_flow)
if (VTKm_ENABLE_MPI)
target_link_libraries(vtkm_filter_flow PUBLIC MPI::MPI_CXX)
endif()
add_subdirectory(internal)
add_subdirectory(worklet)
#-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -
if (VTKm_ENABLE_TESTING)
add_subdirectory(testing)
endif ()

@ -0,0 +1,42 @@
//============================================================================
// 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_flow_FlowTypes_h
#define vtk_m_filter_flow_FlowTypes_h
namespace vtkm
{
namespace filter
{
namespace flow
{
enum class IntegrationSolverType
{
RK4_TYPE = 0,
EULER_TYPE,
};
enum class VectorFieldType
{
VELOCITY_FIELD_TYPE = 0,
ELECTRO_MAGNETIC_FIELD_TYPE,
};
enum FlowResultType
{
UNKNOWN_TYPE = -1,
PARTICLE_ADVECT_TYPE,
STREAMLINE_TYPE,
};
}
}
}
#endif // vtk_m_filter_flow_FlowTypes_h

@ -0,0 +1,54 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
//
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#include <vtkm/Particle.h>
#include <vtkm/cont/ErrorFilterExecution.h>
#include <vtkm/filter/flow/NewFilterParticleAdvection.h>
namespace vtkm
{
namespace filter
{
namespace flow
{
VTKM_CONT vtkm::cont::DataSet NewFilterParticleAdvection::DoExecute(
const vtkm::cont::DataSet& inData)
{
auto out = this->DoExecutePartitions(inData);
if (out.GetNumberOfPartitions() != 1)
throw vtkm::cont::ErrorFilterExecution("Wrong number of results");
return out.GetPartition(0);
}
VTKM_CONT void NewFilterParticleAdvection::ValidateOptions() const
{
if (this->GetUseCoordinateSystemAsField())
throw vtkm::cont::ErrorFilterExecution("Coordinate system as field not supported");
if (this->Seeds.GetNumberOfValues() == 0)
throw vtkm::cont::ErrorFilterExecution("No seeds provided.");
if (!this->Seeds.IsBaseComponentType<vtkm::Particle>() &&
this->Seeds.IsBaseComponentType<vtkm::ChargedParticle>())
throw vtkm::cont::ErrorFilterExecution("Unsupported particle type in seed array.");
if (this->NumberOfSteps == 0)
throw vtkm::cont::ErrorFilterExecution("Number of steps not specified.");
if (this->StepSize == 0)
throw vtkm::cont::ErrorFilterExecution("Step size not specified.");
if (this->NumberOfSteps < 0)
throw vtkm::cont::ErrorFilterExecution("NumberOfSteps cannot be negative");
if (this->StepSize < 0)
throw vtkm::cont::ErrorFilterExecution("StepSize cannot be negative");
}
}
}
} // namespace vtkm::filter::flow

@ -0,0 +1,88 @@
//============================================================================
// 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_flow_NewFilterParticleAdvection_h
#define vtk_m_filter_flow_NewFilterParticleAdvection_h
#include <vtkm/filter/NewFilterField.h>
#include <vtkm/filter/flow/FlowTypes.h>
#include <vtkm/filter/flow/vtkm_filter_flow_export.h>
namespace vtkm
{
namespace filter
{
namespace flow
{
/// \brief base class for advecting particles in a vector field.
/// Takes as input a vector field and seed locations and advects the seeds
/// through the flow field.
class VTKM_FILTER_FLOW_EXPORT NewFilterParticleAdvection : public vtkm::filter::NewFilterField
{
public:
VTKM_CONT
void SetStepSize(vtkm::FloatDefault s) { this->StepSize = s; }
VTKM_CONT
void SetNumberOfSteps(vtkm::Id n) { this->NumberOfSteps = n; }
template <typename ParticleType>
VTKM_CONT void SetSeeds(vtkm::cont::ArrayHandle<ParticleType>& seeds)
{
this->Seeds = seeds;
}
template <typename ParticleType>
VTKM_CONT void SetSeeds(const std::vector<ParticleType>& seeds,
vtkm::CopyFlag copyFlag = vtkm::CopyFlag::On)
{
this->Seeds = vtkm::cont::make_ArrayHandle(seeds, copyFlag);
}
VTKM_CONT
void SetSolverRK4() { this->SolverType = vtkm::filter::flow::IntegrationSolverType::RK4_TYPE; }
VTKM_CONT
void SetSolverEuler()
{
this->SolverType = vtkm::filter::flow::IntegrationSolverType::EULER_TYPE;
}
VTKM_CONT
bool GetUseThreadedAlgorithm() { return this->UseThreadedAlgorithm; }
VTKM_CONT
void SetUseThreadedAlgorithm(bool val) { this->UseThreadedAlgorithm = val; }
protected:
VTKM_CONT virtual void ValidateOptions() const;
VTKM_CONT virtual vtkm::filter::flow::FlowResultType GetResultType() const = 0;
vtkm::Id NumberOfSteps = 0;
vtkm::cont::UnknownArrayHandle Seeds;
vtkm::filter::flow::IntegrationSolverType SolverType =
vtkm::filter::flow::IntegrationSolverType::RK4_TYPE;
vtkm::FloatDefault StepSize = 0;
bool UseThreadedAlgorithm = false;
vtkm::filter::flow::VectorFieldType VecFieldType =
vtkm::filter::flow::VectorFieldType::VELOCITY_FIELD_TYPE;
private:
VTKM_CONT vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet& inData) override;
};
}
}
} // namespace vtkm::filter::flow
#endif // vtk_m_filter_flow_NewFilterParticleAdvection_h

@ -0,0 +1,73 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
//
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#include <vtkm/filter/flow/NewFilterParticleAdvectionSteadyState.h>
#include <vtkm/filter/flow/internal/DataSetIntegratorSteadyState.h>
#include <vtkm/filter/flow/internal/BoundsMap.h>
#include <vtkm/filter/flow/internal/ParticleAdvector.h>
namespace vtkm
{
namespace filter
{
namespace flow
{
namespace
{
VTKM_CONT std::vector<vtkm::filter::flow::internal::DataSetIntegratorSteadyState>
CreateDataSetIntegrators(const vtkm::cont::PartitionedDataSet& input,
const std::string& activeField,
const vtkm::filter::flow::internal::BoundsMap& boundsMap,
const vtkm::filter::flow::IntegrationSolverType solverType,
const vtkm::filter::flow::VectorFieldType vecFieldType,
const vtkm::filter::flow::FlowResultType resultType)
{
using DSIType = vtkm::filter::flow::internal::DataSetIntegratorSteadyState;
std::vector<DSIType> dsi;
for (vtkm::Id i = 0; i < input.GetNumberOfPartitions(); i++)
{
vtkm::Id blockId = boundsMap.GetLocalBlockId(i);
auto ds = input.GetPartition(i);
if (!ds.HasPointField(activeField) && !ds.HasCellField(activeField))
throw vtkm::cont::ErrorFilterExecution("Unsupported field assocation");
dsi.emplace_back(ds, blockId, activeField, solverType, vecFieldType, resultType);
}
return dsi;
}
} //anonymous namespace
VTKM_CONT vtkm::cont::PartitionedDataSet NewFilterParticleAdvectionSteadyState::DoExecutePartitions(
const vtkm::cont::PartitionedDataSet& input)
{
using DSIType = vtkm::filter::flow::internal::DataSetIntegratorSteadyState;
this->ValidateOptions();
vtkm::filter::flow::internal::BoundsMap boundsMap(input);
auto dsi = CreateDataSetIntegrators(input,
this->GetActiveFieldName(),
boundsMap,
this->SolverType,
this->VecFieldType,
this->GetResultType());
vtkm::filter::flow::internal::ParticleAdvector<DSIType> pav(
boundsMap, dsi, this->UseThreadedAlgorithm, this->GetResultType());
return pav.Execute(this->NumberOfSteps, this->StepSize, this->Seeds);
}
}
}
} // namespace vtkm::filter::flow

@ -0,0 +1,36 @@
//============================================================================
// 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_flow_NewFilterParticleAdvectionSteadyState_h
#define vtk_m_filter_flow_NewFilterParticleAdvectionSteadyState_h
#include <vtkm/filter/flow/FlowTypes.h>
#include <vtkm/filter/flow/NewFilterParticleAdvection.h>
#include <vtkm/filter/flow/vtkm_filter_flow_export.h>
namespace vtkm
{
namespace filter
{
namespace flow
{
class VTKM_FILTER_FLOW_EXPORT NewFilterParticleAdvectionSteadyState
: public NewFilterParticleAdvection
{
private:
VTKM_CONT vtkm::cont::PartitionedDataSet DoExecutePartitions(
const vtkm::cont::PartitionedDataSet& inData) override;
};
}
}
} // namespace vtkm::filter::flow
#endif // vtk_m_filter_flow_NewFilterParticleAdvectionSteadyState_h

@ -0,0 +1,87 @@
//============================================================================
// 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/filter/flow/NewFilterParticleAdvectionUnsteadyState.h>
#include <vtkm/filter/flow/internal/DataSetIntegratorUnsteadyState.h>
#include <vtkm/filter/flow/internal/ParticleAdvector.h>
namespace vtkm
{
namespace filter
{
namespace flow
{
namespace
{
VTKM_CONT std::vector<vtkm::filter::flow::internal::DataSetIntegratorUnsteadyState>
CreateDataSetIntegrators(const vtkm::cont::PartitionedDataSet& input,
const vtkm::cont::PartitionedDataSet& input2,
const std::string& activeField,
const vtkm::FloatDefault timer1,
const vtkm::FloatDefault timer2,
const vtkm::filter::flow::internal::BoundsMap& boundsMap,
const vtkm::filter::flow::IntegrationSolverType solverType,
const vtkm::filter::flow::VectorFieldType vecFieldType,
const vtkm::filter::flow::FlowResultType resultType)
{
using DSIType = vtkm::filter::flow::internal::DataSetIntegratorUnsteadyState;
std::vector<DSIType> dsi;
for (vtkm::Id i = 0; i < input.GetNumberOfPartitions(); i++)
{
vtkm::Id blockId = boundsMap.GetLocalBlockId(i);
auto ds1 = input.GetPartition(i);
auto ds2 = input2.GetPartition(i);
if ((!ds1.HasPointField(activeField) && !ds1.HasCellField(activeField)) ||
(!ds2.HasPointField(activeField) && !ds2.HasCellField(activeField)))
throw vtkm::cont::ErrorFilterExecution("Unsupported field assocation");
dsi.emplace_back(
ds1, ds2, timer1, timer2, blockId, activeField, solverType, vecFieldType, resultType);
}
return dsi;
}
} // anonymous namespace
VTKM_CONT void NewFilterParticleAdvectionUnsteadyState::ValidateOptions() const
{
this->NewFilterParticleAdvection::ValidateOptions();
if (this->Time1 >= this->Time2)
throw vtkm::cont::ErrorFilterExecution("PreviousTime must be less than NextTime");
}
VTKM_CONT vtkm::cont::PartitionedDataSet
NewFilterParticleAdvectionUnsteadyState::DoExecutePartitions(
const vtkm::cont::PartitionedDataSet& input)
{
using DSIType = vtkm::filter::flow::internal::DataSetIntegratorUnsteadyState;
this->ValidateOptions();
vtkm::filter::flow::internal::BoundsMap boundsMap(input);
auto dsi = CreateDataSetIntegrators(input,
this->Input2,
this->GetActiveFieldName(),
this->Time1,
this->Time2,
boundsMap,
this->SolverType,
this->VecFieldType,
this->GetResultType());
vtkm::filter::flow::internal::ParticleAdvector<DSIType> pav(
boundsMap, dsi, this->UseThreadedAlgorithm, this->GetResultType());
return pav.Execute(this->NumberOfSteps, this->StepSize, this->Seeds);
}
}
}
} // namespace vtkm::filter::flow

@ -0,0 +1,52 @@
//============================================================================
// 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_flow_NewFilterParticleAdvectionUnsteadyState_h
#define vtk_m_filter_flow_NewFilterParticleAdvectionUnsteadyState_h
#include <vtkm/filter/flow/NewFilterParticleAdvection.h>
#include <vtkm/filter/flow/vtkm_filter_flow_export.h>
namespace vtkm
{
namespace filter
{
namespace flow
{
class VTKM_FILTER_FLOW_EXPORT NewFilterParticleAdvectionUnsteadyState
: public NewFilterParticleAdvection
{
public:
VTKM_CONT void SetPreviousTime(vtkm::FloatDefault t1) { this->Time1 = t1; }
VTKM_CONT void SetNextTime(vtkm::FloatDefault t2) { this->Time2 = t2; }
VTKM_CONT void SetNextDataSet(const vtkm::cont::DataSet& ds) { this->Input2 = { ds }; }
VTKM_CONT void SetNextDataSet(const vtkm::cont::PartitionedDataSet& pds) { this->Input2 = pds; }
protected:
VTKM_CONT virtual void ValidateOptions() const override;
private:
VTKM_CONT vtkm::cont::PartitionedDataSet DoExecutePartitions(
const vtkm::cont::PartitionedDataSet& inData) override;
vtkm::cont::PartitionedDataSet Input2;
vtkm::FloatDefault Time1 = -1;
vtkm::FloatDefault Time2 = -1;
};
}
}
} // namespace vtkm::filter::flow
#endif // vtk_m_filter_flow_NewFilterParticleAdvectionUnsteadyState_h

@ -0,0 +1,27 @@
//============================================================================
// 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/filter/flow/ParticleAdvection.h>
namespace vtkm
{
namespace filter
{
namespace flow
{
VTKM_CONT vtkm::filter::flow::FlowResultType ParticleAdvection::GetResultType() const
{
return vtkm::filter::flow::FlowResultType::PARTICLE_ADVECT_TYPE;
}
}
}
} // namespace vtkm::filter::flow

@ -0,0 +1,42 @@
//============================================================================
// 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_flow_ParticleAdvection_h
#define vtk_m_filter_flow_ParticleAdvection_h
#include <vtkm/filter/flow/FlowTypes.h>
#include <vtkm/filter/flow/NewFilterParticleAdvectionSteadyState.h>
#include <vtkm/filter/flow/vtkm_filter_flow_export.h>
namespace vtkm
{
namespace filter
{
namespace flow
{
/// \brief Advect particles in a vector field.
/// Takes as input a vector field and seed locations and generates the
/// end points for each seed through the vector field.
class VTKM_FILTER_FLOW_EXPORT ParticleAdvection
: public vtkm::filter::flow::NewFilterParticleAdvectionSteadyState
{
public:
private:
VTKM_CONT vtkm::filter::flow::FlowResultType GetResultType() const override;
};
}
}
} // namespace vtkm::filter::flow
#endif // vtk_m_filter_flow_ParticleAdvection_h

@ -0,0 +1,27 @@
//============================================================================
// 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/filter/flow/PathParticle.h>
namespace vtkm
{
namespace filter
{
namespace flow
{
VTKM_CONT vtkm::filter::flow::FlowResultType PathParticle::GetResultType() const
{
return vtkm::filter::flow::FlowResultType::PARTICLE_ADVECT_TYPE;
}
}
}
} // namespace vtkm::filter::flow

@ -0,0 +1,41 @@
//============================================================================
// 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_PathParticle_h
#define vtk_m_filter_PathParticle_h
#include <vtkm/filter/flow/FlowTypes.h>
#include <vtkm/filter/flow/NewFilterParticleAdvectionUnsteadyState.h>
#include <vtkm/filter/flow/vtkm_filter_flow_export.h>
namespace vtkm
{
namespace filter
{
namespace flow
{
/// \brief Advect particles in a vector field.
/// Takes as input a vector field and seed locations and generates the
/// end points for each seed through the vector field.
class VTKM_FILTER_FLOW_EXPORT PathParticle
: public vtkm::filter::flow::NewFilterParticleAdvectionUnsteadyState
{
private:
VTKM_CONT vtkm::filter::flow::FlowResultType GetResultType() const override;
};
}
}
} // namespace vtkm::filter::flow
#endif // vtk_m_filter_PathParticle_h

@ -0,0 +1,27 @@
//============================================================================
// 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/filter/flow/Pathline.h>
namespace vtkm
{
namespace filter
{
namespace flow
{
VTKM_CONT vtkm::filter::flow::FlowResultType Pathline::GetResultType() const
{
return vtkm::filter::flow::FlowResultType::STREAMLINE_TYPE;
}
}
}
} // namespace vtkm::filter::flow

@ -0,0 +1,41 @@
//============================================================================
// 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_Pathline_h
#define vtk_m_filter_Pathline_h
#include <vtkm/filter/flow/FlowTypes.h>
#include <vtkm/filter/flow/NewFilterParticleAdvectionUnsteadyState.h>
#include <vtkm/filter/flow/vtkm_filter_flow_export.h>
namespace vtkm
{
namespace filter
{
namespace flow
{
/// \brief Advect particles in a vector field.
/// Takes as input a vector field and seed locations and generates the
/// end points for each seed through the vector field.
class VTKM_FILTER_FLOW_EXPORT Pathline
: public vtkm::filter::flow::NewFilterParticleAdvectionUnsteadyState
{
private:
VTKM_CONT vtkm::filter::flow::FlowResultType GetResultType() const override;
};
}
}
} // namespace vtkm::filter::flow
#endif // vtk_m_filter_Pathline_h

@ -0,0 +1,92 @@
//============================================================================
// 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/ErrorFilterExecution.h>
#include <vtkm/filter/flow/StreamSurface.h>
#include <vtkm/filter/flow/worklet/Field.h>
#include <vtkm/filter/flow/worklet/GridEvaluators.h>
#include <vtkm/filter/flow/worklet/ParticleAdvection.h>
#include <vtkm/filter/flow/worklet/RK4Integrator.h>
#include <vtkm/filter/flow/worklet/Stepper.h>
#include <vtkm/filter/flow/worklet/StreamSurface.h>
namespace vtkm
{
namespace filter
{
namespace flow
{
VTKM_CONT vtkm::cont::DataSet StreamSurface::DoExecute(const vtkm::cont::DataSet& input)
{
//Validate inputs.
if (this->GetUseCoordinateSystemAsField())
throw vtkm::cont::ErrorFilterExecution("Coordinate system as field not supported");
if (this->Seeds.GetNumberOfValues() == 0)
throw vtkm::cont::ErrorFilterExecution("No seeds provided.");
if (!this->Seeds.IsBaseComponentType<vtkm::Particle>() &&
this->Seeds.IsBaseComponentType<vtkm::ChargedParticle>())
throw vtkm::cont::ErrorFilterExecution("Unsupported particle type in seed array.");
if (this->NumberOfSteps == 0)
throw vtkm::cont::ErrorFilterExecution("Number of steps not specified.");
if (this->StepSize == 0)
throw vtkm::cont::ErrorFilterExecution("Step size not specified.");
if (this->NumberOfSteps < 0)
throw vtkm::cont::ErrorFilterExecution("NumberOfSteps cannot be negative");
if (this->StepSize < 0)
throw vtkm::cont::ErrorFilterExecution("StepSize cannot be negative");
if (!this->Seeds.IsBaseComponentType<vtkm::Particle>())
throw vtkm::cont::ErrorFilterExecution("Unsupported seed type in StreamSurface filter.");
const vtkm::cont::UnknownCellSet& cells = input.GetCellSet();
const vtkm::cont::CoordinateSystem& coords =
input.GetCoordinateSystem(this->GetActiveCoordinateSystemIndex());
using FieldHandle = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
using FieldType = vtkm::worklet::flow::VelocityField<FieldHandle>;
using GridEvalType = vtkm::worklet::flow::GridEvaluator<FieldType>;
using RK4Type = vtkm::worklet::flow::RK4Integrator<GridEvalType>;
using Stepper = vtkm::worklet::flow::Stepper<RK4Type, GridEvalType>;
//compute streamlines
const auto& field = input.GetField(this->GetActiveFieldName());
FieldHandle arr;
vtkm::cont::ArrayCopyShallowIfPossible(field.GetData(), arr);
FieldType velocities(arr, field.GetAssociation());
GridEvalType eval(coords, cells, velocities);
Stepper rk4(eval, this->StepSize);
vtkm::worklet::flow::Streamline streamline;
using ParticleArray = vtkm::cont::ArrayHandle<vtkm::Particle>;
vtkm::cont::ArrayHandle<vtkm::Particle> seedArray;
vtkm::cont::ArrayCopy(this->Seeds.AsArrayHandle<ParticleArray>(), seedArray);
auto res = streamline.Run(rk4, seedArray, this->NumberOfSteps);
//compute surface from streamlines
vtkm::worklet::flow::StreamSurface streamSurface;
vtkm::cont::ArrayHandle<vtkm::Vec3f> srfPoints;
vtkm::cont::CellSetSingleType<> srfCells;
vtkm::cont::CoordinateSystem slCoords("coordinates", res.Positions);
streamSurface.Run(slCoords, res.PolyLines, srfPoints, srfCells);
vtkm::cont::DataSet outData;
vtkm::cont::CoordinateSystem outputCoords("coordinates", srfPoints);
outData.AddCoordinateSystem(outputCoords);
outData.SetCellSet(srfCells);
return outData;
}
}
}
} // namespace vtkm::filter::flow

@ -0,0 +1,64 @@
//============================================================================
// 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_flow_StreamSurface_h
#define vtk_m_filter_flow_StreamSurface_h
#include <vtkm/filter/NewFilterField.h>
#include <vtkm/filter/flow/FlowTypes.h>
#include <vtkm/filter/flow/vtkm_filter_flow_export.h>
namespace vtkm
{
namespace filter
{
namespace flow
{
/// \brief generate streamlines from a vector field.
/// Takes as input a vector field and seed locations and generates the
/// paths taken by the seeds through the vector field.
class VTKM_FILTER_FLOW_EXPORT StreamSurface : public vtkm::filter::NewFilterField
{
public:
VTKM_CONT
void SetStepSize(vtkm::FloatDefault s) { this->StepSize = s; }
VTKM_CONT
void SetNumberOfSteps(vtkm::Id n) { this->NumberOfSteps = n; }
template <typename ParticleType>
VTKM_CONT void SetSeeds(vtkm::cont::ArrayHandle<ParticleType>& seeds)
{
this->Seeds = seeds;
}
template <typename ParticleType>
VTKM_CONT void SetSeeds(const std::vector<ParticleType>& seeds,
vtkm::CopyFlag copyFlag = vtkm::CopyFlag::On)
{
this->Seeds = vtkm::cont::make_ArrayHandle(seeds, copyFlag);
}
private:
VTKM_CONT vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet& inData) override;
vtkm::Id NumberOfSteps = 0;
vtkm::cont::UnknownArrayHandle Seeds;
vtkm::FloatDefault StepSize = 0;
};
}
}
} // namespace vtkm::filter::flow
#endif // vtk_m_filter_flow_StreamSurface_h

@ -0,0 +1,27 @@
//============================================================================
// 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/filter/flow/Streamline.h>
namespace vtkm
{
namespace filter
{
namespace flow
{
VTKM_CONT vtkm::filter::flow::FlowResultType Streamline::GetResultType() const
{
return vtkm::filter::flow::FlowResultType::STREAMLINE_TYPE;
}
}
}
} // namespace vtkm::filter::flow

@ -0,0 +1,41 @@
//============================================================================
// 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_flow_Streamline_h
#define vtk_m_filter_flow_Streamline_h
#include <vtkm/filter/flow/FlowTypes.h>
#include <vtkm/filter/flow/NewFilterParticleAdvectionSteadyState.h>
#include <vtkm/filter/flow/vtkm_filter_flow_export.h>
namespace vtkm
{
namespace filter
{
namespace flow
{
/// \brief Advect particles in a vector field.
/// Takes as input a vector field and seed locations and generates the
/// end points for each seed through the vector field.
class VTKM_FILTER_FLOW_EXPORT Streamline
: public vtkm::filter::flow::NewFilterParticleAdvectionSteadyState
{
private:
VTKM_CONT vtkm::filter::flow::FlowResultType GetResultType() const override;
};
}
}
} // namespace vtkm::filter::flow
#endif // vtk_m_filter_flow_Streamline_h

@ -0,0 +1,284 @@
//============================================================================
// 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_flow_internal_AdvectAlgorithm_h
#define vtk_m_filter_flow_internal_AdvectAlgorithm_h
#include <vtkm/cont/PartitionedDataSet.h>
#include <vtkm/filter/flow/internal/BoundsMap.h>
#include <vtkm/filter/flow/internal/DataSetIntegrator.h>
#include <vtkm/filter/flow/internal/ParticleMessenger.h>
namespace vtkm
{
namespace filter
{
namespace flow
{
namespace internal
{
template <typename DSIType, template <typename> class ResultType, typename ParticleType>
class AdvectAlgorithm
{
public:
AdvectAlgorithm(const vtkm::filter::flow::internal::BoundsMap& bm, std::vector<DSIType>& blocks)
: Blocks(blocks)
, BoundsMap(bm)
, NumRanks(this->Comm.size())
, Rank(this->Comm.rank())
{
}
void Execute(vtkm::Id numSteps,
vtkm::FloatDefault stepSize,
const vtkm::cont::ArrayHandle<ParticleType>& seeds)
{
this->SetNumberOfSteps(numSteps);
this->SetStepSize(stepSize);
this->SetSeeds(seeds);
this->Go();
}
vtkm::cont::PartitionedDataSet GetOutput() const
{
vtkm::cont::PartitionedDataSet output;
for (const auto& b : this->Blocks)
{
vtkm::cont::DataSet ds;
if (b.template GetOutput<ParticleType>(ds))
output.AppendPartition(ds);
}
return output;
}
void SetStepSize(vtkm::FloatDefault stepSize) { this->StepSize = stepSize; }
void SetNumberOfSteps(vtkm::Id numSteps) { this->NumberOfSteps = numSteps; }
void SetSeeds(const vtkm::cont::ArrayHandle<ParticleType>& seeds)
{
this->ClearParticles();
vtkm::Id n = seeds.GetNumberOfValues();
auto portal = seeds.ReadPortal();
std::vector<std::vector<vtkm::Id>> blockIDs;
std::vector<ParticleType> particles;
for (vtkm::Id i = 0; i < n; i++)
{
const ParticleType p = portal.Get(i);
std::vector<vtkm::Id> ids = this->BoundsMap.FindBlocks(p.Pos);
if (!ids.empty() && this->BoundsMap.FindRank(ids[0]) == this->Rank)
{
particles.emplace_back(p);
blockIDs.emplace_back(ids);
}
}
this->SetSeedArray(particles, blockIDs);
}
//Advect all the particles.
virtual void Go()
{
vtkm::filter::flow::internal::ParticleMessenger<ParticleType> messenger(
this->Comm, this->BoundsMap, 1, 128);
vtkm::Id nLocal = static_cast<vtkm::Id>(this->Active.size() + this->Inactive.size());
this->ComputeTotalNumParticles(nLocal);
while (this->TotalNumTerminatedParticles < this->TotalNumParticles)
{
std::vector<ParticleType> v;
vtkm::Id numTerm = 0, blockId = -1;
if (this->GetActiveParticles(v, blockId))
{
//make this a pointer to avoid the copy?
auto& block = this->GetDataSet(blockId);
DSIHelperInfoType bb =
DSIHelperInfo<ParticleType>(v, this->BoundsMap, this->ParticleBlockIDsMap);
block.Advect(bb, this->StepSize, this->NumberOfSteps);
numTerm = this->UpdateResult(bb.Get<DSIHelperInfo<ParticleType>>());
}
vtkm::Id numTermMessages = 0;
this->Communicate(messenger, numTerm, numTermMessages);
this->TotalNumTerminatedParticles += (numTerm + numTermMessages);
if (this->TotalNumTerminatedParticles > this->TotalNumParticles)
throw vtkm::cont::ErrorFilterExecution("Particle count error");
}
}
virtual void ClearParticles()
{
this->Active.clear();
this->Inactive.clear();
this->ParticleBlockIDsMap.clear();
}
void ComputeTotalNumParticles(const vtkm::Id& numLocal)
{
long long total = static_cast<long long>(numLocal);
#ifdef VTKM_ENABLE_MPI
MPI_Comm mpiComm = vtkmdiy::mpi::mpi_cast(this->Comm.handle());
MPI_Allreduce(MPI_IN_PLACE, &total, 1, MPI_LONG_LONG, MPI_SUM, mpiComm);
#endif
this->TotalNumParticles = static_cast<vtkm::Id>(total);
}
DataSetIntegrator<DSIType>& GetDataSet(vtkm::Id id)
{
for (auto& it : this->Blocks)
if (it.GetID() == id)
return it;
throw vtkm::cont::ErrorFilterExecution("Bad block");
}
virtual void SetSeedArray(const std::vector<ParticleType>& particles,
const std::vector<std::vector<vtkm::Id>>& blockIds)
{
VTKM_ASSERT(particles.size() == blockIds.size());
auto pit = particles.begin();
auto bit = blockIds.begin();
while (pit != particles.end() && bit != blockIds.end())
{
this->ParticleBlockIDsMap[pit->ID] = *bit;
pit++;
bit++;
}
this->Active.insert(this->Active.end(), particles.begin(), particles.end());
}
virtual bool GetActiveParticles(std::vector<ParticleType>& particles, vtkm::Id& blockId)
{
particles.clear();
blockId = -1;
if (this->Active.empty())
return false;
blockId = this->ParticleBlockIDsMap[this->Active.front().ID][0];
auto it = this->Active.begin();
while (it != this->Active.end())
{
auto p = *it;
if (blockId == this->ParticleBlockIDsMap[p.ID][0])
{
particles.emplace_back(p);
it = this->Active.erase(it);
}
else
it++;
}
return !particles.empty();
}
virtual void Communicate(vtkm::filter::flow::internal::ParticleMessenger<ParticleType>& messenger,
vtkm::Id numLocalTerminations,
vtkm::Id& numTermMessages)
{
std::vector<ParticleType> incoming;
std::unordered_map<vtkm::Id, std::vector<vtkm::Id>> incomingIDs;
numTermMessages = 0;
messenger.Exchange(this->Inactive,
this->ParticleBlockIDsMap,
numLocalTerminations,
incoming,
incomingIDs,
numTermMessages,
this->GetBlockAndWait(numLocalTerminations));
this->Inactive.clear();
this->UpdateActive(incoming, incomingIDs);
}
virtual void UpdateActive(const std::vector<ParticleType>& particles,
const std::unordered_map<vtkm::Id, std::vector<vtkm::Id>>& idsMap)
{
this->Update(this->Active, particles, idsMap);
}
virtual void UpdateInactive(const std::vector<ParticleType>& particles,
const std::unordered_map<vtkm::Id, std::vector<vtkm::Id>>& idsMap)
{
this->Update(this->Inactive, particles, idsMap);
}
void Update(std::vector<ParticleType>& arr,
const std::vector<ParticleType>& particles,
const std::unordered_map<vtkm::Id, std::vector<vtkm::Id>>& idsMap)
{
VTKM_ASSERT(particles.size() == idsMap.size());
arr.insert(arr.end(), particles.begin(), particles.end());
for (const auto& it : idsMap)
this->ParticleBlockIDsMap[it.first] = it.second;
}
vtkm::Id UpdateResult(const DSIHelperInfo<ParticleType>& stuff)
{
this->UpdateActive(stuff.A, stuff.IdMapA);
this->UpdateInactive(stuff.I, stuff.IdMapI);
vtkm::Id numTerm = static_cast<vtkm::Id>(stuff.TermID.size());
//Update terminated particles.
if (numTerm > 0)
{
for (const auto& id : stuff.TermID)
this->ParticleBlockIDsMap.erase(id);
}
return numTerm;
}
virtual bool GetBlockAndWait(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.
if (this->Active.empty() && this->Inactive.empty() &&
(numLocalTerm + this->TotalNumTerminatedParticles < this->TotalNumParticles))
{
return true;
}
return false;
}
//Member data
std::vector<ParticleType> Active;
std::vector<DSIType> Blocks;
vtkm::filter::flow::internal::BoundsMap BoundsMap;
vtkmdiy::mpi::communicator Comm = vtkm::cont::EnvironmentTracker::GetCommunicator();
std::vector<ParticleType> Inactive;
vtkm::Id NumberOfSteps;
vtkm::Id NumRanks;
std::unordered_map<vtkm::Id, std::vector<vtkm::Id>> ParticleBlockIDsMap;
vtkm::Id Rank;
vtkm::FloatDefault StepSize;
vtkm::Id TotalNumParticles = 0;
vtkm::Id TotalNumTerminatedParticles = 0;
};
}
}
}
} //vtkm::filter::flow::internal
#endif //vtk_m_filter_flow_internal_AdvectAlgorithm_h

@ -8,10 +8,14 @@
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#ifndef vtk_m_filter_particleadvection_AdvectorBaseThreadedAlgorithm_h
#define vtk_m_filter_particleadvection_AdvectorBaseThreadedAlgorithm_h
#ifndef vtk_m_filter_flow_internal_AdvectAlgorithmThreaded_h
#define vtk_m_filter_flow_internal_AdvectAlgorithmThreaded_h
#include <vtkm/filter/particleadvection/AdvectorBaseAlgorithm.h>
#include <vtkm/cont/PartitionedDataSet.h>
#include <vtkm/filter/flow/internal/AdvectAlgorithm.h>
#include <vtkm/filter/flow/internal/BoundsMap.h>
#include <vtkm/filter/flow/internal/DataSetIntegrator.h>
#include <vtkm/filter/flow/internal/ParticleMessenger.h>
#include <thread>
@ -19,17 +23,18 @@ namespace vtkm
{
namespace filter
{
namespace particleadvection
namespace flow
{
namespace internal
{
template <typename DataSetIntegratorType, typename ResultType>
class VTKM_ALWAYS_EXPORT AdvectorBaseThreadedAlgorithm
: public AdvectorBaseAlgorithm<DataSetIntegratorType, ResultType>
template <typename DSIType, template <typename> class ResultType, typename ParticleType>
class AdvectAlgorithmThreaded : public AdvectAlgorithm<DSIType, ResultType, ParticleType>
{
public:
AdvectorBaseThreadedAlgorithm(const vtkm::filter::particleadvection::BoundsMap& bm,
const std::vector<DataSetIntegratorType>& blocks)
: AdvectorBaseAlgorithm<DataSetIntegratorType, ResultType>(bm, blocks)
AdvectAlgorithmThreaded(const vtkm::filter::flow::internal::BoundsMap& bm,
std::vector<DSIType>& blocks)
: AdvectAlgorithm<DSIType, ResultType, ParticleType>(bm, blocks)
, Done(false)
, WorkerActivate(false)
{
@ -46,30 +51,32 @@ public:
this->ComputeTotalNumParticles(nLocal);
std::vector<std::thread> workerThreads;
workerThreads.push_back(std::thread(AdvectorBaseThreadedAlgorithm::Worker, this));
workerThreads.emplace_back(std::thread(AdvectAlgorithmThreaded::Worker, this));
this->Manage();
//This will only work for 1 thread. For > 1, the Blocks will need a mutex.
VTKM_ASSERT(workerThreads.size() == 1);
for (auto& t : workerThreads)
t.join();
}
protected:
bool GetActiveParticles(std::vector<vtkm::Particle>& particles, vtkm::Id& blockId) override
bool GetActiveParticles(std::vector<ParticleType>& particles, vtkm::Id& blockId) override
{
std::lock_guard<std::mutex> lock(this->Mutex);
bool val = this->AdvectorBaseAlgorithm<DataSetIntegratorType, ResultType>::GetActiveParticles(
bool val = this->AdvectAlgorithm<DSIType, ResultType, ParticleType>::GetActiveParticles(
particles, blockId);
this->WorkerActivate = val;
return val;
}
void UpdateActive(const std::vector<vtkm::Particle>& particles,
void UpdateActive(const std::vector<ParticleType>& particles,
const std::unordered_map<vtkm::Id, std::vector<vtkm::Id>>& idsMap) override
{
if (!particles.empty())
{
std::lock_guard<std::mutex> lock(this->Mutex);
this->AdvectorBaseAlgorithm<DataSetIntegratorType, ResultType>::UpdateActive(particles,
idsMap);
this->AdvectAlgorithm<DSIType, ResultType, ParticleType>::UpdateActive(particles, idsMap);
//Let workers know there is new work
this->WorkerActivateCondition.notify_all();
@ -90,7 +97,7 @@ protected:
this->WorkerActivateCondition.notify_all();
}
static void Worker(AdvectorBaseThreadedAlgorithm* algo) { algo->Work(); }
static void Worker(AdvectAlgorithmThreaded* algo) { algo->Work(); }
void WorkerWait()
{
@ -98,18 +105,26 @@ protected:
this->WorkerActivateCondition.wait(lock, [this] { return WorkerActivate || Done; });
}
void UpdateWorkerResult(vtkm::Id blockId, DSIHelperInfoType& b)
{
std::lock_guard<std::mutex> lock(this->Mutex);
auto& it = this->WorkerResults[blockId];
it.emplace_back(b);
}
void Work()
{
while (!this->CheckDone())
{
std::vector<vtkm::Particle> v;
std::vector<ParticleType> v;
vtkm::Id blockId = -1;
if (this->GetActiveParticles(v, blockId))
{
const auto& block = this->GetDataSet(blockId);
ResultType r;
block.Advect(v, this->StepSize, this->NumberOfSteps, r);
this->UpdateWorkerResult(blockId, r);
auto& block = this->GetDataSet(blockId);
DSIHelperInfoType bb =
DSIHelperInfo<ParticleType>(v, this->BoundsMap, this->ParticleBlockIDsMap);
block.Advect(bb, this->StepSize, this->NumberOfSteps);
this->UpdateWorkerResult(blockId, bb);
}
else
this->WorkerWait();
@ -118,20 +133,19 @@ protected:
void Manage()
{
vtkm::filter::particleadvection::ParticleMessenger messenger(
vtkm::filter::flow::internal::ParticleMessenger<ParticleType> messenger(
this->Comm, this->BoundsMap, 1, 128);
while (this->TotalNumTerminatedParticles < this->TotalNumParticles)
{
std::unordered_map<vtkm::Id, std::vector<ResultType>> workerResults;
std::unordered_map<vtkm::Id, std::vector<DSIHelperInfoType>> workerResults;
this->GetWorkerResults(workerResults);
vtkm::Id numTerm = 0;
for (auto& it : workerResults)
{
vtkm::Id blockId = it.first;
for (auto& r : it.second)
numTerm += this->UpdateResult(r, blockId);
numTerm += this->UpdateResult(r.Get<DSIHelperInfo<ParticleType>>());
}
vtkm::Id numTermMessages = 0;
@ -150,12 +164,12 @@ protected:
{
std::lock_guard<std::mutex> lock(this->Mutex);
return (this->AdvectorBaseAlgorithm<DataSetIntegratorType, ResultType>::GetBlockAndWait(
numLocalTerm) &&
!this->WorkerActivate && this->WorkerResults.empty());
return (
this->AdvectAlgorithm<DSIType, ResultType, ParticleType>::GetBlockAndWait(numLocalTerm) &&
!this->WorkerActivate && this->WorkerResults.empty());
}
void GetWorkerResults(std::unordered_map<vtkm::Id, std::vector<ResultType>>& results)
void GetWorkerResults(std::unordered_map<vtkm::Id, std::vector<DSIHelperInfoType>>& results)
{
results.clear();
@ -167,23 +181,16 @@ protected:
}
}
void UpdateWorkerResult(vtkm::Id blockId, const ResultType& result)
{
std::lock_guard<std::mutex> lock(this->Mutex);
auto& it = this->WorkerResults[blockId];
it.push_back(result);
}
std::atomic<bool> Done;
std::mutex Mutex;
bool WorkerActivate;
std::condition_variable WorkerActivateCondition;
std::unordered_map<vtkm::Id, std::vector<ResultType>> WorkerResults;
std::unordered_map<vtkm::Id, std::vector<DSIHelperInfoType>> WorkerResults;
};
}
}
} // namespace vtkm::filter::particleadvection
}
} //vtkm::filter::flow::internal
#endif
#endif //vtk_m_filter_flow_internal_AdvectAlgorithmThreaded_h

@ -8,8 +8,8 @@
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#ifndef vtk_m_filter_BoundsMap_h
#define vtk_m_filter_BoundsMap_h
#ifndef vtk_m_filter_flow_internal_BoundsMap_h
#define vtk_m_filter_flow_internal_BoundsMap_h
#include <vtkm/Bounds.h>
#include <vtkm/cont/AssignerPartitionedDataSet.h>
@ -31,7 +31,9 @@ namespace vtkm
{
namespace filter
{
namespace particleadvection
namespace flow
{
namespace internal
{
class VTKM_ALWAYS_EXPORT BoundsMap
@ -96,7 +98,7 @@ public:
for (auto& it : this->BlockBounds)
{
if (blockId != ignoreBlock && it.Contains(p))
blockIDs.push_back(blockId);
blockIDs.emplace_back(blockId);
blockId++;
}
}
@ -117,7 +119,7 @@ private:
vtkmdiy::mpi::communicator Comm = vtkm::cont::EnvironmentTracker::GetCommunicator();
assigner.local_gids(Comm.rank(), ids);
for (const auto& i : ids)
this->LocalIDs.push_back(static_cast<vtkm::Id>(i));
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));
@ -174,8 +176,10 @@ private:
std::vector<vtkm::Bounds> BlockBounds;
vtkm::Bounds GlobalBounds;
};
}
}
} // namespace vtkm::filter::particleadvection
#endif
}
}
}
} // namespace vtkm::filter::flow::internal
#endif //vtk_m_filter_flow_internal_BoundsMap_h

@ -9,16 +9,18 @@
##============================================================================
set(headers
AdvectorBaseAlgorithm.h
AdvectorBaseAlgorithm.hxx
AdvectorBaseThreadedAlgorithm.h
AdvectAlgorithm.h
AdvectAlgorithmThreaded.h
BoundsMap.h
DataSetIntegrator.h
DataSetIntegrator.hxx
DataSetIntegratorSteadyState.h
DataSetIntegratorUnsteadyState.h
Messenger.h
ParticleAdvectionAlgorithm.h
ParticleAdvector.h
ParticleMessenger.h
)
)
# Note: The C++ source files are added to the flow library
# in the CMakeLists.txt in the parent directory.
#-----------------------------------------------------------------------------
vtkm_declare_headers(${headers})

@ -0,0 +1,375 @@
//============================================================================
// 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_flow_internal_DataSetIntegrator_h
#define vtk_m_filter_flow_internal_DataSetIntegrator_h
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/ErrorFilterExecution.h>
#include <vtkm/cont/ParticleArrayCopy.h>
#include <vtkm/filter/flow/FlowTypes.h>
#include <vtkm/filter/flow/internal/BoundsMap.h>
#include <vtkm/filter/flow/worklet/EulerIntegrator.h>
#include <vtkm/filter/flow/worklet/IntegratorStatus.h>
#include <vtkm/filter/flow/worklet/ParticleAdvection.h>
#include <vtkm/filter/flow/worklet/RK4Integrator.h>
#include <vtkm/filter/flow/worklet/Stepper.h>
#include <vtkm/cont/internal/Variant.h>
namespace vtkm
{
namespace filter
{
namespace flow
{
namespace internal
{
template <typename ParticleType>
class DSIHelperInfo
{
public:
DSIHelperInfo(const std::vector<ParticleType>& v,
const vtkm::filter::flow::internal::BoundsMap& boundsMap,
const std::unordered_map<vtkm::Id, std::vector<vtkm::Id>>& particleBlockIDsMap)
: BoundsMap(boundsMap)
, ParticleBlockIDsMap(particleBlockIDsMap)
, V(v)
{
}
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;
};
using DSIHelperInfoType = vtkm::cont::internal::Variant<DSIHelperInfo<vtkm::Particle>,
DSIHelperInfo<vtkm::ChargedParticle>>;
template <typename Derived>
class DataSetIntegrator
{
protected:
using VelocityFieldNameType = std::string;
using ElectroMagneticFieldNameType = std::pair<std::string, std::string>;
using FieldNameType =
vtkm::cont::internal::Variant<VelocityFieldNameType, ElectroMagneticFieldNameType>;
using RType = vtkm::cont::internal::Variant<
vtkm::worklet::flow::ParticleAdvectionResult<vtkm::Particle>,
vtkm::worklet::flow::ParticleAdvectionResult<vtkm::ChargedParticle>,
vtkm::worklet::flow::StreamlineResult<vtkm::Particle>,
vtkm::worklet::flow::StreamlineResult<vtkm::ChargedParticle>>;
public:
DataSetIntegrator(vtkm::Id id,
const FieldNameType& fieldName,
vtkm::filter::flow::IntegrationSolverType solverType,
vtkm::filter::flow::VectorFieldType vecFieldType,
vtkm::filter::flow::FlowResultType resultType)
: FieldName(fieldName)
, Id(id)
, SolverType(solverType)
, VecFieldType(vecFieldType)
, AdvectionResType(resultType)
, Rank(this->Comm.rank())
{
//check that things are valid.
}
VTKM_CONT vtkm::Id GetID() const { return this->Id; }
VTKM_CONT void SetCopySeedFlag(bool val) { this->CopySeedArray = val; }
VTKM_CONT
void Advect(DSIHelperInfoType& b,
vtkm::FloatDefault stepSize, //move these to member data(?)
vtkm::Id maxSteps)
{
Derived* inst = static_cast<Derived*>(this);
//Cast the DSIHelperInfo<ParticleType> to the concrete type and call DoAdvect.
b.CastAndCall([&](auto& concrete) { inst->DoAdvect(concrete, stepSize, maxSteps); });
}
template <typename ParticleType>
VTKM_CONT bool GetOutput(vtkm::cont::DataSet& ds) const;
protected:
template <typename ParticleType, template <typename> class ResultType>
VTKM_CONT void UpdateResult(const ResultType<ParticleType>& result,
DSIHelperInfo<ParticleType>& dsiInfo);
VTKM_CONT bool IsParticleAdvectionResult() const
{
return this->AdvectionResType == FlowResultType::PARTICLE_ADVECT_TYPE;
}
VTKM_CONT bool IsStreamlineResult() const
{
return this->AdvectionResType == FlowResultType::STREAMLINE_TYPE;
}
template <typename ParticleType>
VTKM_CONT inline void ClassifyParticles(const vtkm::cont::ArrayHandle<ParticleType>& particles,
DSIHelperInfo<ParticleType>& dsiInfo) const;
//Data members.
vtkm::cont::internal::Variant<VelocityFieldNameType, ElectroMagneticFieldNameType> FieldName;
vtkm::Id Id;
vtkm::filter::flow::IntegrationSolverType SolverType;
vtkm::filter::flow::VectorFieldType VecFieldType;
vtkm::filter::flow::FlowResultType AdvectionResType =
vtkm::filter::flow::FlowResultType::UNKNOWN_TYPE;
vtkmdiy::mpi::communicator Comm = vtkm::cont::EnvironmentTracker::GetCommunicator();
vtkm::Id Rank;
bool CopySeedArray = false;
std::vector<RType> Results;
};
template <typename Derived>
template <typename ParticleType>
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();
auto portal = particles.WritePortal();
vtkm::Id n = portal.GetNumberOfValues();
for (vtkm::Id i = 0; i < n; i++)
{
auto p = portal.Get(i);
if (p.Status.CheckTerminate())
{
dsiInfo.TermIdx.emplace_back(i);
dsiInfo.TermID.emplace_back(p.ID);
}
else
{
const auto& it = dsiInfo.ParticleBlockIDsMap.find(p.ID);
VTKM_ASSERT(it != dsiInfo.ParticleBlockIDsMap.end());
auto currBIDs = it->second;
VTKM_ASSERT(!currBIDs.empty());
std::vector<vtkm::Id> newIDs;
if (p.Status.CheckSpatialBounds() && !p.Status.CheckTookAnySteps())
newIDs.assign(std::next(currBIDs.begin(), 1), currBIDs.end());
else
newIDs = dsiInfo.BoundsMap.FindBlocks(p.Pos, currBIDs);
//reset the particle status.
p.Status = vtkm::ParticleStatus();
if (newIDs.empty()) //No blocks, we're done.
{
p.Status.SetTerminate();
dsiInfo.TermIdx.emplace_back(i);
dsiInfo.TermID.emplace_back(p.ID);
}
else
{
//If we have more than 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)
{
newIDs.erase(idit);
newIDs.insert(newIDs.begin(), bid);
break;
}
}
}
int dstRank = dsiInfo.BoundsMap.FindRank(newIDs[0]);
if (dstRank == this->Rank)
{
dsiInfo.A.emplace_back(p);
dsiInfo.IdMapA[p.ID] = newIDs;
}
else
{
dsiInfo.I.emplace_back(p);
dsiInfo.IdMapI[p.ID] = newIDs;
}
}
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());
}
template <typename Derived>
template <typename ParticleType, template <typename> class ResultType>
VTKM_CONT inline void DataSetIntegrator<Derived>::UpdateResult(
const ResultType<ParticleType>& result,
DSIHelperInfo<ParticleType>& dsiInfo)
{
this->ClassifyParticles(result.Particles, dsiInfo);
if (this->IsParticleAdvectionResult())
{
if (dsiInfo.TermIdx.empty())
return;
using ResType = vtkm::worklet::flow::ParticleAdvectionResult<ParticleType>;
auto indicesAH = vtkm::cont::make_ArrayHandle(dsiInfo.TermIdx, vtkm::CopyFlag::Off);
auto termPerm = vtkm::cont::make_ArrayHandlePermutation(indicesAH, result.Particles);
vtkm::cont::ArrayHandle<ParticleType> termParticles;
vtkm::cont::Algorithm::Copy(termPerm, termParticles);
ResType termRes(termParticles);
this->Results.emplace_back(termRes);
}
else if (this->IsStreamlineResult())
this->Results.emplace_back(result);
}
template <typename Derived>
template <typename ParticleType>
VTKM_CONT inline bool DataSetIntegrator<Derived>::GetOutput(vtkm::cont::DataSet& ds) const
{
std::size_t nResults = this->Results.size();
if (nResults == 0)
return false;
if (this->IsParticleAdvectionResult())
{
using ResType = vtkm::worklet::flow::ParticleAdvectionResult<ParticleType>;
std::vector<vtkm::cont::ArrayHandle<ParticleType>> allParticles;
allParticles.reserve(nResults);
for (const auto& vres : this->Results)
allParticles.emplace_back(vres.template Get<ResType>().Particles);
vtkm::cont::ArrayHandle<vtkm::Vec3f> pts;
vtkm::cont::ParticleArrayCopy(allParticles, pts);
vtkm::Id numPoints = pts.GetNumberOfValues();
if (numPoints > 0)
{
//Create coordinate system and vertex cell set.
ds.AddCoordinateSystem(vtkm::cont::CoordinateSystem("coordinates", pts));
vtkm::cont::CellSetSingleType<> cells;
vtkm::cont::ArrayHandleIndex conn(numPoints);
vtkm::cont::ArrayHandle<vtkm::Id> connectivity;
vtkm::cont::ArrayCopy(conn, connectivity);
cells.Fill(numPoints, vtkm::CELL_SHAPE_VERTEX, 1, connectivity);
ds.SetCellSet(cells);
}
}
else if (this->IsStreamlineResult())
{
using ResType = vtkm::worklet::flow::StreamlineResult<ParticleType>;
//Easy case with one result.
if (nResults == 1)
{
const auto& res = this->Results[0].template Get<ResType>();
ds.AddCoordinateSystem(vtkm::cont::CoordinateSystem("coordinates", res.Positions));
ds.SetCellSet(res.PolyLines);
}
else
{
std::vector<vtkm::Id> posOffsets(nResults, 0);
vtkm::Id totalNumCells = 0, totalNumPts = 0;
for (std::size_t i = 0; i < nResults; i++)
{
const auto& res = this->Results[i].template Get<ResType>();
if (i == 0)
posOffsets[i] = 0;
else
posOffsets[i] = totalNumPts;
totalNumPts += res.Positions.GetNumberOfValues();
totalNumCells += res.PolyLines.GetNumberOfCells();
}
//Append all the points together.
vtkm::cont::ArrayHandle<vtkm::Vec3f> appendPts;
appendPts.Allocate(totalNumPts);
for (std::size_t i = 0; i < nResults; i++)
{
const auto& res = this->Results[i].template Get<ResType>();
// copy all values into appendPts starting at offset.
vtkm::cont::Algorithm::CopySubRange(
res.Positions, 0, res.Positions.GetNumberOfValues(), appendPts, posOffsets[i]);
}
ds.AddCoordinateSystem(vtkm::cont::CoordinateSystem("coordinates", appendPts));
//Create polylines.
std::vector<vtkm::Id> numPtsPerCell(static_cast<std::size_t>(totalNumCells));
std::size_t off = 0;
for (std::size_t i = 0; i < nResults; i++)
{
const auto& res = this->Results[i].template Get<ResType>();
vtkm::Id nCells = res.PolyLines.GetNumberOfCells();
for (vtkm::Id j = 0; j < nCells; j++)
numPtsPerCell[off++] = static_cast<vtkm::Id>(res.PolyLines.GetNumberOfPointsInCell(j));
}
auto numPointsPerCellArray = vtkm::cont::make_ArrayHandle(numPtsPerCell, vtkm::CopyFlag::Off);
vtkm::cont::ArrayHandle<vtkm::Id> cellIndex;
vtkm::Id connectivityLen =
vtkm::cont::Algorithm::ScanExclusive(numPointsPerCellArray, cellIndex);
vtkm::cont::ArrayHandleIndex connCount(connectivityLen);
vtkm::cont::ArrayHandle<vtkm::Id> connectivity;
vtkm::cont::ArrayCopy(connCount, connectivity);
vtkm::cont::ArrayHandle<vtkm::UInt8> cellTypes;
auto polyLineShape = vtkm::cont::make_ArrayHandleConstant<vtkm::UInt8>(
vtkm::CELL_SHAPE_POLY_LINE, totalNumCells);
vtkm::cont::ArrayCopy(polyLineShape, cellTypes);
auto offsets = vtkm::cont::ConvertNumComponentsToOffsets(numPointsPerCellArray);
vtkm::cont::CellSetExplicit<> polyLines;
polyLines.Fill(totalNumPts, cellTypes, connectivity, offsets);
ds.SetCellSet(polyLines);
}
}
else
{
throw vtkm::cont::ErrorFilterExecution("Unsupported ParticleAdvectionResultType");
}
return true;
}
}
}
}
} //vtkm::filter::flow::internal
#endif //vtk_m_filter_flow_internal_DataSetIntegrator_h

@ -0,0 +1,214 @@
//============================================================================
// 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_flow_internal_DataSetIntegratorSteadyState_h
#define vtk_m_filter_flow_internal_DataSetIntegratorSteadyState_h
#include <vtkm/cont/ArrayCopy.h>
#include <vtkm/filter/flow/internal/DataSetIntegrator.h>
namespace vtkm
{
namespace filter
{
namespace flow
{
namespace internal
{
class DataSetIntegratorSteadyState
: public vtkm::filter::flow::internal::DataSetIntegrator<DataSetIntegratorSteadyState>
{
public:
DataSetIntegratorSteadyState(const vtkm::cont::DataSet& ds,
vtkm::Id id,
const FieldNameType& fieldName,
vtkm::filter::flow::IntegrationSolverType solverType,
vtkm::filter::flow::VectorFieldType vecFieldType,
vtkm::filter::flow::FlowResultType resultType)
: vtkm::filter::flow::internal::DataSetIntegrator<DataSetIntegratorSteadyState>(id,
fieldName,
solverType,
vecFieldType,
resultType)
, DataSet(ds)
{
}
VTKM_CONT inline void DoAdvect(DSIHelperInfo<vtkm::Particle>& b,
vtkm::FloatDefault stepSize,
vtkm::Id maxSteps);
VTKM_CONT inline void DoAdvect(DSIHelperInfo<vtkm::ChargedParticle>& b,
vtkm::FloatDefault stepSize,
vtkm::Id maxSteps);
protected:
template <typename ArrayType>
VTKM_CONT void GetVelocityField(
vtkm::worklet::flow::VelocityField<ArrayType>& velocityField) const
{
if (this->FieldName.GetIndex() == this->FieldName.GetIndexOf<VelocityFieldNameType>())
{
const auto& fieldNm = this->FieldName.Get<VelocityFieldNameType>();
auto assoc = this->DataSet.GetField(fieldNm).GetAssociation();
ArrayType arr;
vtkm::cont::ArrayCopyShallowIfPossible(this->DataSet.GetField(fieldNm).GetData(), arr);
velocityField = vtkm::worklet::flow::VelocityField<ArrayType>(arr, assoc);
}
else
throw vtkm::cont::ErrorFilterExecution("Velocity field vector type not available");
}
private:
vtkm::cont::DataSet DataSet;
};
namespace internal
{
using ArrayType = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
using VelocityFieldType = vtkm::worklet::flow::VelocityField<ArrayType>;
using SteadyStateGridEvalType = vtkm::worklet::flow::GridEvaluator<VelocityFieldType>;
template <typename GridEvalType, typename ParticleType>
class AdvectHelper;
template <typename ParticleType>
class AdvectHelper<SteadyStateGridEvalType, ParticleType>
{
public:
static void Advect(const VelocityFieldType& velField,
const vtkm::cont::DataSet& ds,
vtkm::cont::ArrayHandle<ParticleType>& seedArray,
vtkm::FloatDefault stepSize,
vtkm::Id maxSteps,
const IntegrationSolverType& solverType,
vtkm::worklet::flow::ParticleAdvectionResult<ParticleType>& result)
{
if (solverType == IntegrationSolverType::RK4_TYPE)
{
DoAdvect<vtkm::worklet::flow::ParticleAdvection,
vtkm::worklet::flow::ParticleAdvectionResult,
vtkm::worklet::flow::RK4Integrator>(
velField, ds, seedArray, stepSize, maxSteps, result);
}
else if (solverType == IntegrationSolverType::EULER_TYPE)
{
DoAdvect<vtkm::worklet::flow::ParticleAdvection,
vtkm::worklet::flow::ParticleAdvectionResult,
vtkm::worklet::flow::EulerIntegrator>(
velField, ds, seedArray, stepSize, maxSteps, result);
}
else
throw vtkm::cont::ErrorFilterExecution("Unsupported Integrator type");
}
static void Advect(const VelocityFieldType& velField,
const vtkm::cont::DataSet& ds,
vtkm::cont::ArrayHandle<ParticleType>& seedArray,
vtkm::FloatDefault stepSize,
vtkm::Id maxSteps,
const IntegrationSolverType& solverType,
vtkm::worklet::flow::StreamlineResult<ParticleType>& result)
{
if (solverType == IntegrationSolverType::RK4_TYPE)
{
DoAdvect<vtkm::worklet::flow::Streamline,
vtkm::worklet::flow::StreamlineResult,
vtkm::worklet::flow::RK4Integrator>(
velField, ds, seedArray, stepSize, maxSteps, result);
}
else if (solverType == IntegrationSolverType::EULER_TYPE)
{
DoAdvect<vtkm::worklet::flow::Streamline,
vtkm::worklet::flow::StreamlineResult,
vtkm::worklet::flow::EulerIntegrator>(
velField, ds, seedArray, stepSize, maxSteps, result);
}
else
throw vtkm::cont::ErrorFilterExecution("Unsupported Integrator type");
}
template <typename WorkletType,
template <typename>
class ResultType,
template <typename>
class SolverType>
static void DoAdvect(const VelocityFieldType& velField,
const vtkm::cont::DataSet& ds,
vtkm::cont::ArrayHandle<ParticleType>& seedArray,
vtkm::FloatDefault stepSize,
vtkm::Id maxSteps,
ResultType<ParticleType>& result)
{
using StepperType =
vtkm::worklet::flow::Stepper<SolverType<SteadyStateGridEvalType>, SteadyStateGridEvalType>;
WorkletType worklet;
SteadyStateGridEvalType eval(ds, velField);
StepperType stepper(eval, stepSize);
result = worklet.Run(stepper, seedArray, maxSteps);
}
};
}
VTKM_CONT inline void DataSetIntegratorSteadyState::DoAdvect(DSIHelperInfo<vtkm::Particle>& b,
vtkm::FloatDefault stepSize,
vtkm::Id maxSteps)
{
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);
if (this->VecFieldType == VectorFieldType::VELOCITY_FIELD_TYPE)
{
using FieldType = vtkm::worklet::flow::VelocityField<ArrayType>;
FieldType velField;
this->GetVelocityField(velField);
using AHType = internal::AdvectHelper<internal::SteadyStateGridEvalType, vtkm::Particle>;
if (this->IsParticleAdvectionResult())
{
vtkm::worklet::flow::ParticleAdvectionResult<vtkm::Particle> result;
AHType::Advect(
velField, this->DataSet, seedArray, stepSize, maxSteps, this->SolverType, result);
this->UpdateResult(result, b);
}
else if (this->IsStreamlineResult())
{
vtkm::worklet::flow::StreamlineResult<vtkm::Particle> result;
AHType::Advect(
velField, this->DataSet, seedArray, stepSize, maxSteps, this->SolverType, result);
this->UpdateResult(result, b);
}
else
throw vtkm::cont::ErrorFilterExecution("Unsupported result type");
}
else
throw vtkm::cont::ErrorFilterExecution("Unsupported vector field type");
}
VTKM_CONT inline void DataSetIntegratorSteadyState::DoAdvect(
DSIHelperInfo<vtkm::ChargedParticle>& vtkmNotUsed(b),
vtkm::FloatDefault vtkmNotUsed(stepSize),
vtkm::Id vtkmNotUsed(maxSteps))
{
}
}
}
}
} //vtkm::filter::flow::internal
#endif //vtk_m_filter_flow_internal_DataSetIntegratorSteadyState_h

@ -0,0 +1,265 @@
//============================================================================
// 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_flow_internal_DataSetIntegratorUnsteadyState_h
#define vtk_m_filter_flow_internal_DataSetIntegratorUnsteadyState_h
#include <vtkm/cont/ArrayCopy.h>
#include <vtkm/filter/flow/internal/DataSetIntegrator.h>
#include <vtkm/filter/flow/worklet/TemporalGridEvaluators.h>
namespace vtkm
{
namespace filter
{
namespace flow
{
namespace internal
{
class DataSetIntegratorUnsteadyState
: public vtkm::filter::flow::internal::DataSetIntegrator<DataSetIntegratorUnsteadyState>
{
public:
DataSetIntegratorUnsteadyState(const vtkm::cont::DataSet& ds1,
const vtkm::cont::DataSet& ds2,
vtkm::FloatDefault t1,
vtkm::FloatDefault t2,
vtkm::Id id,
const vtkm::filter::flow::internal::DataSetIntegrator<
DataSetIntegratorUnsteadyState>::FieldNameType& fieldName,
vtkm::filter::flow::IntegrationSolverType solverType,
vtkm::filter::flow::VectorFieldType vecFieldType,
vtkm::filter::flow::FlowResultType resultType)
: vtkm::filter::flow::internal::DataSetIntegrator<DataSetIntegratorUnsteadyState>(id,
fieldName,
solverType,
vecFieldType,
resultType)
, DataSet1(ds1)
, DataSet2(ds2)
, Time1(t1)
, Time2(t2)
{
}
VTKM_CONT inline void DoAdvect(DSIHelperInfo<vtkm::Particle>& b,
vtkm::FloatDefault stepSize,
vtkm::Id maxSteps);
VTKM_CONT inline void DoAdvect(DSIHelperInfo<vtkm::ChargedParticle>& b,
vtkm::FloatDefault stepSize,
vtkm::Id maxSteps);
protected:
template <typename ArrayType>
VTKM_CONT void GetVelocityFields(
vtkm::worklet::flow::VelocityField<ArrayType>& velocityField1,
vtkm::worklet::flow::VelocityField<ArrayType>& velocityField2) const
{
if (this->FieldName.GetIndex() == this->FieldName.GetIndexOf<VelocityFieldNameType>())
{
const auto& fieldNm = this->FieldName.Get<VelocityFieldNameType>();
auto assoc = this->DataSet1.GetField(fieldNm).GetAssociation();
if (assoc != this->DataSet2.GetField(fieldNm).GetAssociation())
throw vtkm::cont::ErrorFilterExecution(
"Unsteady state velocity fields have differnt associations");
ArrayType arr1, arr2;
vtkm::cont::ArrayCopyShallowIfPossible(this->DataSet1.GetField(fieldNm).GetData(), arr1);
vtkm::cont::ArrayCopyShallowIfPossible(this->DataSet2.GetField(fieldNm).GetData(), arr2);
velocityField1 = vtkm::worklet::flow::VelocityField<ArrayType>(arr1, assoc);
velocityField2 = vtkm::worklet::flow::VelocityField<ArrayType>(arr2, assoc);
}
else
throw vtkm::cont::ErrorFilterExecution("Velocity field vector type not available");
}
private:
vtkm::cont::DataSet DataSet1;
vtkm::cont::DataSet DataSet2;
vtkm::FloatDefault Time1;
vtkm::FloatDefault Time2;
};
namespace internal
{
using ArrayType = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
using VelocityFieldType = vtkm::worklet::flow::VelocityField<ArrayType>;
using UnsteadyStateGridEvalType = vtkm::worklet::flow::TemporalGridEvaluator<VelocityFieldType>;
template <typename GridEvalType, typename ParticleType>
class AdvectHelper;
template <typename ParticleType>
class AdvectHelper<UnsteadyStateGridEvalType, ParticleType>
{
public:
static void Advect(const VelocityFieldType& velField1,
const vtkm::cont::DataSet& ds1,
vtkm::FloatDefault t1,
const VelocityFieldType& velField2,
const vtkm::cont::DataSet& ds2,
vtkm::FloatDefault t2,
vtkm::cont::ArrayHandle<ParticleType>& seedArray,
vtkm::FloatDefault stepSize,
vtkm::Id maxSteps,
const IntegrationSolverType& solverType,
vtkm::worklet::flow::ParticleAdvectionResult<ParticleType>& result)
{
if (solverType == IntegrationSolverType::RK4_TYPE)
{
DoAdvect<vtkm::worklet::flow::ParticleAdvection,
vtkm::worklet::flow::ParticleAdvectionResult,
vtkm::worklet::flow::RK4Integrator>(
velField1, ds1, t1, velField2, ds2, t2, seedArray, stepSize, maxSteps, result);
}
else if (solverType == IntegrationSolverType::EULER_TYPE)
{
DoAdvect<vtkm::worklet::flow::ParticleAdvection,
vtkm::worklet::flow::ParticleAdvectionResult,
vtkm::worklet::flow::EulerIntegrator>(
velField1, ds1, t1, velField2, ds2, t2, seedArray, stepSize, maxSteps, result);
}
else
throw vtkm::cont::ErrorFilterExecution("Unsupported Integrator type");
}
static void Advect(const VelocityFieldType& velField1,
const vtkm::cont::DataSet& ds1,
vtkm::FloatDefault t1,
const VelocityFieldType& velField2,
const vtkm::cont::DataSet& ds2,
vtkm::FloatDefault t2,
vtkm::cont::ArrayHandle<ParticleType>& seedArray,
vtkm::FloatDefault stepSize,
vtkm::Id maxSteps,
const IntegrationSolverType& solverType,
vtkm::worklet::flow::StreamlineResult<ParticleType>& result)
{
if (solverType == IntegrationSolverType::RK4_TYPE)
{
DoAdvect<vtkm::worklet::flow::Streamline,
vtkm::worklet::flow::StreamlineResult,
vtkm::worklet::flow::RK4Integrator>(
velField1, ds1, t1, velField2, ds2, t2, seedArray, stepSize, maxSteps, result);
}
else if (solverType == IntegrationSolverType::EULER_TYPE)
{
DoAdvect<vtkm::worklet::flow::Streamline,
vtkm::worklet::flow::StreamlineResult,
vtkm::worklet::flow::EulerIntegrator>(
velField1, ds1, t1, velField2, ds2, t2, seedArray, stepSize, maxSteps, result);
}
else
throw vtkm::cont::ErrorFilterExecution("Unsupported Integrator type");
}
template <typename WorkletType,
template <typename>
class ResultType,
template <typename>
class SolverType>
static void DoAdvect(const VelocityFieldType& velField1,
const vtkm::cont::DataSet& ds1,
vtkm::FloatDefault t1,
const VelocityFieldType& velField2,
const vtkm::cont::DataSet& ds2,
vtkm::FloatDefault t2,
vtkm::cont::ArrayHandle<ParticleType>& seedArray,
vtkm::FloatDefault stepSize,
vtkm::Id maxSteps,
ResultType<ParticleType>& result)
{
using StepperType = vtkm::worklet::flow::Stepper<SolverType<UnsteadyStateGridEvalType>,
UnsteadyStateGridEvalType>;
WorkletType worklet;
UnsteadyStateGridEvalType eval(ds1, t1, velField1, ds2, t2, velField2);
StepperType stepper(eval, stepSize);
result = worklet.Run(stepper, seedArray, maxSteps);
}
};
}
VTKM_CONT inline void DataSetIntegratorUnsteadyState::DoAdvect(DSIHelperInfo<vtkm::Particle>& b,
vtkm::FloatDefault stepSize,
vtkm::Id maxSteps)
{
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);
using AHType = internal::AdvectHelper<internal::UnsteadyStateGridEvalType, vtkm::Particle>;
if (this->VecFieldType == VectorFieldType::VELOCITY_FIELD_TYPE)
{
using FieldType = vtkm::worklet::flow::VelocityField<ArrayType>;
FieldType velField1, velField2;
this->GetVelocityFields(velField1, velField2);
if (this->IsParticleAdvectionResult())
{
vtkm::worklet::flow::ParticleAdvectionResult<vtkm::Particle> result;
AHType::Advect(velField1,
this->DataSet1,
this->Time1,
velField2,
this->DataSet2,
this->Time2,
seedArray,
stepSize,
maxSteps,
this->SolverType,
result);
this->UpdateResult(result, b);
}
else if (this->IsStreamlineResult())
{
vtkm::worklet::flow::StreamlineResult<vtkm::Particle> result;
AHType::Advect(velField1,
this->DataSet1,
this->Time1,
velField2,
this->DataSet2,
this->Time2,
seedArray,
stepSize,
maxSteps,
this->SolverType,
result);
this->UpdateResult(result, b);
}
else
throw vtkm::cont::ErrorFilterExecution("Unsupported result type");
}
else
throw vtkm::cont::ErrorFilterExecution("Unsupported vector field type");
}
VTKM_CONT inline void DataSetIntegratorUnsteadyState::DoAdvect(
DSIHelperInfo<vtkm::ChargedParticle>& vtkmNotUsed(b),
vtkm::FloatDefault vtkmNotUsed(stepSize),
vtkm::Id vtkmNotUsed(maxSteps))
{
}
}
}
}
} //vtkm::filter::flow::internal
#endif //vtk_m_filter_flow_internal_DataSetIntegratorUnsteadyState_h

@ -8,11 +8,12 @@
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#include <vtkm/cont/ErrorFilterExecution.h>
#include <vtkm/filter/flow/internal/Messenger.h>
#include <iostream>
#include <sstream>
#include <string.h>
#include <vtkm/cont/ErrorFilterExecution.h>
#include <vtkm/filter/particleadvection/Messenger.h>
#ifdef VTKM_ENABLE_MPI
#include <vtkm/thirdparty/diy/mpi-cast.h>
@ -22,7 +23,9 @@ namespace vtkm
{
namespace filter
{
namespace particleadvection
namespace flow
{
namespace internal
{
VTKM_CONT
@ -80,7 +83,7 @@ void Messenger::CleanupRequests(int tag)
for (const auto& i : this->RecvBuffers)
{
if (tag == TAG_ANY || tag == i.first.second)
delKeys.push_back(i.first);
delKeys.emplace_back(i.first);
}
if (!delKeys.empty())
@ -150,9 +153,9 @@ void Messenger::CheckRequests(const std::map<RequestTagPair, char*>& buffers,
{
if (tagsToCheck.empty() || tagsToCheck.find(it.first.second) != tagsToCheck.end())
{
req.push_back(it.first.first);
copy.push_back(it.first.first);
tags.push_back(it.first.second);
req.emplace_back(it.first.first);
copy.emplace_back(it.first.first);
tags.emplace_back(it.first.second);
}
}
@ -176,7 +179,7 @@ void Messenger::CheckRequests(const std::map<RequestTagPair, char*>& buffers,
//Add the req/tag to the return vector.
reqTags.reserve(static_cast<std::size_t>(num));
for (int i = 0; i < num; i++)
reqTags.push_back(RequestTagPair(copy[indices[i]], tags[indices[i]]));
reqTags.emplace_back(RequestTagPair(copy[indices[i]], tags[indices[i]]));
}
bool Messenger::PacketCompare(const char* a, const char* b)
@ -300,7 +303,7 @@ bool Messenger::RecvData(const std::set<int>& tags,
if (it == this->RecvBuffers.end())
throw vtkm::cont::ErrorFilterExecution("receive buffer not found");
incomingBuffers.push_back(it->second);
incomingBuffers.emplace_back(it->second);
this->RecvBuffers.erase(it);
}
@ -329,7 +332,7 @@ void Messenger::ProcessReceivedBuffers(std::vector<char*>& incomingBuffers,
entry.first = header.tag;
entry.second.save_binary((char*)(buff + sizeof(header)), header.dataSz);
entry.second.reset();
buffers.push_back(std::move(entry));
buffers.emplace_back(std::move(entry));
delete[] buff;
}
@ -344,12 +347,12 @@ void Messenger::ProcessReceivedBuffers(std::vector<char*>& incomingBuffers,
if (i2 == this->RecvPackets.end())
{
std::list<char*> l;
l.push_back(buff);
l.emplace_back(buff);
this->RecvPackets[k] = l;
}
else
{
i2->second.push_back(buff);
i2->second.emplace_back(buff);
// The last packet came in, merge into one MemStream.
if (i2->second.size() == header.numPackets)
@ -370,7 +373,7 @@ void Messenger::ProcessReceivedBuffers(std::vector<char*>& incomingBuffers,
}
entry.second.reset();
buffers.push_back(std::move(entry));
buffers.emplace_back(std::move(entry));
this->RecvPackets.erase(i2);
}
}
@ -379,6 +382,8 @@ void Messenger::ProcessReceivedBuffers(std::vector<char*>& incomingBuffers,
}
#endif
}
}
} // namespace vtkm::filter::particleadvection
}
} // vtkm::filter::flow::internal

@ -8,13 +8,11 @@
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#ifndef vtk_m_filter_Messenger_h
#define vtk_m_filter_Messenger_h
#include <vtkm/filter/vtkm_filter_extra_export.h>
#ifndef vtk_m_filter_flow_internal_Messenger_h
#define vtk_m_filter_flow_internal_Messenger_h
#include <vtkm/Types.h>
#include <vtkm/cont/Serialization.h>
#include <vtkm/filter/flow/vtkm_filter_flow_export.h>
#include <vtkm/thirdparty/diy/diy.h>
#include <list>
@ -30,10 +28,12 @@ namespace vtkm
{
namespace filter
{
namespace particleadvection
namespace flow
{
namespace internal
{
class VTKM_FILTER_EXTRA_EXPORT Messenger
class VTKM_FILTER_FLOW_EXPORT Messenger
{
public:
VTKM_CONT Messenger(vtkmdiy::mpi::communicator& comm);
@ -106,9 +106,10 @@ protected:
static constexpr int Rank = 0;
#endif
};
}
}
} // namespace vtkm::filter::particleadvection
}
} // vtkm::filter::flow::internal
#endif
#endif // vtk_m_filter_flow_internal_Messenger_h

@ -0,0 +1,125 @@
//============================================================================
// 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_flow_internal_ParticleAdvector_h
#define vtk_m_filter_flow_internal_ParticleAdvector_h
#include <vtkm/filter/flow/internal/AdvectAlgorithm.h>
#include <vtkm/filter/flow/internal/AdvectAlgorithmThreaded.h>
#include <vtkm/filter/flow/internal/BoundsMap.h>
#include <vtkm/filter/flow/internal/DataSetIntegrator.h>
namespace vtkm
{
namespace filter
{
namespace flow
{
namespace internal
{
template <typename DSIType>
class ParticleAdvector
{
public:
ParticleAdvector(const vtkm::filter::flow::internal::BoundsMap& bm,
const std::vector<DSIType>& blocks,
const bool& useThreaded,
const vtkm::filter::flow::FlowResultType& parType)
: Blocks(blocks)
, BoundsMap(bm)
, ResultType(parType)
, UseThreadedAlgorithm(useThreaded)
{
}
vtkm::cont::PartitionedDataSet Execute(vtkm::Id numSteps,
vtkm::FloatDefault stepSize,
const vtkm::cont::UnknownArrayHandle& seeds)
{
using ParticleTypeList = vtkm::List<vtkm::Particle, vtkm::ChargedParticle>;
vtkm::cont::PartitionedDataSet result;
seeds.CastAndCallForTypes<ParticleTypeList, VTKM_DEFAULT_STORAGE_LIST>(
[&](const auto& concreteSeeds) {
result = this->Execute(numSteps, stepSize, concreteSeeds);
});
return result;
}
private:
template <typename AlgorithmType, typename ParticleType>
vtkm::cont::PartitionedDataSet RunAlgo(vtkm::Id numSteps,
vtkm::FloatDefault stepSize,
const vtkm::cont::ArrayHandle<ParticleType>& seeds)
{
AlgorithmType algo(this->BoundsMap, this->Blocks);
algo.Execute(numSteps, stepSize, seeds);
return algo.GetOutput();
}
template <typename ParticleType>
vtkm::cont::PartitionedDataSet Execute(vtkm::Id numSteps,
vtkm::FloatDefault stepSize,
const vtkm::cont::ArrayHandle<ParticleType>& seeds)
{
if (!this->UseThreadedAlgorithm)
{
if (this->ResultType == vtkm::filter::flow::FlowResultType::PARTICLE_ADVECT_TYPE)
{
using AlgorithmType = vtkm::filter::flow::internal::
AdvectAlgorithm<DSIType, vtkm::worklet::flow::ParticleAdvectionResult, ParticleType>;
return this->RunAlgo<AlgorithmType, ParticleType>(numSteps, stepSize, seeds);
}
else
{
using AlgorithmType = vtkm::filter::flow::internal::
AdvectAlgorithm<DSIType, vtkm::worklet::flow::StreamlineResult, ParticleType>;
return this->RunAlgo<AlgorithmType, ParticleType>(numSteps, stepSize, seeds);
}
}
else
{
if (this->ResultType == vtkm::filter::flow::FlowResultType::PARTICLE_ADVECT_TYPE)
{
using AlgorithmType = vtkm::filter::flow::internal::AdvectAlgorithmThreaded<
DSIType,
vtkm::worklet::flow::ParticleAdvectionResult,
ParticleType>;
return this->RunAlgo<AlgorithmType, ParticleType>(numSteps, stepSize, seeds);
}
else
{
using AlgorithmType = vtkm::filter::flow::internal::
AdvectAlgorithmThreaded<DSIType, vtkm::worklet::flow::StreamlineResult, ParticleType>;
return this->RunAlgo<AlgorithmType, ParticleType>(numSteps, stepSize, seeds);
}
}
}
std::vector<DSIType> Blocks;
vtkm::filter::flow::internal::BoundsMap BoundsMap;
FlowResultType ResultType;
bool UseThreadedAlgorithm;
};
}
}
}
} //vtkm::filter::flow::internal
#endif //vtk_m_filter_flow_internal_ParticleAdvector_h

@ -0,0 +1,368 @@
//============================================================================
// 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_flow_internal_ParticleMessenger_h
#define vtk_m_filter_flow_internal_ParticleMessenger_h
#include <vtkm/Particle.h>
#include <vtkm/filter/flow/internal/BoundsMap.h>
#include <vtkm/filter/flow/internal/Messenger.h>
#include <vtkm/filter/flow/vtkm_filter_flow_export.h>
#include <list>
#include <map>
#include <set>
#include <vector>
namespace vtkm
{
namespace filter
{
namespace flow
{
namespace internal
{
template <typename ParticleType>
class VTKM_FILTER_FLOW_EXPORT ParticleMessenger : public vtkm::filter::flow::internal::Messenger
{
//sendRank, message
using MsgCommType = std::pair<int, std::vector<int>>;
//particle + blockIDs.
using ParticleCommType = std::pair<ParticleType, std::vector<vtkm::Id>>;
//sendRank, vector of ParticleCommType.
using ParticleRecvCommType = std::pair<int, std::vector<ParticleCommType>>;
public:
VTKM_CONT ParticleMessenger(vtkmdiy::mpi::communicator& comm,
const vtkm::filter::flow::internal::BoundsMap& bm,
int msgSz = 1,
int numParticles = 128,
int numBlockIds = 2);
VTKM_CONT ~ParticleMessenger() {}
VTKM_CONT void Exchange(const std::vector<ParticleType>& outData,
const std::unordered_map<vtkm::Id, std::vector<vtkm::Id>>& outBlockIDsMap,
vtkm::Id numLocalTerm,
std::vector<ParticleType>& inData,
std::unordered_map<vtkm::Id, std::vector<vtkm::Id>>& inDataBlockIDsMap,
vtkm::Id& numTerminateMessages,
bool blockAndWait = false);
protected:
#ifdef VTKM_ENABLE_MPI
static constexpr int MSG_TERMINATE = 1;
enum { MESSAGE_TAG = 0x42000, PARTICLE_TAG = 0x42001 };
VTKM_CONT void RegisterMessages(int msgSz, int nParticles, int numBlockIds);
// Send/Recv particles
VTKM_CONT
template <typename P,
template <typename, typename>
class Container,
typename Allocator = std::allocator<P>>
inline void SendParticles(int dst, const Container<P, Allocator>& c);
VTKM_CONT
template <typename P,
template <typename, typename>
class Container,
typename Allocator = std::allocator<P>>
inline void SendParticles(const std::unordered_map<int, Container<P, Allocator>>& m);
// Send/Recv messages.
VTKM_CONT void SendMsg(int dst, const std::vector<int>& msg);
VTKM_CONT void SendAllMsg(const std::vector<int>& msg);
VTKM_CONT bool RecvMsg(std::vector<MsgCommType>& msgs) { return RecvAny(&msgs, NULL, false); }
// Send/Recv datasets.
VTKM_CONT bool RecvAny(std::vector<MsgCommType>* msgs,
std::vector<ParticleRecvCommType>* recvParticles,
bool blockAndWait);
const vtkm::filter::flow::internal::BoundsMap& BoundsMap;
#endif
VTKM_CONT void SerialExchange(
const std::vector<ParticleType>& outData,
const std::unordered_map<vtkm::Id, std::vector<vtkm::Id>>& outBlockIDsMap,
vtkm::Id numLocalTerm,
std::vector<ParticleType>& inData,
std::unordered_map<vtkm::Id, std::vector<vtkm::Id>>& inDataBlockIDsMap,
bool blockAndWait) const;
static std::size_t CalcParticleBufferSize(std::size_t nParticles, std::size_t numBlockIds = 2);
};
//methods
VTKM_CONT
template <typename ParticleType>
ParticleMessenger<ParticleType>::ParticleMessenger(
vtkmdiy::mpi::communicator& comm,
const vtkm::filter::flow::internal::BoundsMap& boundsMap,
int msgSz,
int numParticles,
int numBlockIds)
: Messenger(comm)
#ifdef VTKM_ENABLE_MPI
, BoundsMap(boundsMap)
#endif
{
#ifdef VTKM_ENABLE_MPI
this->RegisterMessages(msgSz, numParticles, numBlockIds);
#else
(void)(boundsMap);
(void)(msgSz);
(void)(numParticles);
(void)(numBlockIds);
#endif
}
template <typename ParticleType>
std::size_t ParticleMessenger<ParticleType>::CalcParticleBufferSize(std::size_t nParticles,
std::size_t nBlockIds)
{
ParticleType pTmp;
std::size_t pSize = ParticleType::Sizeof();
#ifndef NDEBUG
vtkmdiy::MemoryBuffer buff;
ParticleType p;
vtkmdiy::save(buff, p);
//Make sure the buffer size is correct.
//If this fires, then the size of the class has changed.
VTKM_ASSERT(pSize == buff.size());
#endif
return
// rank
sizeof(int)
//std::vector<ParticleType> p;
//p.size()
+ sizeof(std::size_t)
//nParticles of ParticleType
+ nParticles * pSize
// std::vector<vtkm::Id> blockIDs for each particle.
// blockIDs.size() for each particle
+ nParticles * sizeof(std::size_t)
// nBlockIDs of vtkm::Id for each particle.
+ nParticles * nBlockIds * sizeof(vtkm::Id);
}
VTKM_CONT
template <typename ParticleType>
void ParticleMessenger<ParticleType>::SerialExchange(
const std::vector<ParticleType>& outData,
const std::unordered_map<vtkm::Id, std::vector<vtkm::Id>>& outBlockIDsMap,
vtkm::Id vtkmNotUsed(numLocalTerm),
std::vector<ParticleType>& inData,
std::unordered_map<vtkm::Id, std::vector<vtkm::Id>>& inDataBlockIDsMap,
bool vtkmNotUsed(blockAndWait)) const
{
for (auto& p : outData)
{
const auto& bids = outBlockIDsMap.find(p.ID)->second;
inData.emplace_back(p);
inDataBlockIDsMap[p.ID] = bids;
}
}
VTKM_CONT
template <typename ParticleType>
void ParticleMessenger<ParticleType>::Exchange(
const std::vector<ParticleType>& outData,
const std::unordered_map<vtkm::Id, std::vector<vtkm::Id>>& outBlockIDsMap,
vtkm::Id numLocalTerm,
std::vector<ParticleType>& inData,
std::unordered_map<vtkm::Id, std::vector<vtkm::Id>>& inDataBlockIDsMap,
vtkm::Id& numTerminateMessages,
bool blockAndWait)
{
numTerminateMessages = 0;
inDataBlockIDsMap.clear();
if (this->GetNumRanks() == 1)
return this->SerialExchange(
outData, 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)
{
const auto& bids = outBlockIDsMap.find(p.ID)->second;
int dstRank = this->BoundsMap.FindRank(bids[0]);
sendData[dstRank].emplace_back(std::make_pair(p, bids));
}
//Do all the sends first.
if (numLocalTerm > 0)
this->SendAllMsg({ MSG_TERMINATE, static_cast<int>(numLocalTerm) });
this->SendParticles(sendData);
this->CheckPendingSendRequests();
//Check if we have anything coming in.
std::vector<ParticleRecvCommType> particleData;
std::vector<MsgCommType> msgData;
if (RecvAny(&msgData, &particleData, blockAndWait))
{
for (const auto& it : particleData)
for (const auto& v : it.second)
{
const auto& p = v.first;
const auto& bids = v.second;
inData.emplace_back(p);
inDataBlockIDsMap[p.ID] = bids;
}
for (const auto& m : msgData)
{
if (m.second[0] == MSG_TERMINATE)
numTerminateMessages += static_cast<vtkm::Id>(m.second[1]);
}
}
#endif
}
#ifdef VTKM_ENABLE_MPI
VTKM_CONT
template <typename ParticleType>
void ParticleMessenger<ParticleType>::RegisterMessages(int msgSz, int nParticles, int numBlockIds)
{
//Determine buffer size for msg and particle tags.
std::size_t messageBuffSz = CalcMessageBufferSize(msgSz + 1);
std::size_t particleBuffSz = CalcParticleBufferSize(nParticles, numBlockIds);
int numRecvs = std::min(64, this->GetNumRanks() - 1);
this->RegisterTag(ParticleMessenger::MESSAGE_TAG, numRecvs, messageBuffSz);
this->RegisterTag(ParticleMessenger::PARTICLE_TAG, numRecvs, particleBuffSz);
this->InitializeBuffers();
}
VTKM_CONT
template <typename ParticleType>
void ParticleMessenger<ParticleType>::SendMsg(int dst, const std::vector<int>& msg)
{
vtkmdiy::MemoryBuffer buff;
//Write data.
vtkmdiy::save(buff, this->GetRank());
vtkmdiy::save(buff, msg);
this->SendData(dst, ParticleMessenger::MESSAGE_TAG, buff);
}
VTKM_CONT
template <typename ParticleType>
void ParticleMessenger<ParticleType>::SendAllMsg(const std::vector<int>& msg)
{
for (int i = 0; i < this->GetNumRanks(); i++)
if (i != this->GetRank())
this->SendMsg(i, msg);
}
VTKM_CONT
template <typename ParticleType>
bool ParticleMessenger<ParticleType>::RecvAny(std::vector<MsgCommType>* msgs,
std::vector<ParticleRecvCommType>* recvParticles,
bool blockAndWait)
{
std::set<int> tags;
if (msgs)
{
tags.insert(ParticleMessenger::MESSAGE_TAG);
msgs->resize(0);
}
if (recvParticles)
{
tags.insert(ParticleMessenger::PARTICLE_TAG);
recvParticles->resize(0);
}
if (tags.empty())
return false;
std::vector<std::pair<int, vtkmdiy::MemoryBuffer>> buffers;
if (!this->RecvData(tags, buffers, blockAndWait))
return false;
for (auto& buff : buffers)
{
if (buff.first == ParticleMessenger::MESSAGE_TAG)
{
int sendRank;
std::vector<int> m;
vtkmdiy::load(buff.second, sendRank);
vtkmdiy::load(buff.second, m);
msgs->emplace_back(std::make_pair(sendRank, m));
}
else if (buff.first == ParticleMessenger::PARTICLE_TAG)
{
int sendRank;
std::vector<ParticleCommType> particles;
vtkmdiy::load(buff.second, sendRank);
vtkmdiy::load(buff.second, particles);
recvParticles->emplace_back(std::make_pair(sendRank, particles));
}
}
return true;
}
VTKM_CONT
template <typename ParticleType>
template <typename P, template <typename, typename> class Container, typename Allocator>
inline void ParticleMessenger<ParticleType>::SendParticles(int dst,
const Container<P, Allocator>& c)
{
if (dst == this->GetRank())
{
VTKM_LOG_S(vtkm::cont::LogLevel::Error, "Error. Sending a particle to yourself.");
return;
}
if (c.empty())
return;
vtkmdiy::MemoryBuffer bb;
vtkmdiy::save(bb, this->GetRank());
vtkmdiy::save(bb, c);
this->SendData(dst, ParticleMessenger::PARTICLE_TAG, bb);
}
VTKM_CONT
template <typename ParticleType>
template <typename P, template <typename, typename> class Container, typename Allocator>
inline void ParticleMessenger<ParticleType>::SendParticles(
const std::unordered_map<int, Container<P, Allocator>>& m)
{
for (const auto& mit : m)
if (!mit.second.empty())
this->SendParticles(mit.first, mit.second);
}
#endif
}
}
}
} // vtkm::filter::flow::internal
#endif // vtk_m_filter_flow_internal_ParticleMessenger_h

@ -0,0 +1,55 @@
##============================================================================
## Copyright (c) Kitware, Inc.
## All rights reserved.
## See LICENSE.txt for details.
##
## This software is distributed WITHOUT ANY WARRANTY; without even
## the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
## PURPOSE. See the above copyright notice for more information.
##============================================================================
set(filter_unit_tests
UnitTestStreamlineFilter.cxx
UnitTestStreamSurfaceFilter.cxx
)
set(worklet_unit_tests
UnitTestWorkletParticleAdvection.cxx
UnitTestWorkletTemporalAdvection.cxx
UnitTestStreamSurfaceWorklet.cxx
)
set(libraries
vtkm_filter_mesh_info
vtkm_filter_flow
vtkm_io
)
if (VTKm_ENABLE_RENDERING)
list(APPEND libraries vtkm_filter_geometry_refinement vtkm_rendering vtkm_rendering_testing)
list(APPEND filter_unit_tests
RenderTestStreamline.cxx
)
endif()
vtkm_unit_tests(
SOURCES ${filter_unit_tests}
DEVICE_SOURCES ${worklet_unit_tests}
LIBRARIES ${libraries}
USE_VTKM_JOB_POOL
)
#add distributed tests i.e.test to run with MPI
#if MPI is enabled.
if (VTKm_ENABLE_MPI)
set(mpi_unit_tests
UnitTestParticleMessengerMPI.cxx
UnitTestStreamlineFilterMPI.cxx
)
vtkm_unit_tests(
MPI DEVICE_SOURCES ${mpi_unit_tests}
LIBRARIES ${libraries}
ALL_BACKENDS
USE_VTKM_JOB_POOL
)
endif()

@ -0,0 +1,74 @@
//============================================================================
// 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/Math.h>
#include <vtkm/Particle.h>
#include <vtkm/cont/ArrayHandleConstant.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/io/VTKDataSetReader.h>
#include <vtkm/filter/flow/Streamline.h>
#include <vtkm/filter/geometry_refinement/Tube.h>
#include <vtkm/rendering/testing/RenderTest.h>
#include <vtkm/rendering/testing/Testing.h>
namespace
{
void TestStreamline()
{
std::cout << "Generate Image for Streamline filter" << std::endl;
auto pathname = vtkm::cont::testing::Testing::DataPath("uniform/StreamlineTestDataSet.vtk");
vtkm::io::VTKDataSetReader reader(pathname);
vtkm::cont::DataSet dataSet = reader.ReadDataSet();
vtkm::cont::ArrayHandle<vtkm::Particle> seedArray =
vtkm::cont::make_ArrayHandle({ vtkm::Particle(vtkm::Vec3f(.2f, 1.0f, .2f), 0),
vtkm::Particle(vtkm::Vec3f(.2f, 2.0f, .2f), 1),
vtkm::Particle(vtkm::Vec3f(.2f, 3.0f, .2f), 2) });
vtkm::filter::flow::Streamline streamline;
streamline.SetStepSize(0.1f);
streamline.SetNumberOfSteps(20);
streamline.SetSeeds(seedArray);
streamline.SetActiveField("vector");
auto result = streamline.Execute(dataSet);
// Some sort of color map is needed when rendering the coordinates of a dataset
// so create a zeroed array for the coordinates.
std::vector<vtkm::FloatDefault> colorMap(static_cast<std::vector<vtkm::FloatDefault>::size_type>(
result.GetCoordinateSystem().GetNumberOfPoints()));
for (std::vector<vtkm::FloatDefault>::size_type i = 0; i < colorMap.size(); i++)
{
colorMap[i] = static_cast<vtkm::FloatDefault>(i);
}
result.AddPointField("pointvar", colorMap);
// The streamline by itself doesn't generate renderable geometry, so surround the
// streamlines in tubes.
vtkm::filter::geometry_refinement::Tube tube;
tube.SetCapping(true);
tube.SetNumberOfSides(3);
tube.SetRadius(static_cast<vtkm::FloatDefault>(0.2));
result = tube.Execute(result);
result.PrintSummary(std::cout);
vtkm::rendering::testing::RenderTestOptions testOptions;
testOptions.ColorTable = vtkm::cont::ColorTable::Preset::Inferno;
testOptions.EnableAnnotations = false;
vtkm::rendering::testing::RenderTest(result, "pointvar", "filter/streamline.png", testOptions);
}
} // namespace
int RenderTestStreamline(int argc, char* argv[])
{
return vtkm::cont::testing::Testing::Run(TestStreamline, argc, argv);
}

@ -10,7 +10,7 @@
#include <vtkm/cont/Serialization.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/filter/particleadvection/ParticleMessenger.h>
#include <vtkm/filter/flow/internal/ParticleMessenger.h>
#include <vtkm/thirdparty/diy/diy.h>
@ -23,11 +23,11 @@ using MCommType = std::pair<int, std::vector<int>>;
using PCommType = std::pair<vtkm::Particle, std::vector<vtkm::Id>>;
using PRecvCommType = std::pair<int, std::vector<PCommType>>;
class TestMessenger : public vtkm::filter::particleadvection::ParticleMessenger
class TestMessenger : public vtkm::filter::flow::internal::ParticleMessenger<vtkm::Particle>
{
public:
TestMessenger(vtkmdiy::mpi::communicator& comm,
const vtkm::filter::particleadvection::BoundsMap& bm,
const vtkm::filter::flow::internal::BoundsMap& bm,
int msgSz = 1,
int numParticles = 1,
int numBlockIds = 1)
@ -120,7 +120,7 @@ void TestParticleMessenger()
//Only works for 2 or more ranks.
if (comm.size() == 1)
return;
vtkm::filter::particleadvection::BoundsMap boundsMap;
vtkm::filter::flow::internal::BoundsMap boundsMap;
int maxMsgSz = 100;
int maxNumParticles = 128;
@ -244,7 +244,7 @@ void TestBufferSizes()
{
//Make sure the buffer sizes are correct.
auto comm = vtkm::cont::EnvironmentTracker::GetCommunicator();
vtkm::filter::particleadvection::BoundsMap boundsMap;
vtkm::filter::flow::internal::BoundsMap boundsMap;
std::vector<int> mSzs = { 1, 2, 3, 4, 5 };
std::vector<int> numPs = { 1, 2, 3, 4, 5 };

@ -8,9 +8,10 @@
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#include <vtkm/Particle.h>
#include <vtkm/cont/DataSetBuilderUniform.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/filter/StreamSurface.h>
#include <vtkm/filter/flow/StreamSurface.h>
namespace
{
@ -44,7 +45,7 @@ void TestStreamSurface()
vtkm::Particle(vtkm::Vec3f(.1f, 3.0f, .3f), 2),
vtkm::Particle(vtkm::Vec3f(.1f, 3.5f, .2f), 3) });
vtkm::filter::StreamSurface streamSrf;
vtkm::filter::flow::StreamSurface streamSrf;
streamSrf.SetStepSize(0.1f);
streamSrf.SetNumberOfSteps(20);

@ -8,11 +8,11 @@
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#include <vtkm/cont/Algorithm.h>
#include <vtkm/cont/DataSetBuilderExplicit.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/filter/flow/worklet/StreamSurface.h>
#include <vtkm/io/VTKDataSetWriter.h>
#include <vtkm/worklet/DispatcherMapField.h>
#include <vtkm/worklet/StreamSurface.h>
namespace
{
@ -53,7 +53,7 @@ void TestSameNumPolylines()
dsb.AddCell(vtkm::CELL_SHAPE_POLY_LINE, ids);
vtkm::cont::DataSet ds = dsb.Create();
vtkm::worklet::StreamSurface streamSurfaceWorklet;
vtkm::worklet::flow::StreamSurface streamSurfaceWorklet;
vtkm::cont::ArrayHandle<vtkm::Vec3f> newPoints;
vtkm::cont::CellSetSingleType<> newCells;
streamSurfaceWorklet.Run(ds.GetCoordinateSystem(0), ds.GetCellSet(), newPoints, newCells);
@ -111,7 +111,7 @@ void TestUnequalNumPolylines(int unequalType)
dsb.AddCell(vtkm::CELL_SHAPE_POLY_LINE, ids);
vtkm::cont::DataSet ds = dsb.Create();
vtkm::worklet::StreamSurface streamSurfaceWorklet;
vtkm::worklet::flow::StreamSurface streamSurfaceWorklet;
vtkm::cont::ArrayHandle<vtkm::Vec3f> newPoints;
vtkm::cont::CellSetSingleType<> newCells;
streamSurfaceWorklet.Run(ds.GetCoordinateSystem(0), ds.GetCellSet(), newPoints, newCells);
@ -124,7 +124,7 @@ void TestUnequalNumPolylines(int unequalType)
"Wrong number of cells in StreamSurface worklet");
}
void TestStreamSurface()
void TestStreamSurfaceWorklet()
{
std::cout << "Testing Stream Surface Worklet" << std::endl;
TestSameNumPolylines();
@ -134,7 +134,7 @@ void TestStreamSurface()
}
}
int UnitTestStreamSurface(int argc, char* argv[])
int UnitTestStreamSurfaceWorklet(int argc, char* argv[])
{
return vtkm::cont::testing::Testing::Run(TestStreamSurface, argc, argv);
return vtkm::cont::testing::Testing::Run(TestStreamSurfaceWorklet, argc, argv);
}

@ -8,12 +8,15 @@
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#include <vtkm/CellClassification.h>
#include <vtkm/Particle.h>
#include <vtkm/cont/ArrayCopy.h>
#include <vtkm/cont/DataSetBuilderUniform.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/filter/ParticleAdvection.h>
#include <vtkm/filter/PathParticle.h>
#include <vtkm/filter/Pathline.h>
#include <vtkm/filter/Streamline.h>
#include <vtkm/filter/flow/ParticleAdvection.h>
#include <vtkm/filter/flow/PathParticle.h>
#include <vtkm/filter/flow/Pathline.h>
#include <vtkm/filter/flow/Streamline.h>
#include <vtkm/io/VTKDataSetReader.h>
#include <vtkm/worklet/testing/GenerateTestDataSets.h>
@ -63,7 +66,7 @@ void TestStreamline()
vtkm::Particle(vtkm::Vec3f(.2f, 2.0f, .2f), 1),
vtkm::Particle(vtkm::Vec3f(.2f, 3.0f, .2f), 2) });
vtkm::filter::Streamline streamline;
vtkm::filter::flow::Streamline streamline;
streamline.SetStepSize(0.1f);
streamline.SetNumberOfSteps(20);
@ -122,7 +125,7 @@ void TestPathline()
vtkm::Id numExpectedPoints;
if (fType == 0)
{
vtkm::filter::Pathline filt;
vtkm::filter::flow::Pathline filt;
filt.SetActiveField(var);
filt.SetStepSize(stepSize);
filt.SetNumberOfSteps(numSteps);
@ -135,7 +138,7 @@ void TestPathline()
}
else
{
vtkm::filter::PathParticle filt;
vtkm::filter::flow::PathParticle filt;
filt.SetActiveField(var);
filt.SetStepSize(stepSize);
filt.SetNumberOfSteps(numSteps);
@ -192,7 +195,7 @@ void TestAMRStreamline(bool useSL)
ghosts[idx] = vtkm::CellClassification::Normal;
idx++;
}
dsOuter.AddCellField("vtkmGhostCells", ghosts);
dsOuter.AddGhostCellField(vtkm::cont::make_ArrayHandle(ghosts, vtkm::CopyFlag::On));
//Create a partitioned dataset with 1 inner and 1 outer.
vtkm::cont::PartitionedDataSet pds;
@ -221,7 +224,7 @@ void TestAMRStreamline(bool useSL)
if (useSL)
{
vtkm::filter::Streamline filter;
vtkm::filter::flow::Streamline filter;
filter.SetStepSize(0.1f);
filter.SetNumberOfSteps(100000);
filter.SetSeeds(seedArray);
@ -293,7 +296,7 @@ void TestAMRStreamline(bool useSL)
}
else
{
vtkm::filter::ParticleAdvection filter;
vtkm::filter::flow::ParticleAdvection filter;
filter.SetStepSize(0.1f);
filter.SetNumberOfSteps(100000);
filter.SetSeeds(seedArray);
@ -375,7 +378,7 @@ void TestPartitionedDataSet(vtkm::Id num, bool useGhost, FilterType fType)
vtkm::cont::PartitionedDataSet out;
if (fType == FilterType::STREAMLINE)
{
vtkm::filter::Streamline streamline;
vtkm::filter::flow::Streamline streamline;
streamline.SetStepSize(0.1f);
streamline.SetNumberOfSteps(100000);
streamline.SetSeeds(seedArray);
@ -388,7 +391,7 @@ void TestPartitionedDataSet(vtkm::Id num, bool useGhost, FilterType fType)
auto pds2 = allPDs2[idx];
AddVectorFields(pds2, fieldName, vecX);
vtkm::filter::Pathline pathline;
vtkm::filter::flow::Pathline pathline;
pathline.SetPreviousTime(0);
pathline.SetNextTime(1000);
pathline.SetNextDataSet(pds2);
@ -439,7 +442,7 @@ void TestPartitionedDataSet(vtkm::Id num, bool useGhost, FilterType fType)
vtkm::cont::PartitionedDataSet out;
if (fType == FilterType::PARTICLE_ADVECTION)
{
vtkm::filter::ParticleAdvection particleAdvection;
vtkm::filter::flow::ParticleAdvection particleAdvection;
particleAdvection.SetStepSize(0.1f);
particleAdvection.SetNumberOfSteps(100000);
@ -453,7 +456,7 @@ void TestPartitionedDataSet(vtkm::Id num, bool useGhost, FilterType fType)
auto pds2 = allPDs2[idx];
AddVectorFields(pds2, fieldName, vecX);
vtkm::filter::PathParticle pathParticle;
vtkm::filter::flow::PathParticle pathParticle;
pathParticle.SetPreviousTime(0);
pathParticle.SetNextTime(1000);
pathParticle.SetNextDataSet(pds2);
@ -465,6 +468,7 @@ void TestPartitionedDataSet(vtkm::Id num, bool useGhost, FilterType fType)
out = pathParticle.Execute(pds);
}
VTKM_TEST_ASSERT(out.GetNumberOfPartitions() == 1, "Wrong number of partitions in output");
auto ds = out.GetPartition(0);
//Validate the result is correct.
@ -542,7 +546,7 @@ void TestStreamlineFile(const std::string& fname,
vtkm::cont::DataSet output;
if (useSL)
{
vtkm::filter::Streamline streamline;
vtkm::filter::flow::Streamline streamline;
streamline.SetStepSize(stepSize);
streamline.SetNumberOfSteps(maxSteps);
streamline.SetSeeds(seedArray);
@ -551,7 +555,7 @@ void TestStreamlineFile(const std::string& fname,
}
else
{
vtkm::filter::ParticleAdvection particleAdvection;
vtkm::filter::flow::ParticleAdvection particleAdvection;
particleAdvection.SetStepSize(stepSize);
particleAdvection.SetNumberOfSteps(maxSteps);
particleAdvection.SetSeeds(seedArray);
@ -584,13 +588,14 @@ void TestStreamlineFilters()
FilterType::STREAMLINE,
FilterType::PATHLINE,
FilterType::PATH_PARTICLE };
//fTypes = {FilterType::PARTICLE_ADVECTION,FilterType::STREAMLINE};
//fTypes = {FilterType::PATHLINE};
for (int n = 1; n < 3; n++)
{
for (auto useGhost : flags)
for (auto ft : fTypes)
{
TestPartitionedDataSet(n, useGhost, ft);
}
}
TestStreamline();

@ -8,10 +8,12 @@
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#include <vtkm/CellClassification.h>
#include <vtkm/Particle.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/filter/ParticleAdvection.h>
#include <vtkm/filter/Pathline.h>
#include <vtkm/filter/Streamline.h>
#include <vtkm/filter/flow/ParticleAdvection.h>
#include <vtkm/filter/flow/Pathline.h>
#include <vtkm/filter/flow/Streamline.h>
#include <vtkm/thirdparty/diy/diy.h>
#include <vtkm/worklet/testing/GenerateTestDataSets.h>
@ -111,7 +113,7 @@ void TestAMRStreamline(FilterType fType, bool useThreaded)
ghosts[idx] = vtkm::CellClassification::Normal;
idx++;
}
dsOuter.AddCellField("vtkmGhostCells", ghosts);
dsOuter.AddGhostCellField(vtkm::cont::make_ArrayHandle(ghosts, vtkm::CopyFlag::On));
//Create a partitioned dataset with 1 inner and 1 outer.
vtkm::cont::PartitionedDataSet pds;
@ -149,13 +151,13 @@ void TestAMRStreamline(FilterType fType, bool useThreaded)
if (fType == STREAMLINE)
{
vtkm::filter::Streamline streamline;
vtkm::filter::flow::Streamline streamline;
SetFilter(streamline, stepSize, numSteps, fieldName, seedArray, useThreaded);
out = streamline.Execute(pds);
}
else if (fType == PATHLINE)
{
vtkm::filter::Pathline pathline;
vtkm::filter::flow::Pathline pathline;
SetFilter(pathline, stepSize, numSteps, fieldName, seedArray, useThreaded);
//Create timestep 2
auto pds2 = vtkm::cont::PartitionedDataSet(pds);
@ -239,7 +241,7 @@ void TestAMRStreamline(FilterType fType, bool useThreaded)
}
else if (fType == PARTICLE_ADVECTION)
{
vtkm::filter::ParticleAdvection filter;
vtkm::filter::flow::ParticleAdvection filter;
filter.SetUseThreadedAlgorithm(useThreaded);
filter.SetStepSize(0.1f);
filter.SetNumberOfSteps(100000);
@ -283,6 +285,8 @@ void ValidateOutput(const vtkm::cont::DataSet& out,
"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");
auto coords = out.GetCoordinateSystem().GetDataAsMultiplexer();
auto ptPortal = coords.ReadPortal();
@ -392,7 +396,7 @@ void TestPartitionedDataSet(vtkm::Id nPerRank, bool useGhost, FilterType fType,
if (fType == STREAMLINE)
{
vtkm::filter::Streamline streamline;
vtkm::filter::flow::Streamline streamline;
SetFilter(streamline, stepSize, numSteps, fieldName, seedArray, useThreaded);
auto out = streamline.Execute(pds);
@ -401,7 +405,7 @@ void TestPartitionedDataSet(vtkm::Id nPerRank, bool useGhost, FilterType fType,
}
else if (fType == PARTICLE_ADVECTION)
{
vtkm::filter::ParticleAdvection particleAdvection;
vtkm::filter::flow::ParticleAdvection particleAdvection;
SetFilter(particleAdvection, stepSize, numSteps, fieldName, seedArray, useThreaded);
auto out = particleAdvection.Execute(pds);
@ -419,7 +423,7 @@ void TestPartitionedDataSet(vtkm::Id nPerRank, bool useGhost, FilterType fType,
auto pds2 = allPDS2[n];
AddVectorFields(pds2, fieldName, vecX);
vtkm::filter::Pathline pathline;
vtkm::filter::flow::Pathline pathline;
SetFilter(pathline, stepSize, numSteps, fieldName, seedArray, useThreaded);
pathline.SetPreviousTime(time0);

@ -9,19 +9,20 @@
//============================================================================
#include <typeinfo>
#include <vtkm/cont/Algorithm.h>
#include <vtkm/cont/ArrayCopy.h>
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/filter/flow/worklet/EulerIntegrator.h>
#include <vtkm/filter/flow/worklet/Field.h>
#include <vtkm/filter/flow/worklet/GridEvaluators.h>
#include <vtkm/filter/flow/worklet/ParticleAdvection.h>
#include <vtkm/filter/flow/worklet/Particles.h>
#include <vtkm/filter/flow/worklet/RK4Integrator.h>
#include <vtkm/filter/flow/worklet/Stepper.h>
#include <vtkm/filter/mesh_info/GhostCellClassify.h>
#include <vtkm/io/VTKDataSetReader.h>
#include <vtkm/worklet/ParticleAdvection.h>
#include <vtkm/worklet/particleadvection/EulerIntegrator.h>
#include <vtkm/worklet/particleadvection/Field.h>
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
#include <vtkm/worklet/particleadvection/Particles.h>
#include <vtkm/worklet/particleadvection/RK4Integrator.h>
#include <vtkm/worklet/particleadvection/Stepper.h>
#include <vtkm/worklet/testing/GenerateTestDataSets.h>
#include <random>
@ -142,7 +143,7 @@ public:
template <typename EvaluatorType>
VTKM_EXEC void operator()(vtkm::Particle& pointIn,
const EvaluatorType& evaluator,
vtkm::worklet::particleadvection::GridEvaluatorStatus& status,
vtkm::worklet::flow::GridEvaluatorStatus& status,
vtkm::Vec3f& pointOut) const
{
vtkm::VecVariable<vtkm::Vec3f, 2> values;
@ -159,7 +160,7 @@ void ValidateEvaluator(const EvalType& eval,
{
using EvalTester = TestEvaluatorWorklet;
using EvalTesterDispatcher = vtkm::worklet::DispatcherMapField<EvalTester>;
using Status = vtkm::worklet::particleadvection::GridEvaluatorStatus;
using Status = vtkm::worklet::flow::GridEvaluatorStatus;
EvalTester evalTester;
EvalTesterDispatcher evalTesterDispatcher(evalTester);
vtkm::cont::ArrayHandle<vtkm::Particle> pointsHandle =
@ -192,7 +193,7 @@ public:
template <typename Particle, typename IntegratorType>
VTKM_EXEC void operator()(Particle& pointIn,
const IntegratorType integrator,
vtkm::worklet::particleadvection::IntegratorStatus& status,
vtkm::worklet::flow::IntegratorStatus& status,
vtkm::Vec3f& pointOut) const
{
vtkm::FloatDefault time = 0;
@ -211,7 +212,7 @@ void ValidateIntegrator(const IntegratorType& integrator,
{
using IntegratorTester = TestIntegratorWorklet;
using IntegratorTesterDispatcher = vtkm::worklet::DispatcherMapField<IntegratorTester>;
using Status = vtkm::worklet::particleadvection::IntegratorStatus;
using Status = vtkm::worklet::flow::IntegratorStatus;
IntegratorTesterDispatcher integratorTesterDispatcher;
auto pointsHandle = vtkm::cont::make_ArrayHandle(pointIns, vtkm::CopyFlag::Off);
vtkm::Id numPoints = pointsHandle.GetNumberOfValues();
@ -243,7 +244,7 @@ void ValidateIntegratorForBoundary(const vtkm::Bounds& bounds,
{
using IntegratorTester = TestIntegratorWorklet;
using IntegratorTesterDispatcher = vtkm::worklet::DispatcherMapField<IntegratorTester>;
using Status = vtkm::worklet::particleadvection::IntegratorStatus;
using Status = vtkm::worklet::flow::IntegratorStatus;
IntegratorTesterDispatcher integratorTesterDispatcher;
auto pointsHandle = vtkm::cont::make_ArrayHandle(pointIns, vtkm::CopyFlag::Off);
@ -269,10 +270,10 @@ void ValidateIntegratorForBoundary(const vtkm::Bounds& bounds,
void TestEvaluators()
{
using FieldHandle = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
using FieldType = vtkm::worklet::particleadvection::VelocityField<FieldHandle>;
using GridEvalType = vtkm::worklet::particleadvection::GridEvaluator<FieldType>;
using RK4Type = vtkm::worklet::particleadvection::RK4Integrator<GridEvalType>;
using Stepper = vtkm::worklet::particleadvection::Stepper<RK4Type, GridEvalType>;
using FieldType = vtkm::worklet::flow::VelocityField<FieldHandle>;
using GridEvalType = vtkm::worklet::flow::GridEvaluator<FieldType>;
using RK4Type = vtkm::worklet::flow::RK4Integrator<GridEvalType>;
using Stepper = vtkm::worklet::flow::Stepper<RK4Type, GridEvalType>;
std::vector<vtkm::Vec3f> vecs;
vtkm::FloatDefault vals[3] = { -1., 0., 1. };
@ -364,10 +365,10 @@ void TestEvaluators()
void TestGhostCellEvaluators()
{
using FieldHandle = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
using FieldType = vtkm::worklet::particleadvection::VelocityField<FieldHandle>;
using GridEvalType = vtkm::worklet::particleadvection::GridEvaluator<FieldType>;
using RK4Type = vtkm::worklet::particleadvection::RK4Integrator<GridEvalType>;
using Stepper = vtkm::worklet::particleadvection::Stepper<RK4Type, GridEvalType>;
using FieldType = vtkm::worklet::flow::VelocityField<FieldHandle>;
using GridEvalType = vtkm::worklet::flow::GridEvaluator<FieldType>;
using RK4Type = vtkm::worklet::flow::RK4Integrator<GridEvalType>;
using Stepper = vtkm::worklet::flow::Stepper<RK4Type, GridEvalType>;
constexpr vtkm::Id nX = 6;
constexpr vtkm::Id nY = 6;
@ -395,7 +396,7 @@ void TestGhostCellEvaluators()
vtkm::FloatDefault stepSize = static_cast<vtkm::FloatDefault>(0.1);
Stepper rk4(gridEval, stepSize);
vtkm::worklet::ParticleAdvection pa;
vtkm::worklet::flow::ParticleAdvection pa;
std::vector<vtkm::Particle> seeds;
//Points in a ghost cell.
seeds.push_back(vtkm::Particle(vtkm::Vec3f(.5, .5, .5), 0));
@ -426,7 +427,7 @@ void TestGhostCellEvaluators()
}
void ValidateParticleAdvectionResult(
const vtkm::worklet::ParticleAdvectionResult<vtkm::Particle>& res,
const vtkm::worklet::flow::ParticleAdvectionResult<vtkm::Particle>& res,
vtkm::Id nSeeds,
vtkm::Id maxSteps)
{
@ -446,7 +447,7 @@ void ValidateParticleAdvectionResult(
}
}
void ValidateStreamlineResult(const vtkm::worklet::StreamlineResult<vtkm::Particle>& res,
void ValidateStreamlineResult(const vtkm::worklet::flow::StreamlineResult<vtkm::Particle>& res,
vtkm::Id nSeeds,
vtkm::Id maxSteps)
{
@ -465,8 +466,8 @@ void ValidateStreamlineResult(const vtkm::worklet::StreamlineResult<vtkm::Partic
void TestIntegrators()
{
using FieldHandle = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
using FieldType = vtkm::worklet::particleadvection::VelocityField<FieldHandle>;
using GridEvalType = vtkm::worklet::particleadvection::GridEvaluator<FieldType>;
using FieldType = vtkm::worklet::flow::VelocityField<FieldHandle>;
using GridEvalType = vtkm::worklet::flow::GridEvaluator<FieldType>;
const vtkm::Id3 dims(5, 5, 5);
const vtkm::Bounds bounds(0., 1., 0., 1., .0, .1);
@ -491,20 +492,20 @@ void TestIntegrators()
std::vector<vtkm::Particle> points;
GenerateRandomParticles(points, 3, bounds);
vtkm::worklet::ParticleAdvection pa;
vtkm::worklet::ParticleAdvectionResult<vtkm::Particle> res;
vtkm::worklet::flow::ParticleAdvection pa;
vtkm::worklet::flow::ParticleAdvectionResult<vtkm::Particle> res;
{
auto seeds = vtkm::cont::make_ArrayHandle(points, vtkm::CopyFlag::On);
using IntegratorType = vtkm::worklet::particleadvection::RK4Integrator<GridEvalType>;
using Stepper = vtkm::worklet::particleadvection::Stepper<IntegratorType, GridEvalType>;
using IntegratorType = vtkm::worklet::flow::RK4Integrator<GridEvalType>;
using Stepper = vtkm::worklet::flow::Stepper<IntegratorType, GridEvalType>;
Stepper rk4(eval, stepSize);
res = pa.Run(rk4, seeds, maxSteps);
ValidateParticleAdvectionResult(res, nSeeds, maxSteps);
}
{
auto seeds = vtkm::cont::make_ArrayHandle(points, vtkm::CopyFlag::On);
using IntegratorType = vtkm::worklet::particleadvection::EulerIntegrator<GridEvalType>;
using Stepper = vtkm::worklet::particleadvection::Stepper<IntegratorType, GridEvalType>;
using IntegratorType = vtkm::worklet::flow::EulerIntegrator<GridEvalType>;
using Stepper = vtkm::worklet::flow::Stepper<IntegratorType, GridEvalType>;
Stepper euler(eval, stepSize);
res = pa.Run(euler, seeds, maxSteps);
ValidateParticleAdvectionResult(res, nSeeds, maxSteps);
@ -515,10 +516,10 @@ void TestIntegrators()
void TestParticleWorkletsWithDataSetTypes()
{
using FieldHandle = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
using FieldType = vtkm::worklet::particleadvection::VelocityField<FieldHandle>;
using GridEvalType = vtkm::worklet::particleadvection::GridEvaluator<FieldType>;
using RK4Type = vtkm::worklet::particleadvection::RK4Integrator<GridEvalType>;
using Stepper = vtkm::worklet::particleadvection::Stepper<RK4Type, GridEvalType>;
using FieldType = vtkm::worklet::flow::VelocityField<FieldHandle>;
using GridEvalType = vtkm::worklet::flow::GridEvaluator<FieldType>;
using RK4Type = vtkm::worklet::flow::RK4Integrator<GridEvalType>;
using Stepper = vtkm::worklet::flow::Stepper<RK4Type, GridEvalType>;
vtkm::FloatDefault stepSize = 0.01f;
const vtkm::Id3 dims(5, 5, 5);
@ -569,8 +570,8 @@ void TestParticleWorkletsWithDataSetTypes()
{
if (i < 2)
{
vtkm::worklet::ParticleAdvection pa;
vtkm::worklet::ParticleAdvectionResult<vtkm::Particle> res;
vtkm::worklet::flow::ParticleAdvection pa;
vtkm::worklet::flow::ParticleAdvectionResult<vtkm::Particle> res;
if (i == 0)
{
auto seeds = vtkm::cont::make_ArrayHandle(pts, vtkm::CopyFlag::On);
@ -585,8 +586,8 @@ void TestParticleWorkletsWithDataSetTypes()
}
else
{
vtkm::worklet::Streamline s;
vtkm::worklet::StreamlineResult<vtkm::Particle> res;
vtkm::worklet::flow::Streamline s;
vtkm::worklet::flow::StreamlineResult<vtkm::Particle> res;
if (i == 2)
{
auto seeds = vtkm::cont::make_ArrayHandle(pts, vtkm::CopyFlag::On);
@ -618,10 +619,10 @@ void TestParticleStatus()
auto dataSets = vtkm::worklet::testing::CreateAllDataSets(bounds, dims, false);
for (auto& ds : dataSets)
{
using FieldType = vtkm::worklet::particleadvection::VelocityField<FieldHandle>;
using GridEvalType = vtkm::worklet::particleadvection::GridEvaluator<FieldType>;
using RK4Type = vtkm::worklet::particleadvection::RK4Integrator<GridEvalType>;
using Stepper = vtkm::worklet::particleadvection::Stepper<RK4Type, GridEvalType>;
using FieldType = vtkm::worklet::flow::VelocityField<FieldHandle>;
using GridEvalType = vtkm::worklet::flow::GridEvaluator<FieldType>;
using RK4Type = vtkm::worklet::flow::RK4Integrator<GridEvalType>;
using Stepper = vtkm::worklet::flow::Stepper<RK4Type, GridEvalType>;
vtkm::Id maxSteps = 1000;
vtkm::FloatDefault stepSize = 0.01f;
@ -631,7 +632,7 @@ void TestParticleStatus()
GridEvalType eval(ds, velocities);
Stepper rk4(eval, stepSize);
vtkm::worklet::ParticleAdvection pa;
vtkm::worklet::flow::ParticleAdvection pa;
std::vector<vtkm::Particle> pts;
pts.push_back(vtkm::Particle(vtkm::Vec3f(.5, .5, .5), 0));
pts.push_back(vtkm::Particle(vtkm::Vec3f(-1, -1, -1), 1));
@ -649,10 +650,10 @@ void TestParticleStatus()
void TestWorkletsBasic()
{
using FieldHandle = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
using FieldType = vtkm::worklet::particleadvection::VelocityField<FieldHandle>;
using GridEvalType = vtkm::worklet::particleadvection::GridEvaluator<FieldType>;
using RK4Type = vtkm::worklet::particleadvection::RK4Integrator<GridEvalType>;
using Stepper = vtkm::worklet::particleadvection::Stepper<RK4Type, GridEvalType>;
using FieldType = vtkm::worklet::flow::VelocityField<FieldHandle>;
using GridEvalType = vtkm::worklet::flow::GridEvaluator<FieldType>;
using RK4Type = vtkm::worklet::flow::RK4Integrator<GridEvalType>;
using Stepper = vtkm::worklet::flow::Stepper<RK4Type, GridEvalType>;
vtkm::FloatDefault stepSize = 0.01f;
const vtkm::Id3 dims(5, 5, 5);
@ -711,8 +712,8 @@ void TestWorkletsBasic()
if (w == "particleAdvection")
{
vtkm::worklet::ParticleAdvection pa;
vtkm::worklet::ParticleAdvectionResult<vtkm::Particle> res;
vtkm::worklet::flow::ParticleAdvection pa;
vtkm::worklet::flow::ParticleAdvectionResult<vtkm::Particle> res;
res = pa.Run(rk4, seedsArray, maxSteps);
@ -735,8 +736,8 @@ void TestWorkletsBasic()
}
else if (w == "streamline")
{
vtkm::worklet::Streamline s;
vtkm::worklet::StreamlineResult<vtkm::Particle> res;
vtkm::worklet::flow::Streamline s;
vtkm::worklet::flow::StreamlineResult<vtkm::Particle> res;
res = s.Run(rk4, seedsArray, maxSteps);
@ -830,10 +831,10 @@ void TestParticleAdvectionFile(const std::string& fname,
}
using FieldHandle = vtkm::cont::ArrayHandle<vtkm::Vec3f>;
using FieldType = vtkm::worklet::particleadvection::VelocityField<FieldHandle>;
using GridEvalType = vtkm::worklet::particleadvection::GridEvaluator<FieldType>;
using RK4Type = vtkm::worklet::particleadvection::RK4Integrator<GridEvalType>;
using Stepper = vtkm::worklet::particleadvection::Stepper<RK4Type, GridEvalType>;
using FieldType = vtkm::worklet::flow::VelocityField<FieldHandle>;
using GridEvalType = vtkm::worklet::flow::GridEvaluator<FieldType>;
using RK4Type = vtkm::worklet::flow::RK4Integrator<GridEvalType>;
using Stepper = vtkm::worklet::flow::Stepper<RK4Type, GridEvalType>;
VTKM_TEST_ASSERT(ds.HasField("vec"), "Data set missing a field named 'vec'");
vtkm::cont::Field& field = ds.GetField("vec");
@ -862,16 +863,16 @@ void TestParticleAdvectionFile(const std::string& fname,
if (i == 0)
{
vtkm::worklet::ParticleAdvection pa;
vtkm::worklet::ParticleAdvectionResult<vtkm::Particle> res;
vtkm::worklet::flow::ParticleAdvection pa;
vtkm::worklet::flow::ParticleAdvectionResult<vtkm::Particle> res;
res = pa.Run(rk4, seedArray, maxSteps);
ValidateResult(res, maxSteps, endPts);
}
else if (i == 1)
{
vtkm::worklet::Streamline s;
vtkm::worklet::StreamlineResult<vtkm::Particle> res;
vtkm::worklet::flow::Streamline s;
vtkm::worklet::flow::StreamlineResult<vtkm::Particle> res;
res = s.Run(rk4, seedArray, maxSteps);
ValidateResult(res, maxSteps, endPts);
@ -917,7 +918,7 @@ void TestParticleAdvection()
TestParticleAdvectionFile(fishFile, fishPts, fishStep, 100, fishEndPts);
}
int UnitTestParticleAdvection(int argc, char* argv[])
int UnitTestWorkletParticleAdvection(int argc, char* argv[])
{
return vtkm::cont::testing::Testing::Run(TestParticleAdvection, argc, argv);
}

@ -10,6 +10,7 @@
#include <typeinfo>
#include <vtkm/VecVariable.h>
#include <vtkm/cont/Algorithm.h>
#include <vtkm/cont/ArrayCopy.h>
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/DataSet.h>
@ -17,11 +18,11 @@
#include <vtkm/cont/DataSetBuilderRectilinear.h>
#include <vtkm/cont/DataSetBuilderUniform.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/worklet/ParticleAdvection.h>
#include <vtkm/worklet/particleadvection/Field.h>
#include <vtkm/worklet/particleadvection/Particles.h>
#include <vtkm/worklet/particleadvection/Stepper.h>
#include <vtkm/worklet/particleadvection/TemporalGridEvaluators.h>
#include <vtkm/filter/flow/worklet/Field.h>
#include <vtkm/filter/flow/worklet/ParticleAdvection.h>
#include <vtkm/filter/flow/worklet/Particles.h>
#include <vtkm/filter/flow/worklet/Stepper.h>
#include <vtkm/filter/flow/worklet/TemporalGridEvaluators.h>
template <typename ScalarType>
vtkm::cont::DataSet CreateUniformDataSet(const vtkm::Bounds& bounds, const vtkm::Id3& dims)
@ -52,7 +53,7 @@ public:
template <typename EvaluatorType>
VTKM_EXEC void operator()(vtkm::Particle& pointIn,
const EvaluatorType& evaluator,
vtkm::worklet::particleadvection::GridEvaluatorStatus& status,
vtkm::worklet::flow::GridEvaluatorStatus& status,
vtkm::Vec3f& pointOut) const
{
vtkm::VecVariable<vtkm::Vec3f, 2> values;
@ -70,7 +71,7 @@ void ValidateEvaluator(const EvalType& eval,
{
using EvalTester = TestEvaluatorWorklet;
using EvalTesterDispatcher = vtkm::worklet::DispatcherMapField<EvalTester>;
using Status = vtkm::worklet::particleadvection::GridEvaluatorStatus;
using Status = vtkm::worklet::flow::GridEvaluatorStatus;
EvalTester evalTester;
EvalTesterDispatcher evalTesterDispatcher(evalTester);
@ -149,9 +150,9 @@ void TestTemporalEvaluators()
using ScalarType = vtkm::FloatDefault;
using PointType = vtkm::Vec<ScalarType, 3>;
using FieldHandle = vtkm::cont::ArrayHandle<PointType>;
using FieldType = vtkm::worklet::particleadvection::VelocityField<FieldHandle>;
using EvalType = vtkm::worklet::particleadvection::GridEvaluator<FieldType>;
using TemporalEvalType = vtkm::worklet::particleadvection::TemporalGridEvaluator<FieldType>;
using FieldType = vtkm::worklet::flow::VelocityField<FieldHandle>;
using EvalType = vtkm::worklet::flow::GridEvaluator<FieldType>;
using TemporalEvalType = vtkm::worklet::flow::TemporalGridEvaluator<FieldType>;
// Create Datasets
vtkm::Id3 dims(5, 5, 5);
@ -189,7 +190,7 @@ void TestTemporalAdvection()
TestTemporalEvaluators();
}
int UnitTestTemporalAdvection(int argc, char* argv[])
int UnitTestWorkletTemporalAdvection(int argc, char* argv[])
{
return vtkm::cont::testing::Testing::Run(TestTemporalAdvection, argc, argv);
}

@ -9,6 +9,7 @@
##============================================================================
set(headers
ParticleAdvection.h
CellInterpolationHelper.h
EulerIntegrator.h
Field.h
@ -20,6 +21,7 @@ set(headers
ParticleAdvectionWorklets.h
RK4Integrator.h
TemporalGridEvaluators.h
StreamSurface.h
)
#-----------------------------------------------------------------------------

@ -8,8 +8,8 @@
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#ifndef vtk_m_worklet_particle_advection_cell_interpolation_helper
#define vtk_m_worklet_particle_advection_cell_interpolation_helper
#ifndef vtk_m_filter_flow_worklet_CellInterpolationHelper_h
#define vtk_m_filter_flow_worklet_CellInterpolationHelper_h
#include <vtkm/CellShape.h>
#include <vtkm/Types.h>
@ -282,4 +282,4 @@ private:
} //namespace cont
} //namespace vtkm
#endif //vtk_m_worklet_particle_advection_cell_interpolation_helper
#endif //vtk_m_filter_flow_worklet_CellInterpolationHelper_h

@ -10,14 +10,17 @@
//
//=============================================================================
#ifndef vtk_m_worklet_particleadvection_EulerIntegrator_h
#define vtk_m_worklet_particleadvection_EulerIntegrator_h
#ifndef vtk_m_filter_flow_worklet_EulerIntegrator_h
#define vtk_m_filter_flow_worklet_EulerIntegrator_h
#include <vtkm/filter/flow/worklet/GridEvaluatorStatus.h>
#include <vtkm/filter/flow/worklet/IntegratorStatus.h>
namespace vtkm
{
namespace worklet
{
namespace particleadvection
namespace flow
{
template <typename EvaluatorType>
@ -73,8 +76,8 @@ public:
}
}; //EulerIntegrator
} //namespace particleadvection
} //namespace worklet
} //namespace vtkm
}
}
} //vtkm::worklet::flow
#endif // vtk_m_worklet_particleadvection_EulerIntegrator_h
#endif // vtk_m_filter_flow_worklet_EulerIntegrator_h

@ -8,8 +8,8 @@
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#ifndef vtkm_woklet_particleadvection_field_h
#define vtkm_woklet_particleadvection_field_h
#ifndef vtkm_filter_flow_worklet_Field_h
#define vtkm_filter_flow_worklet_Field_h
#include <vtkm/Types.h>
@ -22,7 +22,7 @@ namespace vtkm
{
namespace worklet
{
namespace particleadvection
namespace flow
{
template <typename FieldArrayType>
@ -209,7 +209,8 @@ private:
Association Assoc;
};
} // namespace particleadvection
} // namespace worklet
} // namespace
#endif //vtkm_woklet_particleadvection_field_h
}
}
} //vtkm::worklet::flow
#endif //vtkm_filter_flow_worklet_Field_h

@ -8,25 +8,16 @@
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#ifndef vtk_m_worklet_particleadvection_GridEvaluatorStatus_h
#define vtk_m_worklet_particleadvection_GridEvaluatorStatus_h
#ifndef vtk_m_filter_flow_worklet_GridEvaluatorStatus_h
#define vtk_m_filter_flow_worklet_GridEvaluatorStatus_h
#include <vtkm/Bitset.h>
#include <vtkm/Types.h>
#include <vtkm/VectorAnalysis.h>
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/CellLocatorRectilinearGrid.h>
#include <vtkm/cont/CellLocatorTwoLevel.h>
#include <vtkm/cont/CellLocatorUniformGrid.h>
#include <vtkm/cont/CellSetStructured.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/DeviceAdapter.h>
namespace vtkm
{
namespace worklet
{
namespace particleadvection
namespace flow
{
class GridEvaluatorStatus : public vtkm::Bitset<vtkm::UInt8>
@ -62,8 +53,9 @@ private:
static constexpr vtkm::Id TEMPORAL_BOUNDS_BIT = 2;
static constexpr vtkm::Id IN_GHOST_CELL_BIT = 3;
};
}
}
}
#endif // vtk_m_worklet_particleadvection_GridEvaluatorStatus_h
#endif // vtk_m_filter_flow_worklet_GridEvaluatorStatus_h

@ -8,14 +8,10 @@
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#ifndef vtk_m_worklet_particleadvection_GridEvaluators_h
#define vtk_m_worklet_particleadvection_GridEvaluators_h
#ifndef vtk_m_filter_flow_worklet_GridEvaluators_h
#define vtk_m_filter_flow_worklet_GridEvaluators_h
#include <vtkm/Bitset.h>
#include <vtkm/CellClassification.h>
#include <vtkm/Types.h>
#include <vtkm/VectorAnalysis.h>
#include <vtkm/cont/ArrayCopy.h>
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/CellLocatorGeneral.h>
#include <vtkm/cont/CellLocatorRectilinearGrid.h>
@ -24,15 +20,15 @@
#include <vtkm/cont/CellSetStructured.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/worklet/particleadvection/CellInterpolationHelper.h>
#include <vtkm/worklet/particleadvection/Field.h>
#include <vtkm/worklet/particleadvection/GridEvaluatorStatus.h>
#include <vtkm/filter/flow/worklet/CellInterpolationHelper.h>
#include <vtkm/filter/flow/worklet/Field.h>
#include <vtkm/filter/flow/worklet/GridEvaluatorStatus.h>
namespace vtkm
{
namespace worklet
{
namespace particleadvection
namespace flow
{
template <typename FieldType>
@ -184,9 +180,9 @@ public:
{
this->InitializeLocator(dataSet.GetCoordinateSystem(), dataSet.GetCellSet());
if (dataSet.HasCellField("vtkmGhostCells"))
if (dataSet.HasGhostCellField())
{
auto arr = dataSet.GetCellField("vtkmGhostCells").GetData();
auto arr = dataSet.GetGhostCellField().GetData();
vtkm::cont::ArrayCopyShallowIfPossible(arr, this->GhostCellArray);
}
}
@ -232,8 +228,8 @@ private:
vtkm::cont::CellLocatorGeneral Locator;
};
} //namespace particleadvection
} //namespace worklet
} //namespace vtkm
}
}
} //vtkm::worklet::flow
#endif // vtk_m_worklet_particleadvection_GridEvaluators_h
#endif // vtk_m_filter_flow_worklet_GridEvaluators_h

@ -10,8 +10,8 @@
//
//=============================================================================
#ifndef vtk_m_worklet_particleadvection_Integrator_Status_h
#define vtk_m_worklet_particleadvection_Integrator_Status_h
#ifndef vtk_m_filter_flow_worklet_IntegratorStatus_h
#define vtk_m_filter_flow_worklet_IntegratorStatus_h
#include <iomanip>
#include <limits>
@ -21,13 +21,13 @@
#include <vtkm/Types.h>
#include <vtkm/VectorAnalysis.h>
#include <vtkm/worklet/particleadvection/GridEvaluatorStatus.h>
#include <vtkm/filter/flow/worklet/GridEvaluatorStatus.h>
namespace vtkm
{
namespace worklet
{
namespace particleadvection
namespace flow
{
class IntegratorStatus : public vtkm::Bitset<vtkm::UInt8>
@ -82,9 +82,9 @@ inline VTKM_CONT std::ostream& operator<<(std::ostream& s, const IntegratorStatu
<< " tm= " << status.CheckTemporalBounds() << " gc= " << status.CheckInGhostCell() << "]";
return s;
}
}
}
}
}
}
} //vtkm::worklet::flow
#endif // vtk_m_worklet_particleadvection_IntegratorStatus_h
#endif // vtk_m_filter_flow_worklet_IntegratorStatus_h

@ -8,19 +8,18 @@
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#ifndef vtk_m_worklet_ParticleAdvection_h
#define vtk_m_worklet_ParticleAdvection_h
#ifndef vtk_m_filter_flow_worklet_ParticleAdvection_h
#define vtk_m_filter_flow_worklet_ParticleAdvection_h
#include <vtkm/Types.h>
#include <vtkm/cont/ArrayCopy.h>
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/Invoker.h>
#include <vtkm/worklet/particleadvection/ParticleAdvectionWorklets.h>
#include <vtkm/filter/flow/worklet/ParticleAdvectionWorklets.h>
namespace vtkm
{
namespace worklet
{
namespace flow
{
namespace detail
{
@ -70,14 +69,25 @@ class ParticleAdvection
public:
ParticleAdvection() {}
template <typename IntegratorType, typename ParticleType, typename ParticleStorage>
void Run2(const IntegratorType& it,
vtkm::cont::ArrayHandle<ParticleType, ParticleStorage>& particles,
vtkm::Id MaxSteps,
ParticleAdvectionResult<ParticleType>& result)
{
vtkm::worklet::flow::ParticleAdvectionWorklet<IntegratorType, ParticleType> worklet;
worklet.Run(it, particles, MaxSteps);
result = ParticleAdvectionResult<ParticleType>(particles);
}
template <typename IntegratorType, typename ParticleType, typename ParticleStorage>
ParticleAdvectionResult<ParticleType> Run(
const IntegratorType& it,
vtkm::cont::ArrayHandle<ParticleType, ParticleStorage>& particles,
vtkm::Id MaxSteps)
{
vtkm::worklet::particleadvection::ParticleAdvectionWorklet<IntegratorType, ParticleType>
worklet;
vtkm::worklet::flow::ParticleAdvectionWorklet<IntegratorType, ParticleType> worklet;
worklet.Run(it, particles, MaxSteps);
return ParticleAdvectionResult<ParticleType>(particles);
@ -89,8 +99,7 @@ public:
const vtkm::cont::ArrayHandle<vtkm::Vec3f, PointStorage>& points,
vtkm::Id MaxSteps)
{
vtkm::worklet::particleadvection::ParticleAdvectionWorklet<IntegratorType, ParticleType>
worklet;
vtkm::worklet::flow::ParticleAdvectionWorklet<IntegratorType, ParticleType> worklet;
vtkm::cont::ArrayHandle<ParticleType> particles;
vtkm::cont::ArrayHandle<vtkm::Id> step, ids;
@ -148,7 +157,7 @@ public:
vtkm::cont::ArrayHandle<ParticleType, ParticleStorage>& particles,
vtkm::Id MaxSteps)
{
vtkm::worklet::particleadvection::StreamlineWorklet<IntegratorType, ParticleType> worklet;
vtkm::worklet::flow::StreamlineWorklet<IntegratorType, ParticleType> worklet;
vtkm::cont::ArrayHandle<vtkm::Vec3f> positions;
vtkm::cont::CellSetExplicit<> polyLines;
@ -158,7 +167,9 @@ public:
return StreamlineResult<ParticleType>(particles, positions, polyLines);
}
};
}
}
#endif // vtk_m_worklet_ParticleAdvection_h
}
}
} // vtkm::worklet::flow
#endif // vtk_m_filter_flow_worklet_ParticleAdvection_h

@ -8,23 +8,17 @@
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#ifndef vtk_m_worklet_particleadvection_ParticleAdvectionWorklets_h
#define vtk_m_worklet_particleadvection_ParticleAdvectionWorklets_h
#ifndef vtk_m_filter_flow_worklet_ParticleAdvectionWorklets_h
#define vtk_m_filter_flow_worklet_ParticleAdvectionWorklets_h
#include <vtkm/Types.h>
#include <vtkm/cont/Algorithm.h>
#include <vtkm/cont/ArrayCopy.h>
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/ArrayHandleCast.h>
#include <vtkm/cont/ArrayHandleCounting.h>
#include <vtkm/cont/CellSetExplicit.h>
#include <vtkm/cont/ConvertNumComponentsToOffsets.h>
#include <vtkm/cont/ExecutionObjectBase.h>
#include <vtkm/Particle.h>
#include <vtkm/filter/flow/worklet/Particles.h>
#include <vtkm/worklet/WorkletMapField.h>
#include <vtkm/worklet/particleadvection/Particles.h>
#include <vtkm/worklet/particleadvection/Stepper.h>
#ifdef VTKM_CUDA
#include <vtkm/cont/cuda/internal/ScopedCudaStackSize.h>
@ -34,7 +28,7 @@ namespace vtkm
{
namespace worklet
{
namespace particleadvection
namespace flow
{
class ParticleAdvectWorklet : public vtkm::worklet::WorkletMapField
@ -106,10 +100,10 @@ public:
vtkm::Id& MaxSteps)
{
using ParticleAdvectWorkletType = vtkm::worklet::particleadvection::ParticleAdvectWorklet;
using ParticleAdvectWorkletType = vtkm::worklet::flow::ParticleAdvectWorklet;
using ParticleWorkletDispatchType =
typename vtkm::worklet::DispatcherMapField<ParticleAdvectWorkletType>;
using ParticleArrayType = vtkm::worklet::particleadvection::Particles<ParticleType>;
using ParticleArrayType = vtkm::worklet::flow::Particles<ParticleType>;
vtkm::Id numSeeds = static_cast<vtkm::Id>(particles.GetNumberOfValues());
//Create and invoke the particle advection.
@ -186,10 +180,9 @@ public:
vtkm::cont::CellSetExplicit<>& polyLines)
{
using ParticleWorkletDispatchType = typename vtkm::worklet::DispatcherMapField<
vtkm::worklet::particleadvection::ParticleAdvectWorklet>;
using StreamlineArrayType =
vtkm::worklet::particleadvection::StateRecordingParticles<ParticleType>;
using ParticleWorkletDispatchType =
typename vtkm::worklet::DispatcherMapField<vtkm::worklet::flow::ParticleAdvectWorklet>;
using StreamlineArrayType = vtkm::worklet::flow::StateRecordingParticles<ParticleType>;
vtkm::cont::ArrayHandle<vtkm::Id> initialStepsTaken;
@ -238,8 +231,9 @@ public:
polyLines.Fill(positions.GetNumberOfValues(), cellTypes, connectivity, offsets);
}
};
}
}
} // namespace vtkm::worklet::particleadvection
#endif // vtk_m_worklet_particleadvection_ParticleAdvectionWorklets_h
}
}
} // namespace vtkm::worklet::flow
#endif // vtk_m_filter_flow_worklet_ParticleAdvectionWorklets_h

@ -8,23 +8,19 @@
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#ifndef vtk_m_worklet_particleadvection_Particles_h
#define vtk_m_worklet_particleadvection_Particles_h
#ifndef vtk_m_filter_flow_worklet_Particles_h
#define vtk_m_filter_flow_worklet_Particles_h
#include <vtkm/Particle.h>
#include <vtkm/Types.h>
#include <vtkm/cont/Algorithm.h>
#include <vtkm/cont/ArrayCopy.h>
#include <vtkm/cont/ArrayHandleConstant.h>
#include <vtkm/cont/ExecutionObjectBase.h>
#include <vtkm/worklet/particleadvection/IntegratorStatus.h>
#include <vtkm/filter/flow/worklet/IntegratorStatus.h>
namespace vtkm
{
namespace worklet
{
namespace particleadvection
namespace flow
{
template <typename ParticleType>
class ParticleExecutionObject
{
@ -66,7 +62,7 @@ public:
VTKM_EXEC
void StatusUpdate(const vtkm::Id& idx,
const vtkm::worklet::particleadvection::IntegratorStatus& status,
const vtkm::worklet::flow::IntegratorStatus& status,
vtkm::Id maxSteps)
{
ParticleType p(this->GetParticle(idx));
@ -116,10 +112,11 @@ template <typename ParticleType>
class Particles : public vtkm::cont::ExecutionObjectBase
{
public:
VTKM_CONT vtkm::worklet::particleadvection::ParticleExecutionObject<ParticleType>
PrepareForExecution(vtkm::cont::DeviceAdapterId device, vtkm::cont::Token& token) const
VTKM_CONT vtkm::worklet::flow::ParticleExecutionObject<ParticleType> PrepareForExecution(
vtkm::cont::DeviceAdapterId device,
vtkm::cont::Token& token) const
{
return vtkm::worklet::particleadvection::ParticleExecutionObject<ParticleType>(
return vtkm::worklet::flow::ParticleExecutionObject<ParticleType>(
this->ParticleArray, this->MaxSteps, device, token);
}
@ -221,10 +218,10 @@ public:
};
VTKM_CONT vtkm::worklet::particleadvection::StateRecordingParticleExecutionObject<ParticleType>
VTKM_CONT vtkm::worklet::flow::StateRecordingParticleExecutionObject<ParticleType>
PrepareForExecution(vtkm::cont::DeviceAdapterId device, vtkm::cont::Token& token) const
{
return vtkm::worklet::particleadvection::StateRecordingParticleExecutionObject<ParticleType>(
return vtkm::worklet::flow::StateRecordingParticleExecutionObject<ParticleType>(
this->ParticleArray,
this->HistoryArray,
this->ValidPointArray,
@ -276,9 +273,9 @@ protected:
};
} //namespace particleadvection
} //namespace worklet
} //namespace vtkm
}
}
} //vtkm::worklet::flow
#endif // vtk_m_worklet_particleadvection_Particles_h
#endif // vtk_m_filter_flow_worklet_Particles_h
//============================================================================

@ -10,14 +10,17 @@
//
//=============================================================================
#ifndef vtk_m_worklet_particleadvection_RK4Integrator_h
#define vtk_m_worklet_particleadvection_RK4Integrator_h
#ifndef vtk_m_filter_flow_worklet_RK4Integrator_h
#define vtk_m_filter_flow_worklet_RK4Integrator_h
#include <vtkm/filter/flow/worklet/GridEvaluatorStatus.h>
#include <vtkm/filter/flow/worklet/IntegratorStatus.h>
namespace vtkm
{
namespace worklet
{
namespace particleadvection
namespace flow
{
template <typename ExecEvaluatorType>
@ -112,8 +115,8 @@ public:
}
};
} //namespace particleadvection
} //namespace worklet
} //namespace vtkm
}
}
} //vtkm::worklet::flow
#endif // vtk_m_worklet_particleadvection_RK4Integrator_h
#endif // vtk_m_filter_flow_worklet_RK4Integrator_h

@ -10,27 +10,20 @@
//
//=============================================================================
#ifndef vtk_m_worklet_particleadvection_Stepper_h
#define vtk_m_worklet_particleadvection_Stepper_h
#ifndef vtk_m_filter_flow_worklet_Stepper_h
#define vtk_m_filter_flow_worklet_Stepper_h
#include <vtkm/filter/flow/worklet/GridEvaluators.h>
#include <vtkm/filter/flow/worklet/IntegratorStatus.h>
#include <vtkm/filter/flow/worklet/Particles.h>
#include <limits>
#include <vtkm/Bitset.h>
#include <vtkm/TypeTraits.h>
#include <vtkm/Types.h>
#include <vtkm/VectorAnalysis.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
#include <vtkm/worklet/particleadvection/IntegratorStatus.h>
#include <vtkm/worklet/particleadvection/Particles.h>
namespace vtkm
{
namespace worklet
{
namespace particleadvection
namespace flow
{
template <typename ExecIntegratorType, typename ExecEvaluatorType>
@ -194,8 +187,8 @@ public:
}
};
} //namespace particleadvection
} //namespace worklet
} //namespace vtkm
}
}
} //vtkm::worklet::flow
#endif // vtk_m_worklet_particleadvection_Stepper_h
#endif // vtk_m_filter_flow_worklet_Stepper_h

@ -7,15 +7,11 @@
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#ifndef vtk_m_worklet_streamsurface_h
#define vtk_m_worklet_streamsurface_h
#ifndef vtk_m_filter_flow_worklet_streamsurface_h
#define vtk_m_filter_flow_worklet_streamsurface_h
#include <typeinfo>
#include <vtkm/VectorAnalysis.h>
#include <vtkm/cont/Algorithm.h>
#include <vtkm/cont/ArrayHandleView.h>
#include <vtkm/cont/CellSetExplicit.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/worklet/DispatcherMapTopology.h>
#include <vtkm/worklet/WorkletMapField.h>
#include <vtkm/worklet/WorkletMapTopology.h>
@ -24,6 +20,8 @@ namespace vtkm
{
namespace worklet
{
namespace flow
{
class StreamSurface
{
@ -276,7 +274,9 @@ public:
private:
};
}
}
#endif // vtk_m_worklet_streamsurface_h
}
}
} //vtkm::worklet::flow
#endif // vtk_m_filter_flow_worklet_streamsurface_h

@ -8,26 +8,25 @@
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#ifndef vtk_m_worklet_particleadvection_TemporalGridEvaluators_h
#define vtk_m_worklet_particleadvection_TemporalGridEvaluators_h
#ifndef vtk_m_filter_flow_worklet_TemporalGridEvaluators_h
#define vtk_m_filter_flow_worklet_TemporalGridEvaluators_h
#include <vtkm/worklet/particleadvection/GridEvaluatorStatus.h>
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
#include <vtkm/filter/flow/worklet/GridEvaluatorStatus.h>
#include <vtkm/filter/flow/worklet/GridEvaluators.h>
namespace vtkm
{
namespace worklet
{
namespace particleadvection
namespace flow
{
template <typename FieldType>
class ExecutionTemporalGridEvaluator
{
private:
using GridEvaluator = vtkm::worklet::particleadvection::GridEvaluator<FieldType>;
using ExecutionGridEvaluator =
vtkm::worklet::particleadvection::ExecutionGridEvaluator<FieldType>;
using GridEvaluator = vtkm::worklet::flow::GridEvaluator<FieldType>;
using ExecutionGridEvaluator = vtkm::worklet::flow::ExecutionGridEvaluator<FieldType>;
public:
VTKM_CONT
@ -117,7 +116,7 @@ template <typename FieldType>
class TemporalGridEvaluator : public vtkm::cont::ExecutionObjectBase
{
private:
using GridEvaluator = vtkm::worklet::particleadvection::GridEvaluator<FieldType>;
using GridEvaluator = vtkm::worklet::flow::GridEvaluator<FieldType>;
public:
VTKM_CONT TemporalGridEvaluator() = default;
@ -178,8 +177,8 @@ private:
vtkm::FloatDefault TimeTwo;
};
} // namespace particleadvection
} // namespace worklet
} // namespace vtkm
}
}
} //vtkm::worklet::flow
#endif
#endif // vtk_m_filter_flow_worklet_TemporalGridEvaluators_h

@ -23,7 +23,6 @@ void appendPts(vtkm::cont::DataSetBuilderExplicitIterative& dsb,
ids.push_back(pid);
}
void TestTubeFilters()
{
using VecType = vtkm::Vec3f;
@ -43,6 +42,18 @@ void TestTubeFilters()
appendPts(dsb, VecType(2, 1, 0), ids);
dsb.AddCell(vtkm::CELL_SHAPE_POLY_LINE, ids);
//add some degenerate polylines.
//polyline with 1 point.
ids.clear();
appendPts(dsb, VecType(0, 0, 0), ids);
dsb.AddCell(vtkm::CELL_SHAPE_POLY_LINE, ids);
//polyline with coincident points.
ids.clear();
appendPts(dsb, VecType(0, 0, 0), ids);
appendPts(dsb, VecType(0, 0, 0), ids);
dsb.AddCell(vtkm::CELL_SHAPE_POLY_LINE, ids);
vtkm::cont::DataSet ds = dsb.Create();
std::vector<vtkm::FloatDefault> ptVar, cellVar;
@ -53,7 +64,6 @@ void TestTubeFilters()
cellVar.push_back(100);
cellVar.push_back(101);
//Polyline 2.
ptVar.push_back(10);
ptVar.push_back(11);
@ -61,6 +71,16 @@ void TestTubeFilters()
cellVar.push_back(110);
cellVar.push_back(111);
//Add some degenerate polylines.
//Polyline 3: (only 1 point)
ptVar.push_back(-1);
cellVar.push_back(-1);
//Polyline 4: (2 coincident points)
ptVar.push_back(-1);
ptVar.push_back(-1);
cellVar.push_back(-1);
cellVar.push_back(-1);
ds.AddPointField("pointVar", ptVar);
ds.AddCellField("cellVar", cellVar);
@ -93,7 +113,6 @@ void TestTubeFilters()
VTKM_TEST_ASSERT(portal.Get(i) == ptVals[static_cast<std::size_t>(i)],
"Wrong value for point field");
//Validate the cell field
vtkm::cont::ArrayHandle<vtkm::FloatDefault> cellArr;
output.GetField("cellVar").GetData().AsArrayHandle(cellArr);

@ -18,6 +18,7 @@
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/UnknownCellSet.h>
#include <vtkm/worklet/DispatcherMapTopology.h>
#include <vtkm/worklet/ScatterCounting.h>
#include <vtkm/worklet/WorkletMapField.h>
#include <vtkm/worklet/WorkletMapTopology.h>
@ -47,7 +48,8 @@ public:
FieldOut ptsPerPolyline,
FieldOut ptsPerTube,
FieldOut numTubeConnIds,
FieldOut linesPerPolyline);
FieldOut linesPerPolyline,
FieldOut validCell);
using ExecutionSignature = void(CellShape shapeType,
PointCount numPoints,
PointIndices ptIndices,
@ -56,7 +58,8 @@ public:
_4 ptsPerPolyline,
_5 ptsPerTube,
_6 numTubeConnIds,
_7 linesPerPolyline);
_7 linesPerPolyline,
_8 validCell);
using InputDomain = _1;
template <typename CellShapeTag, typename PointIndexType, typename InPointsType>
@ -68,11 +71,14 @@ public:
vtkm::Id& ptsPerPolyline,
vtkm::Id& ptsPerTube,
vtkm::Id& numTubeConnIds,
vtkm::Id& linesPerPolyline) const
vtkm::Id& linesPerPolyline,
vtkm::Id& validCell) const
{
// We only support polylines that contain 2 or more points.
vtkm::IdComponent numNonCoincidentPoints = 1;
vtkm::Vec3f p = inPts.Get(ptIndices[0]);
validCell = 0;
for (int i = 1; i < numPoints; ++i)
{
vtkm::Vec3f pNext = inPts.Get(ptIndices[i]);
@ -80,6 +86,7 @@ public:
{
numNonCoincidentPoints++;
p = pNext;
validCell = 1;
}
}
@ -101,6 +108,7 @@ public:
}
else
{
validCell = 0;
ptsPerPolyline = 0;
nonIncidentPtsPerPolyline = 0;
ptsPerTube = 0;
@ -138,7 +146,12 @@ public:
_3 polylineOffset,
_4 outNormals);
using InputDomain = _1;
using ScatterType = vtkm::worklet::ScatterCounting;
VTKM_CONT
static ScatterType MakeScatter(const vtkm::cont::ArrayHandle<vtkm::Id>& validCell)
{
return ScatterType(validCell);
}
template <typename InPointsType, typename PointIndexType>
VTKM_EXEC vtkm::IdComponent FindValidSegment(const InPointsType& inPts,
@ -151,7 +164,7 @@ public:
while (end < numPoints)
{
auto pe = inPts.Get(ptIndices[end]);
if (vtkm::Magnitude(pe - ps) > 0)
if (vtkm::Magnitude(pe - ps) > vtkm::Epsilon<vtkm::FloatDefault>())
return end - 1;
end++;
}
@ -309,6 +322,12 @@ public:
_7 outPts,
_8 outPointSrcIdx);
using InputDomain = _1;
using ScatterType = vtkm::worklet::ScatterCounting;
VTKM_CONT
static ScatterType MakeScatter(const vtkm::cont::ArrayHandle<vtkm::Id>& validCell)
{
return ScatterType(validCell);
}
template <typename CellShapeTag,
typename PointIndexType,
@ -603,7 +622,8 @@ public:
}
//Count number of polyline pts, tube pts and tube cells
vtkm::cont::ArrayHandle<vtkm::Id> ptsPerPolyline, ptsPerTube, numTubeConnIds, segPerPolyline;
vtkm::cont::ArrayHandle<vtkm::Id> ptsPerPolyline, ptsPerTube, numTubeConnIds, segPerPolyline,
validCell;
vtkm::cont::ArrayHandle<vtkm::IdComponent> nonIncidentPtsPerPolyline;
CountSegments countSegs(this->Capping, this->NumSides);
vtkm::worklet::DispatcherMapTopology<CountSegments> countInvoker(countSegs);
@ -613,7 +633,8 @@ public:
ptsPerPolyline,
ptsPerTube,
numTubeConnIds,
segPerPolyline);
segPerPolyline,
validCell);
vtkm::Id totalPolylinePts = vtkm::cont::Algorithm::Reduce(ptsPerPolyline, vtkm::Id(0));
if (totalPolylinePts == 0)
@ -636,14 +657,16 @@ public:
//Generate normals at each point on all polylines
NormalsType normals;
normals.Allocate(totalPolylinePts);
vtkm::worklet::DispatcherMapTopology<GenerateNormals> genNormalsDisp;
vtkm::worklet::DispatcherMapTopology<GenerateNormals> genNormalsDisp(
GenerateNormals::MakeScatter(validCell));
genNormalsDisp.Invoke(cellset, coords, polylinePtOffset, normals);
//Generate the tube points
newPoints.Allocate(totalTubePts);
this->OutputPointSourceIndex.Allocate(totalTubePts);
GeneratePoints genPts(this->Capping, this->NumSides, this->Radius);
vtkm::worklet::DispatcherMapTopology<GeneratePoints> genPtsDisp(genPts);
vtkm::worklet::DispatcherMapTopology<GeneratePoints> genPtsDisp(
genPts, GeneratePoints::MakeScatter(validCell));
genPtsDisp.Invoke(cellset,
coords,
normals,

Some files were not shown because too many files have changed in this diff Show More