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

Merge topic 'PICS'

Browse Source 
8e500a17 Moving PortalType templetization to Step method of Integrators
33150365 Validate the data format for the UniformGridEvaluate.
69c1f946 Merge branch 'PICS' of gitlab.kitware.com:dpugmire/vtk-m into PICS
b87e8881 Merge branch 'PICS' of gitlab.kitware.com:dpugmire/vtk-m into PICS
5a652f41 Fix for MSVC compile isses
50dbd634 Remove unneeded template from grid eval glasses.
6f1e9b40 Fix compile warning.
32a0be8b Fix compile warnings/errors from CDash.
...

Acked-by: Kitware Robot <kwrobot@kitware.com>
Merge-request: !828
This commit is contained in:
Abhishek 2017-07-13 13:41:24 +00:00 committed by Kitware Robot
commit 58a08287bb
16 changed files with 2203 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
examples/CMakeLists.txt
View File

@ -32,6 +32,7 @@ add_subdirectory(isosurface)
add_subdirectory(multi_backend)
add_subdirectory(streamline)
add_subdirectory(tetrahedra)
add_subdirectory(particle_advection)
if(VTKm_ENABLE_RENDERING)
add_subdirectory(rendering)
endif()

47
examples/particle_advection/CMakeLists.txt Normal file
View File

@ -0,0 +1,47 @@
##=============================================================================
##
## Copyright (c) Kitware, Inc.
## All rights reserved.
## See LICENSE.txt for details.
##
## This software is distributed WITHOUT ANY WARRANTY; without even
## the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
## PURPOSE. See the above copyright notice for more information.
##
## Copyright 2015 Sandia Corporation.
## Copyright 2015 UT-Battelle, LLC.
## Copyright 2015 Los Alamos National Security.
##
## Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
## the U.S. Government retains certain rights in this software.
## Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
## Laboratory (LANL), the U.S. Government retains certain rights in
## this software.
##
##=============================================================================
#Find the VTK-m package
find_package(VTKm REQUIRED QUIET
OPTIONAL_COMPONENTS Serial CUDA TBB
)
add_executable(Particle_Advection_SERIAL ParticleAdvection.cxx)
target_include_directories(Particle_Advection_SERIAL PRIVATE ${VTKm_INCLUDE_DIRS})
target_link_libraries(Particle_Advection_SERIAL ${VTKm_LIBRARIES})
target_compile_options(Particle_Advection_SERIAL PRIVATE ${VTKm_COMPILE_OPTIONS})
if(VTKm_TBB_FOUND)
add_executable(Particle_Advection_TBB ParticleAdvectionTBB.cxx)
target_include_directories(Particle_Advection_TBB PRIVATE ${VTKm_INCLUDE_DIRS})
target_link_libraries(Particle_Advection_TBB ${VTKm_LIBRARIES})
target_compile_options(Particle_Advection_TBB PRIVATE ${VTKm_COMPILE_OPTIONS})
endif()
if(VTKm_CUDA_FOUND)
# Cuda compiles do not respect target_include_directories
cuda_include_directories(${VTKm_INCLUDE_DIRS})
cuda_add_executable(Particle_Advection_CUDA ParticleAdvection.cu)
target_include_directories(Particle_Advection_SERIAL PRIVATE ${VTKm_INCLUDE_DIRS})
target_link_libraries(Particle_Advection_CUDA PRIVATE ${VTKm_LIBRARIES})
target_compile_options(Particle_Advection_CUDA PRIVATE ${VTKm_COMPILE_OPTIONS})
endif()

23
examples/particle_advection/ParticleAdvection.cu Normal file
View File

@ -0,0 +1,23 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//
// Copyright 2014 Sandia Corporation.
// Copyright 2014 UT-Battelle, LLC.
// Copyright 2014 Los Alamos National Security.
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#define VTKM_DEVICE_ADAPTER VTKM_DEVICE_ADAPTER_CUDA
#include "ParticleAdvection.cxx"

328
examples/particle_advection/ParticleAdvection.cxx Normal file
View File

@ -0,0 +1,328 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//
// Copyright 2014 Sandia Corporation.
// Copyright 2014 UT-Battelle, LLC.
// Copyright 2014 Los Alamos National Security.
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#ifndef VTKM_DEVICE_ADAPTER
#define VTKM_DEVICE_ADAPTER VTKM_DEVICE_ADAPTER_SERIAL
#endif
#include <vtkm/cont/DataSet.h>
#include <vtkm/worklet/ParticleAdvection.h>
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
#include <vtkm/worklet/particleadvection/Integrators.h>
#include <vtkm/worklet/particleadvection/ParticleAdvectionWorklets.h>
#include <vtkm/worklet/particleadvection/Particles.h>
#include <vtkm/cont/Timer.h>
#include <vtkm/io/reader/BOVDataSetReader.h>
#include <chrono>
#include <vector>
const vtkm::Id SPARSE = 0;
const vtkm::Id DENSE = 1;
const vtkm::Id MEDIUM = 2;
template <typename T>
static vtkm::Range subRange(vtkm::Range& range, T a, T b)
{
vtkm::Float32 arg1, arg2, len;
arg1 = static_cast<vtkm::Float32>(a);
arg2 = static_cast<vtkm::Float32>(b);
len = static_cast<vtkm::Float32>(range.Length());
return vtkm::Range(range.Min + arg1 * len, range.Min + arg2 * len);
}
template <typename T>
void ignore(T&&)
{
}
void RunTest(const std::string& fname,
vtkm::Id numSeeds,
vtkm::Id numSteps,
vtkm::Float32 stepSize,
vtkm::Id numThreads,
vtkm::Id advectType,
vtkm::Id stepsPerRound,
vtkm::Id particlesPerRound,
vtkm::Id seeding)
{
typedef VTKM_DEFAULT_DEVICE_ADAPTER_TAG DeviceAdapter;
typedef vtkm::Float32 FieldType;
typedef vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>> FieldHandle;
typedef
typename FieldHandle::template ExecutionTypes<DeviceAdapter>::PortalConst FieldPortalConstType;
vtkm::io::reader::BOVDataSetReader rdr(fname);
vtkm::cont::DataSet ds = rdr.ReadDataSet();
typedef vtkm::worklet::particleadvection::UniformGridEvaluate<FieldPortalConstType, FieldType>
RGEvalType;
typedef vtkm::worklet::particleadvection::RK4Integrator<RGEvalType, FieldType> RK4RGType;
RGEvalType eval(ds);
RK4RGType rk4(eval, stepSize);
std::vector<vtkm::Vec<FieldType, 3>> seeds;
srand(314);
vtkm::Bounds bounds = ds.GetCoordinateSystem().GetBounds();
if (seeding == SPARSE)
bounds = ds.GetCoordinateSystem().GetBounds();
else if (seeding == DENSE)
{
if (fname.find("astro") != std::string::npos)
{
bounds.X = subRange(bounds.X, .1, .15);
bounds.Y = subRange(bounds.Y, .1, .15);
bounds.Z = subRange(bounds.Z, .1, .15);
}
else if (fname.find("fusion") != std::string::npos)
{
bounds.X = subRange(bounds.X, .8, .85);
bounds.Y = subRange(bounds.Y, .55, .60);
bounds.Z = subRange(bounds.Z, .55, .60);
}
else if (fname.find("fishtank") != std::string::npos)
{
bounds.X = subRange(bounds.X, .1, .15);
bounds.Y = subRange(bounds.Y, .1, .15);
bounds.Z = subRange(bounds.Z, .55, .60);
}
}
else if (seeding == MEDIUM)
{
if (fname.find("astro") != std::string::npos)
{
bounds.X = subRange(bounds.X, .4, .6);
bounds.Y = subRange(bounds.Y, .4, .6);
bounds.Z = subRange(bounds.Z, .4, .6);
}
else if (fname.find("fusion") != std::string::npos)
{
bounds.X = subRange(bounds.X, .01, .99);
bounds.Y = subRange(bounds.Y, .01, .99);
bounds.Z = subRange(bounds.Z, .45, .55);
}
else if (fname.find("fishtank") != std::string::npos)
{
bounds.X = subRange(bounds.X, .4, .6);
bounds.Y = subRange(bounds.Y, .4, .6);
bounds.Z = subRange(bounds.Z, .4, .6);
}
}
for (int i = 0; i < numSeeds; i++)
{
vtkm::Vec<FieldType, 3> p;
vtkm::Float32 rx = (vtkm::Float32)rand() / (vtkm::Float32)RAND_MAX;
vtkm::Float32 ry = (vtkm::Float32)rand() / (vtkm::Float32)RAND_MAX;
vtkm::Float32 rz = (vtkm::Float32)rand() / (vtkm::Float32)RAND_MAX;
p[0] = static_cast<FieldType>(bounds.X.Min + rx * bounds.X.Length());
p[1] = static_cast<FieldType>(bounds.Y.Min + ry * bounds.Y.Length());
p[2] = static_cast<FieldType>(bounds.Z.Min + rz * bounds.Z.Length());
seeds.push_back(p);
}
#ifdef __BUILDING_TBB_VERSION__
int nT = tbb::task_scheduler_init::default_num_threads();
if (numThreads != -1)
nT = (int)numThreads;
//make sure the task_scheduler_init object is in scope when running sth w/ TBB
tbb::task_scheduler_init init(nT);
#else
ignore(numThreads);
#endif
//time only the actual run
auto t0 = std::chrono::high_resolution_clock::now();
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>> seedArray;
seedArray = vtkm::cont::make_ArrayHandle(seeds);
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>> fieldArray;
ds.GetField(0).GetData().CopyTo(fieldArray);
if (advectType == 0)
{
vtkm::worklet::ParticleAdvection particleAdvection;
particleAdvection.Run(rk4, seedArray, fieldArray, numSteps, -1, DeviceAdapter());
}
else
{
vtkm::worklet::Streamline streamline;
streamline.Run(
rk4, seedArray, fieldArray, numSteps, stepsPerRound, particlesPerRound, DeviceAdapter());
}
auto t1 = std::chrono::high_resolution_clock::now() - t0;
auto runtime = std::chrono::duration_cast<std::chrono::milliseconds>(t1).count();
std::cerr << "Runtime = " << runtime << " ms " << std::endl;
}
bool ParseArgs(int argc,
char** argv,
vtkm::Id& numSeeds,
vtkm::Id& numSteps,
vtkm::Float32& stepSize,
vtkm::Id& advectType,
vtkm::Id& stepsPerRound,
vtkm::Id& particlesPerRound,
vtkm::Id& numThreads,
std::string& dataFile,
std::string& pgmType,
bool& dumpOutput,
vtkm::Id& seeding)
{
numSeeds = 100;
numSteps = 100;
stepSize = 0.1f;
advectType = 0;
stepsPerRound = -1;
particlesPerRound = -1;
numThreads = -1;
dataFile = "";
pgmType = "UNKNOWN";
dumpOutput = false;
seeding = SPARSE;
if (argc < 2)
{
std::cerr << "Usage " << argv[0] << std::endl;
std::cerr << " -seeds #seeds" << std::endl;
std::cerr << " -steps maxSteps" << std::endl;
std::cerr << " -h stepSize" << std::endl;
std::cerr << " -particle : particle push" << std::endl;
std::cerr << " -streamline steps_per_round (-1 = 0 rounds): particle history" << std::endl;
std::cerr << " -t #numThreads" << std::endl;
std::cerr << " -file dataFile" << std::endl;
std::cerr << " -dump : dump output points" << std::endl;
return false;
}
std::string pgm = argv[0];
if (pgm.find("SERIAL") != std::string::npos)
pgmType = "SER";
else if (pgm.find("TBB") != std::string::npos)
pgmType = "TBB";
else if (pgm.find("CUDA") != std::string::npos)
pgmType = "CUD";
for (int i = 1; i < argc; i++)
{
std::string arg = argv[i];
if (arg == "-seeds")
numSeeds = static_cast<vtkm::Id>(atoi(argv[++i]));
else if (arg == "-steps")
numSteps = static_cast<vtkm::Id>(atoi(argv[++i]));
else if (arg == "-h")
stepSize = static_cast<vtkm::Float32>(atof(argv[++i]));
else if (arg == "-particle")
advectType = 0;
else if (arg == "-streamline")
{
advectType = 1;
}
else if (arg == "-streamlineS")
{
advectType = 1;
stepsPerRound = static_cast<vtkm::Id>(atoi(argv[++i]));
}
else if (arg == "-streamlineP")
{
advectType = 1;
particlesPerRound = static_cast<vtkm::Id>(atoi(argv[++i]));
}
else if (arg == "-streamlineSP")
{
advectType = 1;
stepsPerRound = static_cast<vtkm::Id>(atoi(argv[++i]));
particlesPerRound = static_cast<vtkm::Id>(atoi(argv[++i]));
}
else if (arg == "-file")
dataFile = argv[++i];
else if (arg == "-t")
numThreads = static_cast<vtkm::Id>(atoi(argv[++i]));
else if (arg == "-dump")
dumpOutput = true;
else if (arg == "-sparse")
seeding = SPARSE;
else if (arg == "-dense")
seeding = DENSE;
else if (arg == "-medium")
seeding = MEDIUM;
else
std::cerr << "Unexpected argument: " << arg << std::endl;
}
if (dataFile.size() == 0)
{
std::cerr << "Error: no data file specified" << std::endl;
return false;
}
//Congratulations user, we have a valid run:
std::cerr << pgmType << ": " << numSeeds << " " << numSteps << " " << stepSize << " ";
if (advectType == 0)
std::cerr << "PP ";
else
std::cerr << "SL ";
std::cerr << numThreads << " ";
std::cerr << dataFile << std::endl;
return true;
}
int main(int argc, char** argv)
{
vtkm::Id numSeeds = 100, numSteps = 100, advectType = 0, numThreads = -1, stepsPerRound = -1,
particlesPerRound = -1;
vtkm::Float32 stepSize = 0.1f;
std::string dataFile, pgmType;
vtkm::Id seeding = SPARSE;
bool dumpOutput = false;
if (!ParseArgs(argc,
argv,
numSeeds,
numSteps,
stepSize,
advectType,
stepsPerRound,
particlesPerRound,
numThreads,
dataFile,
pgmType,
dumpOutput,
seeding))
{
return -1;
}
RunTest(dataFile,
numSeeds,
numSteps,
stepSize,
numThreads,
advectType,
stepsPerRound,
particlesPerRound,
seeding);
return 0;
}

26
examples/particle_advection/ParticleAdvectionTBB.cxx Normal file
View File

@ -0,0 +1,26 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//
// Copyright 2014 Sandia Corporation.
// Copyright 2014 UT-Battelle, LLC.
// Copyright 2014 Los Alamos National Security.
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#define VTKM_DEVICE_ADAPTER VTKM_DEVICE_ADAPTER_TBB
#define __BUILDING_TBB_VERSION__
#include <tbb/task_scheduler_init.h>
#include "ParticleAdvection.cxx"

291
vtkm/io/reader/BOVDataSetReader.h Normal file
View File

@ -0,0 +1,291 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//
// Copyright 2015 Sandia Corporation.
// Copyright 2015 UT-Battelle, LLC.
// Copyright 2015 Los Alamos National Security.
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#ifndef vtk_m_io_reader_BOVDataSetReader_h
#define vtk_m_io_reader_BOVDataSetReader_h
#include <fstream>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/DataSetBuilderUniform.h>
#include <vtkm/cont/DataSetFieldAdd.h>
#include <vtkm/io/ErrorIO.h>
namespace vtkm
{
namespace io
{
namespace reader
{
class BOVDataSetReader
{
public:
BOVDataSetReader(const char* fileName)
: FileName(fileName)
, Loaded(false)
, DataSet()
{
}
BOVDataSetReader(const std::string& fileName)
: FileName(fileName)
, Loaded(false)
, DataSet()
{
}
const vtkm::cont::DataSet& ReadDataSet()
{
try
{
LoadFile();
}
catch (std::ifstream::failure& e)
{
std::string message("IO Error: ");
throw vtkm::io::ErrorIO(message + e.what());
}
return this->DataSet;
}
private:
typedef enum { ByteData, ShortData, IntegerData, FloatData, DoubleData } DataFormat;
void LoadFile()
{
if (this->Loaded)
return;
std::ifstream stream(this->FileName);
if (stream.fail())
throw vtkm::io::ErrorIO("Failed to open file: " + this->FileName);
DataFormat dataFormat;
std::string bovFile, line, token, options, variableName;
vtkm::Id numComponents = 1;
vtkm::Id3 dim;
vtkm::Vec<vtkm::FloatDefault, 3> origin(0, 0, 0);
vtkm::Vec<vtkm::FloatDefault, 3> spacing(1, 1, 1);
bool spacingSet = false;
while (stream.good())
{
std::getline(stream, line);
if (line.size() == 0 || line[0] == '#')
continue;
//std::cout<<"::"<<line<<"::"<<std::endl;
std::size_t pos = line.find(":");
if (pos == std::string::npos)
throw vtkm::io::ErrorIO("Unsupported option: " + line);
token = line.substr(0, pos);
options = line.substr(pos + 1, line.size() - 1);
//std::cout<<token<<"::"<<options<<std::endl;
std::stringstream strStream(options);
//Format supports both space and "_" seperated tokens...
if (token.find("DATA") != std::string::npos && token.find("FILE") != std::string::npos)
{
strStream >> bovFile >> std::ws;
}
else if (token.find("DATA") != std::string::npos && token.find("SIZE") != std::string::npos)
{
strStream >> dim[0] >> dim[1] >> dim[2] >> std::ws;
}
else if (token.find("BRICK") != std::string::npos &&
token.find("ORIGIN") != std::string::npos)
{
strStream >> origin[0] >> origin[1] >> origin[2] >> std::ws;
}
//DRP
else if (token.find("BRICK") != std::string::npos && token.find("SIZE") != std::string::npos)
{
strStream >> spacing[0] >> spacing[1] >> spacing[2] >> std::ws;
spacingSet = true;
}
else if (token.find("DATA") != std::string::npos && token.find("FORMAT") != std::string::npos)
{
std::string opt;
strStream >> opt >> std::ws;
if (opt.find("FLOAT") != std::string::npos || opt.find("REAL") != std::string::npos)
dataFormat = FloatData;
else if (opt.find("DOUBLE") != std::string::npos)
dataFormat = DoubleData;
else
throw vtkm::io::ErrorIO("Unsupported data type: " + token);
}
else if (token.find("DATA") != std::string::npos &&
token.find("COMPONENTS") != std::string::npos)
{
strStream >> numComponents >> std::ws;
if (numComponents != 1 && numComponents != 3)
throw vtkm::io::ErrorIO("Unsupported number of components");
}
else if (token.find("VARIABLE") != std::string::npos &&
token.find("PALETTE") == std::string::npos)
{
strStream >> variableName >> std::ws;
if (variableName[0] == '"')
variableName = variableName.substr(1, variableName.size() - 2);
}
/*
else
std::cerr<<"Unsupported BOV option: "<<token<<std::endl;
*/
}
if (spacingSet)
{
spacing[0] = (spacing[0]) / static_cast<vtkm::FloatDefault>(dim[0] - 1);
spacing[1] = (spacing[1]) / static_cast<vtkm::FloatDefault>(dim[1] - 1);
spacing[2] = (spacing[2]) / static_cast<vtkm::FloatDefault>(dim[2] - 1);
}
std::string fullPathDataFile;
std::size_t pos = FileName.rfind("/");
if (pos != std::string::npos)
{
std::string baseDir;
baseDir = this->FileName.substr(0, pos);
fullPathDataFile = baseDir + "/" + bovFile;
}
else
fullPathDataFile = bovFile;
/*
//Get whole path for data file.
std::string fullPathDataFile;
if (bovFile[0] == '/')
fullPathDataFile = bovFile;
else
{
//Get base dir.
std::string baseDir, baseFile;
std::cout<<FileName<<std::endl;
std::size_t pos = FileName.rfind("/");
if (pos != std::string::npos)
baseDir = this->FileName.substr(0, pos);
std::cout<<"BASE: "<<baseDir<<std::endl;
if (bovFile.substr(0,2) == "./")
{
baseFile = bovFile.substr(2, bovFile.size()-2);
}
if (baseDir.size() == 0)
fullPathDataFile = baseFile;
else
fullPathDataFile = baseDir + "/" + baseFile;
std::cout<<baseDir<<" : "<<baseFile<<std::endl;
std::cout<<fullPathDataFile<<std::endl;
}
*/
vtkm::cont::DataSetBuilderUniform dataSetBuilder;
vtkm::cont::DataSetFieldAdd dsf;
this->DataSet = dataSetBuilder.Create(dim, origin, spacing);
vtkm::Id numTuples = dim[0] * dim[1] * dim[2];
if (numComponents == 1)
{
if (dataFormat == FloatData)
{
vtkm::cont::ArrayHandle<vtkm::Float32> var;
ReadScalar(fullPathDataFile, numTuples, var);
dsf.AddPointField(this->DataSet, variableName, var);
}
else if (dataFormat == DoubleData)
{
vtkm::cont::ArrayHandle<vtkm::Float64> var;
ReadScalar(fullPathDataFile, numTuples, var);
dsf.AddPointField(this->DataSet, variableName, var);
}
}
else if (numComponents == 3)
{
if (dataFormat == FloatData)
{
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float32, 3>> var;
ReadVector(fullPathDataFile, numTuples, var);
dsf.AddPointField(this->DataSet, variableName, var);
}
else if (dataFormat == DoubleData)
{
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float64, 3>> var;
ReadVector(fullPathDataFile, numTuples, var);
dsf.AddPointField(this->DataSet, variableName, var);
}
}
this->Loaded = true;
}
template <typename T>
void ReadBuffer(const std::string& fName, const vtkm::Id& sz, std::vector<T>& buff)
{
FILE* fp = fopen(fName.c_str(), "rb");
size_t readSize = static_cast<size_t>(sz);
if (fp == NULL)
throw vtkm::io::ErrorIO("Unable to open data file: " + fName);
buff.resize(readSize);
size_t nread = fread(&buff[0], sizeof(T), readSize, fp);
if (nread != readSize)
throw vtkm::io::ErrorIO("Data file read failed: " + fName);
fclose(fp);
}
template <typename T>
void ReadScalar(const std::string& fName,
const vtkm::Id& nTuples,
vtkm::cont::ArrayHandle<T>& var)
{
std::vector<T> buff;
ReadBuffer(fName, nTuples, buff);
var.Allocate(nTuples);
for (vtkm::Id i = 0; i < nTuples; i++)
var.GetPortalControl().Set(i, buff[(size_t)i]);
}
template <typename T>
void ReadVector(const std::string& fName,
const vtkm::Id& nTuples,
vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>>& var)
{
std::vector<T> buff;
ReadBuffer(fName, nTuples * 3, buff);
var.Allocate(nTuples);
vtkm::Vec<T, 3> v;
for (vtkm::Id i = 0; i < nTuples; i++)
{
v[0] = buff[static_cast<size_t>(i * 3 + 0)];
v[1] = buff[static_cast<size_t>(i * 3 + 1)];
v[2] = buff[static_cast<size_t>(i * 3 + 2)];
var.GetPortalControl().Set(i, v);
}
}
std::string FileName;
bool Loaded;
vtkm::cont::DataSet DataSet;
};
}
}
} // vtkm::io::reader
#endif // vtk_m_io_reader_BOVReader_h

1
vtkm/io/reader/CMakeLists.txt
View File

@ -19,6 +19,7 @@
##============================================================================
set(headers
BOVDataSetReader.h
VTKDataSetReaderBase.h
VTKPolyDataReader.h
VTKStructuredGridReader.h

2
vtkm/worklet/CMakeLists.txt
View File

@ -44,6 +44,7 @@ set(headers
MarchingCubesDataTables.h
Mask.h
MaskPoints.h
ParticleAdvection.h
PointAverage.h
PointElevation.h
RemoveUnusedPoints.h
@ -73,6 +74,7 @@ add_subdirectory(spatialstructure)
add_subdirectory(tetrahedralize)
add_subdirectory(triangulate)
add_subdirectory(wavelets)
add_subdirectory(particleadvection)
vtkm_declare_headers(${headers})

87
vtkm/worklet/ParticleAdvection.h Normal file
View File

@ -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.
//
// Copyright 2014 Sandia Corporation.
// Copyright 2014 UT-Battelle, LLC.
// Copyright 2014 Los Alamos National Security.
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#ifndef vtk_m_worklet_ParticleAdvection_h
#define vtk_m_worklet_ParticleAdvection_h
#include <vtkm/Types.h>
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/worklet/particleadvection/ParticleAdvectionWorklets.h>
namespace vtkm
{
namespace worklet
{
class ParticleAdvection
{
public:
ParticleAdvection() {}
template <typename IntegratorType,
typename FieldType,
typename PointStorage,
typename FieldStorage,
typename DeviceAdapter>
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, PointStorage> Run(
const IntegratorType& it,
const vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, PointStorage>& pts,
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, FieldStorage> fieldArray,
const vtkm::Id& nSteps,
const vtkm::Id& particlesPerRound,
const DeviceAdapter&)
{
vtkm::worklet::particleadvection::ParticleAdvectionWorklet<IntegratorType,
FieldType,
DeviceAdapter>
worklet;
return worklet.Run(it, pts, fieldArray, nSteps, particlesPerRound);
}
};
class Streamline
{
public:
Streamline() {}
template <typename IntegratorType,
typename FieldType,
typename PointStorage,
typename FieldStorage,
typename DeviceAdapter>
void Run(const IntegratorType& it,
const vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, PointStorage>& pts,
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, FieldStorage> fieldArray,
const vtkm::Id& nSteps,
const vtkm::Id& stepsPerRound,
const vtkm::Id& particlesPerRound,
const DeviceAdapter&)
{
vtkm::worklet::particleadvection::StreamlineWorklet<IntegratorType, FieldType, DeviceAdapter>
worklet;
worklet.Run(it, pts, fieldArray, nSteps, stepsPerRound, particlesPerRound);
}
};
}
}
#endif // vtk_m_worklet_ParticleAdvection_h

29
vtkm/worklet/particleadvection/CMakeLists.txt Normal file
View File

@ -0,0 +1,29 @@
##============================================================================
## Copyright (c) Kitware, Inc.
## All rights reserved.
## See LICENSE.txt for details.
## This software is distributed WITHOUT ANY WARRANTY; without even
## the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
## PURPOSE. See the above copyright notice for more information.
##
## Copyright 2016 Sandia Corporation.
## Copyright 2016 UT-Battelle, LLC.
## Copyright 2016 Los Alamos National Security.
##
## Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
## the U.S. Government retains certain rights in this software.
##
## Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
## Laboratory (LANL), the U.S. Government retains certain rights in
## this software.
##============================================================================
set(headers
GridEvaluators.h
Integrators.h
Particles.h
ParticleAdvectionWorklets.h
)
#-----------------------------------------------------------------------------
vtkm_declare_headers(${headers})

407
vtkm/worklet/particleadvection/GridEvaluators.h Normal file
View File

@ -0,0 +1,407 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//
// Copyright 2014 Sandia Corporation.
// Copyright 2014 UT-Battelle, LLC.
// Copyright 2014 Los Alamos National Security.
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#ifndef vtk_m_worklet_particleadvection_GridEvaluators_h
#define vtk_m_worklet_particleadvection_GridEvaluators_h
#include <vtkm/Types.h>
#include <vtkm/VectorAnalysis.h>
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/CellSetStructured.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/DeviceAdapter.h>
#include <vtkm/cont/DynamicArrayHandle.h>
namespace vtkm
{
namespace worklet
{
namespace particleadvection
{
//Base class for all grid evaluators
class GridEvaluate
{
public:
enum Status
{
OK = 0,
OUTSIDE_SPATIAL,
OUTSIDE_TEMPORAL,
FAIL
};
};
// Constant vector
template <typename PortalType, typename FieldType>
class ConstantField
{
public:
VTKM_CONT
ConstantField(const vtkm::Bounds& bb, const vtkm::Vec<FieldType, 3>& v)
: bounds{ bb }
, vector{ v }
{
}
VTKM_EXEC_CONT
bool Evaluate(const vtkm::Vec<FieldType, 3>& pos,
const PortalType& vtkmNotUsed(vecData),
vtkm::Vec<FieldType, 3>& out) const
{
if (!bounds.Contains(pos))
return false;
out[0] = vector[0];
out[1] = vector[1];
out[2] = vector[2];
return true;
}
private:
vtkm::Bounds bounds;
vtkm::Vec<FieldType, 3> vector;
};
// Circular Orbit.
template <typename PortalType>
class AnalyticalOrbitEvaluate
{
public:
VTKM_CONT
AnalyticalOrbitEvaluate(const vtkm::Bounds& bb)
: bounds{ bb }
{
}
template <typename FieldType>
VTKM_EXEC_CONT bool Evaluate(const vtkm::Vec<FieldType, 3>& pos,
const PortalType& vtkmNotUsed(vecData),
vtkm::Vec<FieldType, 3>& out) const
{
if (!bounds.Contains(pos))
return false;
//statically return a value which is orthogonal to the input pos in the xy plane.
FieldType oneDivLen = 1.0f / Magnitude(pos);
out[0] = -1.0f * pos[1] * oneDivLen;
out[1] = pos[0] * oneDivLen;
out[2] = pos[2] * oneDivLen;
return true;
}
private:
vtkm::Bounds bounds;
};
//Uniform Grid Evaluator
template <typename PortalType, typename FieldType>
class UniformGridEvaluate
{
public:
VTKM_CONT
UniformGridEvaluate() {}
VTKM_CONT
UniformGridEvaluate(const vtkm::cont::DataSet& ds)
{
typedef vtkm::cont::ArrayHandleUniformPointCoordinates UniformType;
vtkm::cont::DynamicArrayHandleCoordinateSystem coordArray = ds.GetCoordinateSystem().GetData();
if (!coordArray.IsSameType(UniformType()))
throw vtkm::cont::ErrorInternal("Given dataset is was not uniform.");
bounds = ds.GetCoordinateSystem(0).GetBounds();
vtkm::cont::CellSetStructured<3> cells;
ds.GetCellSet(0).CopyTo(cells);
dims = cells.GetSchedulingRange(vtkm::TopologyElementTagPoint());
vtkm::Vec<FieldType, 3> castdims;
castdims[0] = static_cast<FieldType>(dims[0]);
castdims[1] = static_cast<FieldType>(dims[1]);
castdims[2] = static_cast<FieldType>(dims[2]);
spacing[0] = static_cast<FieldType>((bounds.X.Max - bounds.X.Min) / (castdims[0] - 1));
spacing[1] = static_cast<FieldType>((bounds.Y.Max - bounds.Y.Min) / (castdims[1] - 1));
spacing[2] = static_cast<FieldType>((bounds.Z.Max - bounds.Z.Min) / (castdims[2] - 1));
oldMin[0] =
static_cast<FieldType>(bounds.X.Min / ((bounds.X.Max - bounds.X.Min) / castdims[0]));
oldMin[1] =
static_cast<FieldType>(bounds.Y.Min / ((bounds.Y.Max - bounds.Y.Min) / castdims[1]));
oldMin[2] =
static_cast<FieldType>(bounds.Z.Min / ((bounds.Z.Max - bounds.Z.Min) / castdims[2]));
planeSize = dims[0] * dims[1];
rowSize = dims[0];
}
VTKM_EXEC_CONT
bool Evaluate(const vtkm::Vec<FieldType, 3>& pos,
const PortalType& vecData,
vtkm::Vec<FieldType, 3>& out) const
{
if (!bounds.Contains(pos))
return false;
// Set the eight corner indices with no wraparound
vtkm::Id3 idx000, idx001, idx010, idx011, idx100, idx101, idx110, idx111;
vtkm::Vec<FieldType, 3> normalizedPos;
normalizedPos[0] = pos[0] / spacing[0];
normalizedPos[1] = pos[1] / spacing[1];
normalizedPos[2] = pos[2] / spacing[2];
idx000[0] = static_cast<vtkm::Id>(floor(normalizedPos[0] - oldMin[0]));
idx000[1] = static_cast<vtkm::Id>(floor(normalizedPos[1] - oldMin[1]));
idx000[2] = static_cast<vtkm::Id>(floor(normalizedPos[2] - oldMin[2]));
idx001 = idx000;
idx001[0] = (idx001[0] + 1) <= dims[0] - 1 ? idx001[0] + 1 : dims[0] - 1;
idx010 = idx000;
idx010[1] = (idx010[1] + 1) <= dims[1] - 1 ? idx010[1] + 1 : dims[1] - 1;
idx011 = idx010;
idx011[0] = (idx011[0] + 1) <= dims[0] - 1 ? idx011[0] + 1 : dims[0] - 1;
idx100 = idx000;
idx100[2] = (idx100[2] + 1) <= dims[2] - 1 ? idx100[2] + 1 : dims[2] - 1;
idx101 = idx100;
idx101[0] = (idx101[0] + 1) <= dims[0] - 1 ? idx101[0] + 1 : dims[0] - 1;
idx110 = idx100;
idx110[1] = (idx110[1] + 1) <= dims[1] - 1 ? idx110[1] + 1 : dims[1] - 1;
idx111 = idx110;
idx111[0] = (idx111[0] + 1) <= dims[0] - 1 ? idx111[0] + 1 : dims[0] - 1;
// Get the vecdata at the eight corners
vtkm::Vec<FieldType, 3> v000, v001, v010, v011, v100, v101, v110, v111;
v000 = vecData.Get(idx000[2] * planeSize + idx000[1] * rowSize + idx000[0]);
v001 = vecData.Get(idx001[2] * planeSize + idx001[1] * rowSize + idx001[0]);
v010 = vecData.Get(idx010[2] * planeSize + idx010[1] * rowSize + idx010[0]);
v011 = vecData.Get(idx011[2] * planeSize + idx011[1] * rowSize + idx011[0]);
v100 = vecData.Get(idx100[2] * planeSize + idx100[1] * rowSize + idx100[0]);
v101 = vecData.Get(idx101[2] * planeSize + idx101[1] * rowSize + idx101[0]);
v110 = vecData.Get(idx110[2] * planeSize + idx110[1] * rowSize + idx110[0]);
v111 = vecData.Get(idx111[2] * planeSize + idx111[1] * rowSize + idx111[0]);
// Interpolation in X
vtkm::Vec<FieldType, 3> v00, v01, v10, v11;
FieldType a = pos[0] - static_cast<FieldType>(floor(pos[0]));
v00[0] = (1.0f - a) * v000[0] + a * v001[0];
v00[1] = (1.0f - a) * v000[1] + a * v001[1];
v00[2] = (1.0f - a) * v000[2] + a * v001[2];
v01[0] = (1.0f - a) * v010[0] + a * v011[0];
v01[1] = (1.0f - a) * v010[1] + a * v011[1];
v01[2] = (1.0f - a) * v010[2] + a * v011[2];
v10[0] = (1.0f - a) * v100[0] + a * v101[0];
v10[1] = (1.0f - a) * v100[1] + a * v101[1];
v10[2] = (1.0f - a) * v100[2] + a * v101[2];
v11[0] = (1.0f - a) * v110[0] + a * v111[0];
v11[1] = (1.0f - a) * v110[1] + a * v111[1];
v11[2] = (1.0f - a) * v110[2] + a * v111[2];
// Interpolation in Y
vtkm::Vec<FieldType, 3> v0, v1;
a = pos[1] - static_cast<FieldType>(floor(pos[1]));
v0[0] = (1.0f - a) * v00[0] + a * v01[0];
v0[1] = (1.0f - a) * v00[1] + a * v01[1];
v0[2] = (1.0f - a) * v00[2] + a * v01[2];
v1[0] = (1.0f - a) * v10[0] + a * v11[0];
v1[1] = (1.0f - a) * v10[1] + a * v11[1];
v1[2] = (1.0f - a) * v10[2] + a * v11[2];
a = pos[2] - static_cast<FieldType>(floor(pos[2]));
out[0] = (1.0f - a) * v0[0] + v1[0];
out[1] = (1.0f - a) * v0[1] + v1[1];
out[2] = (1.0f - a) * v0[2] + v1[2];
return true;
}
private:
vtkm::Bounds bounds;
vtkm::Id3 dims;
vtkm::Id planeSize;
vtkm::Id rowSize;
vtkm::Vec<FieldType, 3> spacing;
vtkm::Vec<FieldType, 3> oldMin;
};
template <typename PortalType, typename FieldType>
class RectilinearGridEvaluate
{
public:
VTKM_CONT
RectilinearGridEvaluate(const vtkm::cont::DataSet& dataset)
{
bounds = dataset.GetCoordinateSystem(0).GetBounds();
vtkm::cont::CellSetStructured<3> cells;
dataset.GetCellSet(0).CopyTo(cells);
dims = cells.GetSchedulingRange(vtkm::TopologyElementTagPoint());
planeSize = dims[0] * dims[1];
rowSize = dims[0];
vtkm::cont::DynamicArrayHandleCoordinateSystem coordArray =
dataset.GetCoordinateSystem().GetData();
if (coordArray.IsSameType(RectilinearType()))
{
RectilinearType gridPoints = coordArray.Cast<RectilinearType>();
xAxis = gridPoints.GetPortalConstControl().GetFirstPortal();
yAxis = gridPoints.GetPortalConstControl().GetSecondPortal();
zAxis = gridPoints.GetPortalConstControl().GetThirdPortal();
}
else
{
/*
* As the data is not in the rectilinear format.
* The code will not be able to continue unless
* the data is in the required format.
*/
throw vtkm::cont::ErrorInternal("Given dataset is was not rectilinear.");
}
}
VTKM_EXEC_CONT
bool Evaluate(const vtkm::Vec<FieldType, 3>& pos,
const PortalType& vecData,
vtkm::Vec<FieldType, 3>& out) const
{
if (!bounds.Contains(pos))
return false;
vtkm::Id3 idx000, idx001, idx010, idx011, idx100, idx101, idx110, idx111;
vtkm::Vec<vtkm::Id, 3> cellPos;
vtkm::Id index;
/*Get floor X location*/
for (index = 0; index < dims[0] - 1; index++)
{
if (xAxis.Get(index) <= pos[0] && pos[0] < xAxis.Get(index + 1))
{
cellPos[0] = index;
break;
}
}
/*Get floor Y location*/
for (index = 0; index < dims[1] - 1; index++)
{
if (yAxis.Get(index) <= pos[1] && pos[1] < yAxis.Get(index + 1))
{
cellPos[1] = index;
break;
}
}
/*Get floor Z location*/
for (index = 0; index < dims[2] - 1; index++)
{
if (zAxis.Get(index) <= pos[2] && pos[2] < zAxis.Get(index + 1))
{
cellPos[2] = index;
break;
}
}
idx000[0] = cellPos[0];
idx000[1] = cellPos[1];
idx000[2] = cellPos[2];
idx001 = idx000;
idx001[0] = (idx001[0] + 1) <= dims[0] - 1 ? idx001[0] + 1 : dims[0] - 1;
idx010 = idx000;
idx010[1] = (idx010[1] + 1) <= dims[1] - 1 ? idx010[1] + 1 : dims[1] - 1;
idx011 = idx010;
idx011[0] = (idx011[0] + 1) <= dims[0] - 1 ? idx011[0] + 1 : dims[0] - 1;
idx100 = idx000;
idx100[2] = (idx100[2] + 1) <= dims[2] - 1 ? idx100[2] + 1 : dims[2] - 1;
idx101 = idx100;
idx101[0] = (idx101[0] + 1) <= dims[0] - 1 ? idx101[0] + 1 : dims[0] - 1;
idx110 = idx100;
idx110[1] = (idx110[1] + 1) <= dims[1] - 1 ? idx110[1] + 1 : dims[1] - 1;
idx111 = idx110;
idx111[0] = (idx111[0] + 1) <= dims[0] - 1 ? idx111[0] + 1 : dims[0] - 1;
// Get the vecdata at the eight corners
vtkm::Vec<FieldType, 3> v000, v001, v010, v011, v100, v101, v110, v111;
v000 = vecData.Get(idx000[2] * planeSize + idx000[1] * rowSize + idx000[0]);
v001 = vecData.Get(idx001[2] * planeSize + idx001[1] * rowSize + idx001[0]);
v010 = vecData.Get(idx010[2] * planeSize + idx010[1] * rowSize + idx010[0]);
v011 = vecData.Get(idx011[2] * planeSize + idx011[1] * rowSize + idx011[0]);
v100 = vecData.Get(idx100[2] * planeSize + idx100[1] * rowSize + idx100[0]);
v101 = vecData.Get(idx101[2] * planeSize + idx101[1] * rowSize + idx101[0]);
v110 = vecData.Get(idx110[2] * planeSize + idx110[1] * rowSize + idx110[0]);
v111 = vecData.Get(idx111[2] * planeSize + idx111[1] * rowSize + idx111[0]);
// Interpolation in X
vtkm::Vec<FieldType, 3> v00, v01, v10, v11;
FieldType a = pos[0] - static_cast<FieldType>(floor(pos[0]));
v00[0] = (1.0f - a) * v000[0] + a * v001[0];
v00[1] = (1.0f - a) * v000[1] + a * v001[1];
v00[2] = (1.0f - a) * v000[2] + a * v001[2];
v01[0] = (1.0f - a) * v010[0] + a * v011[0];
v01[1] = (1.0f - a) * v010[1] + a * v011[1];
v01[2] = (1.0f - a) * v010[2] + a * v011[2];
v10[0] = (1.0f - a) * v100[0] + a * v101[0];
v10[1] = (1.0f - a) * v100[1] + a * v101[1];
v10[2] = (1.0f - a) * v100[2] + a * v101[2];
v11[0] = (1.0f - a) * v110[0] + a * v111[0];
v11[1] = (1.0f - a) * v110[1] + a * v111[1];
v11[2] = (1.0f - a) * v110[2] + a * v111[2];
// Interpolation in Y
vtkm::Vec<FieldType, 3> v0, v1;
a = pos[1] - static_cast<FieldType>(floor(pos[1]));
v0[0] = (1.0f - a) * v00[0] + a * v01[0];
v0[1] = (1.0f - a) * v00[1] + a * v01[1];
v0[2] = (1.0f - a) * v00[2] + a * v01[2];
v1[0] = (1.0f - a) * v10[0] + a * v11[0];
v1[1] = (1.0f - a) * v10[1] + a * v11[1];
v1[2] = (1.0f - a) * v10[2] + a * v11[2];
// Interpolation in Z
a = pos[2] - static_cast<FieldType>(floor(pos[2]));
out[0] = (1.0f - a) * v0[0] + v1[0];
out[1] = (1.0f - a) * v0[1] + v1[1];
out[2] = (1.0f - a) * v0[2] + v1[2];
return true;
}
private:
typedef vtkm::cont::ArrayHandle<FieldType> AxisHandle;
typedef vtkm::cont::ArrayHandleCartesianProduct<AxisHandle, AxisHandle, AxisHandle>
RectilinearType;
typename AxisHandle::PortalConstControl xAxis;
typename AxisHandle::PortalConstControl yAxis;
typename AxisHandle::PortalConstControl zAxis;
vtkm::Bounds bounds;
vtkm::Id3 dims;
vtkm::Id planeSize;
vtkm::Id rowSize;
}; //RectilinearGridEvaluate
} //namespace particleadvection
} //namespace worklet
} //namespace vtkm
#endif // vtk_m_worklet_particleadvection_GridEvaluators_h

104
vtkm/worklet/particleadvection/Integrators.h Normal file
View File

@ -0,0 +1,104 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//
// Copyright 2014 Sandia Corporation.
// Copyright 2014 UT-Battelle, LLC.
// Copyright 2014 Los Alamos National Security.
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#ifndef vtk_m_worklet_particleadvection_Integrators_h
#define vtk_m_worklet_particleadvection_Integrators_h
#include <vtkm/Types.h>
namespace vtkm
{
namespace worklet
{
namespace particleadvection
{
template <typename FieldEvaluateType, typename FieldType>
class RK4Integrator
{
public:
VTKM_EXEC_CONT
RK4Integrator()
: h(0)
, h_2(0)
{
}
VTKM_EXEC_CONT
RK4Integrator(const FieldEvaluateType& field, FieldType _h)
: f(field)
, h(_h)
, h_2(_h / 2.f)
{
}
template <typename PortalType>
VTKM_EXEC bool Step(const vtkm::Vec<FieldType, 3>& pos,
const PortalType& field,
vtkm::Vec<FieldType, 3>& out) const
{
vtkm::Vec<FieldType, 3> k1, k2, k3, k4;
if (f.Evaluate(pos, field, k1) && f.Evaluate(pos + h_2 * k1, field, k2) &&
f.Evaluate(pos + h_2 * k2, field, k3) && f.Evaluate(pos + h * k3, field, k4))
{
out = pos + h / 6.0f * (k1 + 2 * k2 + 2 * k3 + k4);
return true;
}
return false;
}
FieldEvaluateType f;
FieldType h, h_2;
};
template <typename FieldEvaluateType, typename FieldType>
class EulerIntegrator
{
public:
VTKM_EXEC_CONT
EulerIntegrator(const FieldEvaluateType& field, FieldType _h)
: f(field)
, h(_h)
{
}
template <typename PortalType>
VTKM_EXEC bool Step(const vtkm::Vec<FieldType, 3>& pos,
const PortalType& field,
vtkm::Vec<FieldType, 3>& out) const
{
vtkm::Vec<FieldType, 3> vCur;
if (f.Evaluate(pos, field, vCur))
{
out = pos + h * vCur;
return true;
}
return false;
}
FieldEvaluateType f;
FieldType h;
};
}
}
}
#endif // vtk_m_worklet_particleadvection_Integrators_h

326
vtkm/worklet/particleadvection/ParticleAdvectionWorklets.h Normal file
View File

@ -0,0 +1,326 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//
// Copyright 2014 Sandia Corporation.
// Copyright 2014 UT-Battelle, LLC.
// Copyright 2014 Los Alamos National Security.
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#ifndef vtk_m_worklet_particleadvection_ParticleAdvectionWorklets_h
#define vtk_m_worklet_particleadvection_ParticleAdvectionWorklets_h
#include <vtkm/Types.h>
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/ArrayHandleCounting.h>
#include <vtkm/cont/CellSetExplicit.h>
#include <vtkm/cont/CellSetStructured.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/Field.h>
#include <vtkm/exec/ExecutionObjectBase.h>
#include <vtkm/worklet/DispatcherMapField.h>
#include <vtkm/worklet/particleadvection/Particles.h>
namespace vtkm
{
namespace worklet
{
namespace particleadvection
{
template <typename IntegratorType, typename FieldType, typename DeviceAdapterTag>
class ParticleAdvectWorklet : public vtkm::worklet::WorkletMapField
{
public:
typedef vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>> FieldHandle;
typedef typename FieldHandle::template ExecutionTypes<DeviceAdapterTag>::PortalConst
FieldPortalConstType;
typedef void ControlSignature(FieldIn<IdType> idx, ExecObject ic);
typedef void ExecutionSignature(_1, _2);
typedef _1 InputDomain;
template <typename IntegralCurveType>
VTKM_EXEC void operator()(const vtkm::Id& idx, IntegralCurveType& ic) const
{
vtkm::Vec<FieldType, 3> p = ic.GetPos(idx);
vtkm::Vec<FieldType, 3> p2;
while (!ic.Done(idx))
{
if (integrator.Step(p, field, p2))
{
ic.TakeStep(idx, p2);
p = p2;
}
else
{
ic.SetStatusOutOfBounds(idx);
}
}
}
ParticleAdvectWorklet(const IntegratorType& it, const FieldPortalConstType& f)
: integrator(it)
, field(f)
{
}
IntegratorType integrator;
FieldPortalConstType field;
};
template <typename IntegratorType, typename FieldType, typename DeviceAdapterTag>
class ParticleAdvectionWorklet
{
public:
typedef vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>> FieldHandle;
typedef typename FieldHandle::template ExecutionTypes<DeviceAdapterTag>::PortalConst
FieldPortalConstType;
typedef vtkm::worklet::particleadvection::ParticleAdvectWorklet<IntegratorType,
FieldType,
DeviceAdapterTag>
ParticleAdvectWorkletType;
ParticleAdvectionWorklet() {}
template <typename PointStorage, typename FieldStorage>
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, PointStorage> Run(
const IntegratorType& it,
const vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, PointStorage>& pts,
const vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, FieldStorage> fieldArray,
const vtkm::Id& nSteps,
const vtkm::Id& particlesPerRound = -1)
{
integrator = it;
seedArray = pts;
maxSteps = nSteps;
ParticlesPerRound = particlesPerRound;
field = fieldArray.PrepareForInput(DeviceAdapterTag());
return run();
}
~ParticleAdvectionWorklet() {}
private:
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>> run(bool dumpOutput = false)
{
typedef typename vtkm::worklet::DispatcherMapField<ParticleAdvectWorkletType>
ParticleWorkletDispatchType;
typedef vtkm::worklet::particleadvection::Particles<FieldType, DeviceAdapterTag> ParticleType;
vtkm::Id totNumSeeds = static_cast<vtkm::Id>(seedArray.GetNumberOfValues());
vtkm::Id numSeeds = totNumSeeds;
if (ParticlesPerRound == -1 || ParticlesPerRound > totNumSeeds)
numSeeds = totNumSeeds;
else
numSeeds = ParticlesPerRound;
std::vector<vtkm::Id> steps(static_cast<size_t>(numSeeds), 0),
status(static_cast<size_t>(numSeeds), ParticleStatus::OK);
vtkm::cont::ArrayHandle<vtkm::Id> stepArray = vtkm::cont::make_ArrayHandle(&steps[0], numSeeds);
vtkm::cont::ArrayHandle<vtkm::Id> statusArray =
vtkm::cont::make_ArrayHandle(&status[0], numSeeds);
vtkm::cont::ArrayHandleIndex idxArray(numSeeds);
ParticleType particles(seedArray, stepArray, statusArray, maxSteps);
ParticleAdvectWorkletType particleWorklet(integrator, field);
ParticleWorkletDispatchType particleWorkletDispatch(particleWorklet);
particleWorkletDispatch.Invoke(idxArray, particles);
if (dumpOutput)
{
for (vtkm::Id i = 0; i < numSeeds; i++)
{
vtkm::Vec<FieldType, 3> p = particles.GetPos(i);
std::cout << p[0] << " " << p[1] << " " << p[2] << std::endl;
}
}
return seedArray;
}
IntegratorType integrator;
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>> seedArray;
vtkm::cont::DataSet ds;
vtkm::Id maxSteps;
vtkm::Id ParticlesPerRound;
FieldPortalConstType field;
};
template <typename IntegratorType, typename FieldType, typename DeviceAdapterTag>
class StreamlineWorklet
{
public:
typedef vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>> FieldHandle;
typedef typename FieldHandle::template ExecutionTypes<DeviceAdapterTag>::PortalConst
FieldPortalConstType;
typedef vtkm::worklet::particleadvection::ParticleAdvectWorklet<IntegratorType,
FieldType,
DeviceAdapterTag>
ParticleAdvectWorkletType;
StreamlineWorklet() {}
template <typename PointStorage, typename FieldStorage>
void Run(const IntegratorType& it,
const vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, PointStorage>& pts,
const vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>, FieldStorage> fieldArray,
const vtkm::Id& nSteps,
const vtkm::Id stepsPerRound = -1,
const vtkm::Id& particlesPerRound = -1)
{
integrator = it;
seedArray = pts;
maxSteps = nSteps;
StepsPerRound = stepsPerRound;
ParticlesPerRound = particlesPerRound;
field = fieldArray.PrepareForInput(DeviceAdapterTag());
run(true);
}
~StreamlineWorklet() {}
private:
void run(bool dumpOutput = false)
{
typedef typename vtkm::worklet::DispatcherMapField<ParticleAdvectWorkletType>
ParticleWorkletDispatchType;
typedef vtkm::worklet::particleadvection::StateRecordingParticles<FieldType, DeviceAdapterTag>
StreamlineType;
typedef vtkm::worklet::particleadvection::StateRecordingParticlesRound<FieldType,
DeviceAdapterTag>
StreamlineRoundType;
vtkm::Id totNumSeeds = static_cast<vtkm::Id>(seedArray.GetNumberOfValues());
bool NeedParticleRounds = false;
if (!(ParticlesPerRound == -1 || ParticlesPerRound > totNumSeeds))
NeedParticleRounds = true;
ParticleAdvectWorkletType particleWorklet(integrator, field);
ParticleWorkletDispatchType particleWorkletDispatch(particleWorklet);
//Brute force method, or rounds.
if (StepsPerRound == -1)
{
bool particlesDone = false;
vtkm::Id particleOffset = 0;
while (!particlesDone)
{
vtkm::Id num = totNumSeeds - particleOffset;
if (num <= 0)
break;
if (NeedParticleRounds && num > ParticlesPerRound)
num = ParticlesPerRound;
std::vector<vtkm::Id> steps((size_t)num, 0), status((size_t)num, ParticleStatus::OK);
vtkm::cont::ArrayHandle<vtkm::Id> stepArray = vtkm::cont::make_ArrayHandle(&steps[0], num);
vtkm::cont::ArrayHandle<vtkm::Id> statusArray =
vtkm::cont::make_ArrayHandle(&status[0], num);
vtkm::cont::ArrayHandleIndex idxArray(num);
StreamlineType streamlines(seedArray, stepArray, statusArray, maxSteps);
particleWorkletDispatch.Invoke(idxArray, streamlines);
if (dumpOutput)
{
for (vtkm::Id i = 0; i < num; i++)
{
vtkm::Id ns = streamlines.GetStep(i);
for (vtkm::Id j = 0; j < ns; j++)
{
vtkm::Vec<FieldType, 3> p = streamlines.GetHistory(i, j);
std::cout << p[0] << " " << p[1] << " " << p[2] << std::endl;
}
}
}
particleOffset += ParticlesPerRound;
if (!NeedParticleRounds)
particlesDone = true;
}
}
else
{
bool particlesDone = false;
vtkm::Id particleOffset = 0;
while (!particlesDone)
{
vtkm::Id num = totNumSeeds - particleOffset;
if (num <= 0)
break;
if (NeedParticleRounds && num > ParticlesPerRound)
num = ParticlesPerRound;
std::vector<vtkm::Id> steps((size_t)num, 0), status((size_t)num, ParticleStatus::OK);
vtkm::cont::ArrayHandle<vtkm::Id> stepArray = vtkm::cont::make_ArrayHandle(&steps[0], num);
vtkm::cont::ArrayHandle<vtkm::Id> statusArray =
vtkm::cont::make_ArrayHandle(&status[0], num);
vtkm::cont::ArrayHandleIndex idxArray(num);
vtkm::Id numSteps = 0, stepOffset = 0;
bool stepsDone = false;
while (!stepsDone)
{
numSteps += StepsPerRound;
if (numSteps >= maxSteps)
{
numSteps = maxSteps;
stepsDone = true;
}
StreamlineRoundType streamlines(
seedArray, stepArray, statusArray, numSteps, StepsPerRound, stepOffset, maxSteps);
particleWorkletDispatch.Invoke(idxArray, streamlines);
auto historyPortal = streamlines.historyArray.GetPortalConstControl();
if (dumpOutput)
{
for (vtkm::Id i = 0; i < num; i++)
{
vtkm::Id ns = streamlines.GetStep(i);
for (vtkm::Id j = stepOffset; j < ns; j++)
{
vtkm::Vec<FieldType, 3> p = historyPortal.Get(i * StepsPerRound + (j - stepOffset));
std::cout << p[0] << " " << p[1] << " " << p[2] << std::endl;
}
}
}
stepOffset += StepsPerRound;
}
particleOffset += ParticlesPerRound;
if (!NeedParticleRounds)
particlesDone = true;
}
}
}
IntegratorType integrator;
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>> seedArray;
vtkm::cont::DataSet ds;
vtkm::Id maxSteps;
vtkm::Id StepsPerRound, ParticlesPerRound;
FieldPortalConstType field;
};
}
}
}
#endif // vtk_m_worklet_particleadvection_ParticleAdvectionWorklets_h

374
vtkm/worklet/particleadvection/Particles.h Normal file
View File

@ -0,0 +1,374 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//
// Copyright 2014 Sandia Corporation.
// Copyright 2014 UT-Battelle, LLC.
// Copyright 2014 Los Alamos National Security.
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#ifndef vtk_m_worklet_particleadvection_Particles_h
#define vtk_m_worklet_particleadvection_Particles_h
#include <vtkm/Types.h>
#include <vtkm/cont/DeviceAdapter.h>
#include <vtkm/exec/ExecutionObjectBase.h>
namespace vtkm
{
namespace worklet
{
namespace particleadvection
{
enum ParticleStatus
{
OK = 0,
TERMINATE = 1,
OUT_OF_BOUNDS = 2,
};
template <typename T, typename DeviceAdapterTag>
class Particles : public vtkm::exec::ExecutionObjectBase
{
private:
typedef
typename vtkm::cont::ArrayHandle<vtkm::Id>::template ExecutionTypes<DeviceAdapterTag>::Portal
IdPortal;
typedef typename vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>>::template ExecutionTypes<
DeviceAdapterTag>::Portal PosPortal;
public:
VTKM_EXEC_CONT
Particles()
: pos()
, steps()
, status()
, maxSteps(0)
{
}
VTKM_EXEC_CONT
Particles(const Particles& ic)
: pos(ic.pos)
, steps(ic.steps)
, status(ic.status)
, maxSteps(ic.maxSteps)
{
}
VTKM_EXEC_CONT
Particles(const PosPortal& _pos,
const IdPortal& _steps,
const IdPortal& _status,
const vtkm::Id& _maxSteps)
: pos(_pos)
, steps(_steps)
, status(_status)
, maxSteps(_maxSteps)
{
}
VTKM_EXEC_CONT
Particles(vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>>& posArray,
vtkm::cont::ArrayHandle<vtkm::Id>& stepsArray,
vtkm::cont::ArrayHandle<vtkm::Id>& statusArray,
const vtkm::Id& _maxSteps)
: maxSteps(_maxSteps)
{
pos = posArray.PrepareForInPlace(DeviceAdapterTag());
steps = stepsArray.PrepareForInPlace(DeviceAdapterTag());
status = statusArray.PrepareForInPlace(DeviceAdapterTag());
}
VTKM_EXEC_CONT
void TakeStep(const vtkm::Id& idx, const vtkm::Vec<T, 3>& pt)
{
pos.Set(idx, pt);
vtkm::Id nSteps = steps.Get(idx);
nSteps = nSteps + 1;
steps.Set(idx, nSteps);
if (nSteps == maxSteps)
SetStatusTerminate(idx);
}
VTKM_EXEC_CONT
void SetStatusTerminate(const vtkm::Id& idx) { status.Set(idx, TERMINATE); }
VTKM_EXEC_CONT
void SetStatusOutOfBounds(const vtkm::Id& idx) { status.Set(idx, OUT_OF_BOUNDS); }
VTKM_EXEC_CONT
bool Done(const vtkm::Id& idx) { return status.Get(idx) != OK; }
VTKM_EXEC_CONT
vtkm::Vec<T, 3> GetPos(const vtkm::Id& idx) const { return pos.Get(idx); }
VTKM_EXEC_CONT
vtkm::Id GetStep(const vtkm::Id& idx) const { return steps.Get(idx); }
VTKM_EXEC_CONT
vtkm::Id GetStatus(const vtkm::Id& idx) const { return status.Get(idx); }
private:
PosPortal pos;
IdPortal steps, status;
vtkm::Id maxSteps;
};
template <typename T, typename DeviceAdapterTag>
class StateRecordingParticles : public vtkm::exec::ExecutionObjectBase
{
private:
typedef
typename vtkm::cont::ArrayHandle<vtkm::Id>::template ExecutionTypes<DeviceAdapterTag>::Portal
IdPortal;
typedef typename vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>>::template ExecutionTypes<
DeviceAdapterTag>::Portal PosPortal;
public:
VTKM_EXEC_CONT
StateRecordingParticles(const StateRecordingParticles& s)
: pos(s.pos)
, steps(s.steps)
, status(s.status)
, maxSteps(s.maxSteps)
, histSize(s.histSize)
, history(s.history)
{
}
VTKM_EXEC_CONT
StateRecordingParticles()
: pos()
, steps()
, status()
, maxSteps(0)
, histSize(-1)
, history()
{
}
VTKM_EXEC_CONT
StateRecordingParticles(const PosPortal& _pos,
const IdPortal& _steps,
const IdPortal& _status,
const vtkm::Id& _maxSteps)
: pos(_pos)
, steps(_steps)
, status(_status)
, maxSteps(_maxSteps)
, histSize()
, history()
{
}
VTKM_EXEC_CONT
StateRecordingParticles(vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>>& posArray,
vtkm::cont::ArrayHandle<vtkm::Id>& stepsArray,
vtkm::cont::ArrayHandle<vtkm::Id>& statusArray,
const vtkm::Id& _maxSteps)
: maxSteps(_maxSteps)
, histSize(_maxSteps)
{
pos = posArray.PrepareForInPlace(DeviceAdapterTag());
steps = stepsArray.PrepareForInPlace(DeviceAdapterTag());
status = statusArray.PrepareForInPlace(DeviceAdapterTag());
numPos = posArray.GetNumberOfValues();
history = historyArray.PrepareForOutput(numPos * histSize, DeviceAdapterTag());
}
VTKM_EXEC_CONT
StateRecordingParticles(vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>>& posArray,
vtkm::cont::ArrayHandle<vtkm::Id>& stepsArray,
vtkm::cont::ArrayHandle<vtkm::Id>& statusArray,
const vtkm::Id& _maxSteps,
vtkm::Id& _histSize)
: maxSteps(_maxSteps)
, histSize(_histSize)
{
pos = posArray.PrepareForInPlace(DeviceAdapterTag());
steps = stepsArray.PrepareForInPlace(DeviceAdapterTag());
status = statusArray.PrepareForInPlace(DeviceAdapterTag());
numPos = posArray.GetNumberOfValues();
history = historyArray.PrepareForOutput(numPos * histSize, DeviceAdapterTag());
}
VTKM_EXEC_CONT
void TakeStep(const vtkm::Id& idx, const vtkm::Vec<T, 3>& pt)
{
vtkm::Id nSteps = steps.Get(idx);
vtkm::Id loc = idx * histSize + nSteps;
history.Set(loc, pt);
nSteps = nSteps + 1;
steps.Set(idx, nSteps);
if (nSteps == maxSteps)
SetStatusTerminate(idx);
}
VTKM_EXEC_CONT
void SetStatusTerminate(const vtkm::Id& idx) { status.Set(idx, TERMINATE); }
VTKM_EXEC_CONT
void SetStatusOutOfBounds(const vtkm::Id& idx) { status.Set(idx, OUT_OF_BOUNDS); }
VTKM_EXEC_CONT
bool Done(const vtkm::Id& idx) { return status.Get(idx) != OK; }
VTKM_EXEC_CONT
vtkm::Vec<T, 3> GetPos(const vtkm::Id& idx) const { return pos.Get(idx); }
VTKM_EXEC_CONT
vtkm::Id GetStep(const vtkm::Id& idx) const { return steps.Get(idx); }
VTKM_EXEC_CONT
vtkm::Id GetStatus(const vtkm::Id& idx) const { return status.Get(idx); }
VTKM_EXEC_CONT
vtkm::Vec<T, 3> GetHistory(const vtkm::Id& idx, const vtkm::Id& step) const
{
return history.Get(idx * histSize + step);
}
private:
PosPortal pos;
IdPortal steps, status;
vtkm::Id maxSteps, numPos, histSize;
PosPortal history;
public:
vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>> historyArray;
};
template <typename T, typename DeviceAdapterTag>
class StateRecordingParticlesRound : public vtkm::exec::ExecutionObjectBase
{
private:
typedef
typename vtkm::cont::ArrayHandle<vtkm::Id>::template ExecutionTypes<DeviceAdapterTag>::Portal
IdPortal;
typedef typename vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>>::template ExecutionTypes<
DeviceAdapterTag>::Portal PosPortal;
public:
VTKM_EXEC_CONT
StateRecordingParticlesRound(const StateRecordingParticlesRound& s)
: pos(s.pos)
, steps(s.steps)
, status(s.status)
, maxSteps(s.maxSteps)
, numPos(s.numPos)
, histSize(s.histSize)
, offset(s.offset)
, totalMaxSteps(s.totalMaxSteps)
, history(s.history)
{
}
VTKM_EXEC_CONT
StateRecordingParticlesRound()
: pos()
, steps()
, maxSteps(0)
, histSize(-1)
, offset(0)
, totalMaxSteps(0)
{
}
VTKM_EXEC_CONT
StateRecordingParticlesRound(const PosPortal& _pos,
const IdPortal& _steps,
const IdPortal& _status,
const vtkm::Id& _maxSteps,
const vtkm::Id& _histSize,
const vtkm::Id& _offset,
const vtkm::Id& _totalMaxSteps)
: pos(_pos)
, steps(_steps)
, status(_status)
, maxSteps(_maxSteps)
, histSize(_histSize)
, offset(_offset)
, totalMaxSteps(_totalMaxSteps)
{
}
VTKM_EXEC_CONT
StateRecordingParticlesRound(vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>>& posArray,
vtkm::cont::ArrayHandle<vtkm::Id>& stepsArray,
vtkm::cont::ArrayHandle<vtkm::Id>& statusArray,
const vtkm::Id& _maxSteps,
const vtkm::Id& _histSize,
const vtkm::Id& _offset,
const vtkm::Id& _totalMaxSteps)
: maxSteps(_maxSteps)
, histSize(_histSize)
, offset(_offset)
, totalMaxSteps(_totalMaxSteps)
{
pos = posArray.PrepareForInPlace(DeviceAdapterTag());
steps = stepsArray.PrepareForInPlace(DeviceAdapterTag());
status = statusArray.PrepareForInPlace(DeviceAdapterTag());
numPos = posArray.GetNumberOfValues();
history = historyArray.PrepareForOutput(numPos * histSize, DeviceAdapterTag());
}
VTKM_EXEC_CONT
void TakeStep(const vtkm::Id& idx, const vtkm::Vec<T, 3>& pt)
{
vtkm::Id nSteps = steps.Get(idx);
vtkm::Id loc = idx * histSize + (nSteps - offset);
history.Set(loc, pt);
nSteps = nSteps + 1;
steps.Set(idx, nSteps);
if (nSteps == totalMaxSteps)
SetStatusTerminate(idx);
pos.Set(idx, pt);
}
VTKM_EXEC_CONT
void SetStatusTerminate(const vtkm::Id& idx) { status.Set(idx, TERMINATE); }
VTKM_EXEC_CONT
void SetStatusOutOfBounds(const vtkm::Id& idx) { status.Set(idx, OUT_OF_BOUNDS); }
VTKM_EXEC_CONT
bool Done(const vtkm::Id& idx)
{
vtkm::Id nSteps = steps.Get(idx);
return (nSteps - offset == histSize) || status.Get(idx) != OK;
}
VTKM_EXEC_CONT
vtkm::Vec<T, 3> GetPos(const vtkm::Id& idx) const { return pos.Get(idx); }
VTKM_EXEC_CONT
vtkm::Id GetStep(const vtkm::Id& idx) const { return steps.Get(idx); }
VTKM_EXEC_CONT
vtkm::Id GetStatus(const vtkm::Id& idx) const { return status.Get(idx); }
VTKM_EXEC_CONT
vtkm::Vec<T, 3> GetHistory(const vtkm::Id& idx, const vtkm::Id& step) const
{
return history.Get(idx * histSize + step);
}
private:
PosPortal pos;
IdPortal steps, status;
vtkm::Id maxSteps, numPos, histSize, offset, totalMaxSteps;
PosPortal history;
public:
vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>> historyArray;
};
}
}
}
#endif // vtk_m_worklet_particleadvection_Particles_h

1
vtkm/worklet/testing/CMakeLists.txt
View File

@ -36,6 +36,7 @@ set(unit_tests
UnitTestMarchingCubes.cxx
UnitTestMask.cxx
UnitTestMaskPoints.cxx
UnitTestParticleAdvection.cxx
UnitTestPointElevation.cxx
UnitTestPointGradient.cxx
UnitTestRemoveUnusedPoints.cxx

156
vtkm/worklet/testing/UnitTestParticleAdvection.cxx Normal file
View File

@ -0,0 +1,156 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//
// Copyright 2014 Sandia Corporation.
// Copyright 2014 UT-Battelle, LLC.
// Copyright 2014 Los Alamos National Security.
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/DataSetBuilderUniform.h>
#include <vtkm/cont/DeviceAdapter.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/exec/ExecutionWholeArray.h>
#include <vtkm/worklet/ParticleAdvection.h>
#include <vtkm/worklet/particleadvection/GridEvaluators.h>
#include <vtkm/worklet/particleadvection/Integrators.h>
#include <vtkm/worklet/particleadvection/Particles.h>
namespace
{
vtkm::Float32 vecData[125 * 3] = {
-0.00603248f, -0.0966396f, -0.000732792f, 0.000530014f, -0.0986189f, -0.000806706f,
0.00684929f, -0.100098f, -0.000876566f, 0.0129235f, -0.101102f, -0.000942341f,
0.0187515f, -0.101656f, -0.00100401f, 0.0706091f, -0.083023f, -0.00144278f,
0.0736404f, -0.0801616f, -0.00145784f, 0.0765194f, -0.0772063f, -0.00147036f,
0.0792559f, -0.0741751f, -0.00148051f, 0.0818589f, -0.071084f, -0.00148843f,
0.103585f, -0.0342287f, -0.001425f, 0.104472f, -0.0316147f, -0.00140433f,
0.105175f, -0.0291574f, -0.00138057f, 0.105682f, -0.0268808f, -0.00135357f,
0.105985f, -0.0248099f, -0.00132315f, -0.00244603f, -0.0989576f, -0.000821705f,
0.00389525f, -0.100695f, -0.000894513f, 0.00999301f, -0.10193f, -0.000963114f,
0.0158452f, -0.102688f, -0.00102747f, 0.0214509f, -0.102995f, -0.00108757f,
0.0708166f, -0.081799f, -0.00149941f, 0.0736939f, -0.0787879f, -0.00151236f,
0.0764359f, -0.0756944f, -0.00152297f, 0.0790546f, -0.0725352f, -0.00153146f,
0.0815609f, -0.0693255f, -0.001538f, -0.00914287f, -0.104658f, -0.001574f,
-0.00642891f, -0.10239f, -0.00159659f, -0.00402289f, -0.0994835f, -0.00160731f,
-0.00194792f, -0.0959752f, -0.00160528f, -0.00022818f, -0.0919077f, -0.00158957f,
-0.0134913f, -0.0274735f, -9.50056e-05f, -0.0188683f, -0.023273f, 0.000194107f,
-0.0254516f, -0.0197589f, 0.000529693f, -0.0312798f, -0.0179514f, 0.00083619f,
-0.0360426f, -0.0177537f, 0.00110164f, 0.0259929f, -0.0204479f, -0.000304646f,
0.033336f, -0.0157385f, -0.000505569f, 0.0403427f, -0.0104637f, -0.000693529f,
0.0469371f, -0.00477766f, -0.000865609f, 0.0530722f, 0.0011701f, -0.00102f,
-0.0121869f, -0.10317f, -0.0015868f, -0.0096549f, -0.100606f, -0.00160377f,
-0.00743038f, -0.0973796f, -0.00160783f, -0.00553901f, -0.0935261f, -0.00159792f,
-0.00400821f, -0.0890871f, -0.00157287f, -0.0267803f, -0.0165823f, 0.000454173f,
-0.0348303f, -0.011642f, 0.000881271f, -0.0424964f, -0.00870761f, 0.00129226f,
-0.049437f, -0.00781358f, 0.0016728f, -0.0552635f, -0.00888708f, 0.00200659f,
-0.0629746f, -0.0721524f, -0.00160475f, -0.0606813f, -0.0677576f, -0.00158427f,
-0.0582203f, -0.0625009f, -0.00154304f, -0.0555686f, -0.0563905f, -0.00147822f,
-0.0526988f, -0.0494369f, -0.00138643f, 0.0385695f, 0.115704f, 0.00674413f,
0.056434f, 0.128273f, 0.00869052f, 0.0775564f, 0.137275f, 0.0110399f,
0.102515f, 0.140823f, 0.0138637f, 0.131458f, 0.136024f, 0.0171804f,
0.0595175f, -0.0845927f, 0.00512454f, 0.0506615f, -0.0680369f, 0.00376604f,
0.0434904f, -0.0503557f, 0.00261592f, 0.0376711f, -0.0318716f, 0.00163301f,
0.0329454f, -0.0128019f, 0.000785352f, -0.0664062f, -0.0701094f, -0.00160644f,
-0.0641074f, -0.0658893f, -0.00158969f, -0.0616054f, -0.0608302f, -0.00155303f,
-0.0588734f, -0.0549447f, -0.00149385f, -0.0558797f, -0.0482482f, -0.00140906f,
0.0434062f, 0.102969f, 0.00581269f, 0.0619547f, 0.112838f, 0.00742057f,
0.0830229f, 0.118752f, 0.00927516f, 0.106603f, 0.119129f, 0.0113757f,
0.132073f, 0.111946f, 0.0136613f, -0.0135758f, -0.0934604f, -0.000533868f,
-0.00690763f, -0.0958773f, -0.000598878f, -0.000475275f, -0.0977838f, -0.000660985f,
0.00571866f, -0.0992032f, -0.0007201f, 0.0116724f, -0.10016f, -0.000776144f,
0.0651428f, -0.0850475f, -0.00120243f, 0.0682895f, -0.0823666f, -0.00121889f,
0.0712792f, -0.0795772f, -0.00123291f, 0.0741224f, -0.0766981f, -0.00124462f,
0.076829f, -0.0737465f, -0.00125416f, 0.10019f, -0.0375515f, -0.00121866f,
0.101296f, -0.0348723f, -0.00120216f, 0.102235f, -0.0323223f, -0.00118309f,
0.102994f, -0.0299234f, -0.00116131f, 0.103563f, -0.0276989f, -0.0011367f,
-0.00989236f, -0.0958821f, -0.000608883f, -0.00344154f, -0.0980645f, -0.000673641f,
0.00277318f, -0.0997337f, -0.000735354f, 0.00874908f, -0.100914f, -0.000793927f,
0.0144843f, -0.101629f, -0.000849279f, 0.0654428f, -0.0839355f, -0.00125739f,
0.0684225f, -0.0810989f, -0.00127208f, 0.0712599f, -0.0781657f, -0.00128444f,
0.0739678f, -0.0751541f, -0.00129465f, 0.076558f, -0.0720804f, -0.00130286f,
-0.0132841f, -0.103948f, -0.00131159f, -0.010344f, -0.102328f, -0.0013452f,
-0.00768637f, -0.100054f, -0.00136938f, -0.00533293f, -0.0971572f, -0.00138324f,
-0.00330643f, -0.0936735f, -0.00138586f, -0.0116984f, -0.0303752f, -0.000229102f,
-0.0149879f, -0.0265231f, -3.43823e-05f, -0.0212917f, -0.0219544f, 0.000270283f,
-0.0277756f, -0.0186879f, 0.000582781f, -0.0335115f, -0.0171098f, 0.00086919f,
0.0170095f, -0.025299f, -3.73557e-05f, 0.024552f, -0.0214351f, -0.000231975f,
0.0318714f, -0.0168568f, -0.000417463f, 0.0388586f, -0.0117131f, -0.000589883f,
0.0454388f, -0.00615626f, -0.000746594f, -0.0160785f, -0.102675f, -0.00132891f,
-0.0133174f, -0.100785f, -0.00135859f, -0.0108365f, -0.0982184f, -0.00137801f,
-0.00865931f, -0.0950053f, -0.00138614f, -0.00681126f, -0.0911806f, -0.00138185f,
-0.0208973f, -0.0216631f, 0.000111231f, -0.0289373f, -0.0151081f, 0.000512553f,
-0.0368736f, -0.0104306f, 0.000911793f, -0.0444294f, -0.00773838f, 0.00129762f,
-0.0512663f, -0.00706554f, 0.00165611f
};
}
void TestParticleAdvection()
{
typedef VTKM_DEFAULT_DEVICE_ADAPTER_TAG DeviceAdapter;
typedef vtkm::Float32 FieldType;
typedef vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>> FieldHandle;
typedef FieldHandle::template ExecutionTypes<DeviceAdapter>::PortalConst FieldPortalConstType;
std::cout << "Testing Integrators for ParticleAdvection Worklet" << std::endl;
FieldType stepSize = 0.01f;
vtkm::cont::DataSetBuilderUniform dataSetBuilder;
const vtkm::Id3 dims(5, 5, 5);
vtkm::Id nElements = dims[0] * dims[1] * dims[2] * 3;
std::vector<vtkm::Vec<vtkm::Float32, 3>> field;
for (vtkm::Id i = 0; i < nElements; i++)
{
FieldType x = vecData[i];
FieldType y = vecData[++i];
FieldType z = vecData[++i];
vtkm::Vec<FieldType, 3> vec(x, y, z);
field.push_back(vtkm::Normal(vec));
}
vtkm::cont::DataSet ds = dataSetBuilder.Create(dims);
typedef vtkm::worklet::particleadvection::UniformGridEvaluate<FieldPortalConstType, FieldType>
RGEvalType;
typedef vtkm::worklet::particleadvection::RK4Integrator<RGEvalType, FieldType> RK4RGType;
RGEvalType eval(ds);
RK4RGType rk4(eval, stepSize);
std::vector<vtkm::Vec<FieldType, 3>> pts;
pts.push_back(vtkm::Vec<FieldType, 3>(1, 1, 1));
pts.push_back(vtkm::Vec<FieldType, 3>(2, 2, 2));
pts.push_back(vtkm::Vec<FieldType, 3>(3, 3, 3));
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>> seeds;
seeds = vtkm::cont::make_ArrayHandle(pts);
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3>> fieldArray;
fieldArray = vtkm::cont::make_ArrayHandle(field);
vtkm::worklet::ParticleAdvection particleAdvection;
auto endSeeds = particleAdvection.Run(rk4, seeds, fieldArray, 100, -1, DeviceAdapter());
VTKM_TEST_ASSERT(seeds.GetNumberOfValues() == endSeeds.GetNumberOfValues(),
"Number of output particles does not match input.");
}
int UnitTestParticleAdvection(int, char* [])
{
return vtkm::cont::testing::Testing::Run(TestParticleAdvection);
}