Shuenhoy/vtk-m
1
0
Fork 0
mirror of https://gitlab.kitware.com/vtk/vtk-m synced 2024-09-19 18:45:43 +00:00
Code Issues 2 Actions Packages Projects Releases Wiki Activity

Merge topic 'streamline'

Browse Source 
6c4fb856 Add example data file for streamline.
cba0e218 Verify unit test results
3946b0c4 Add unit test, pass all args in Run()
09e43778 Merge branch 'master' of gitlab.kitware.com:Fasel/vtk-m into streamline
02f84a19 StreamLine filter outputs dataset.  Add example.
34cc9719 StreamLine particle tracing returns output dataset.
e34aaa02 Merge branch 'master' of gitlab.kitware.com:Fasel/vtk-m into streamline
82d9c102 Pass input and output datasets and not arrays. Support FORWARD, BACKWARD and BOTH.
...

Acked-by: Kitware Robot <kwrobot@kitware.com>
Merge-request: !267
This commit is contained in:
Patricia 2015-11-24 10:34:10 -05:00 committed by Kitware Robot
commit e033d5a391
10 changed files with 1074 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
examples/CMakeLists.txt
View File

@ -24,4 +24,5 @@ add_subdirectory(clipping)
add_subdirectory(hello_world)
add_subdirectory(isosurface)
add_subdirectory(multi_backend)
add_subdirectory(streamline)
add_subdirectory(tetrahedra)

39
examples/streamline/CMakeLists.txt Normal file
View File

@ -0,0 +1,39 @@
##=============================================================================
##
## 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.
##
##=============================================================================
if(OPENGL_FOUND AND GLUT_FOUND)
add_executable(StreamLineUniformGrid_SERIAL StreamLineUniformGrid.cxx)
target_include_directories(StreamLineUniformGrid_SERIAL PRIVATE ${GLUT_INCLUDE_DIR} ${OPENGL_INCLUDE_DIR})
target_link_libraries(StreamLineUniformGrid_SERIAL ${OPENGL_LIBRARIES} ${GLUT_LIBRARIES})
if(VTKm_Cuda_FOUND)
cuda_add_executable(StreamLineUniformGrid_CUDA StreamLineUniformGrid.cu)
target_link_libraries(StreamLineUniformGrid_CUDA ${OPENGL_LIBRARIES} ${GLUT_LIBRARIES})
endif()
if(VTKm_ENABLE_TBB)
add_executable(StreamLineUniformGrid_TBB StreamLineUniformGridTBB.cxx)
target_include_directories(StreamLineUniformGrid_TBB PRIVATE ${GLUT_INCLUDE_DIR} ${OPENGL_INCLUDE_DIR} ${TBB_INCLUDE_DIRS})
target_link_libraries(StreamLineUniformGrid_TBB ${OPENGL_LIBRARIES} ${GLUT_LIBRARIES} ${TBB_LIBRARIES})
endif()
endif()

24
examples/streamline/StreamLineUniformGrid.cu Normal file
View File

@ -0,0 +1,24 @@
//============================================================================
// 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
#define BOOST_SP_DISABLE_THREADS
#include "StreamLineUniformGrid.cxx"

304
examples/streamline/StreamLineUniformGrid.cxx Normal file
View File

@ -0,0 +1,304 @@
//============================================================================
// 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/worklet/StreamLineUniformGrid.h>
#include <vtkm/worklet/DispatcherMapField.h>
#include <vtkm/Math.h>
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/testing/Testing.h>
#include <fstream>
#include <vector>
#include <math.h>
//Suppress warnings about glut being deprecated on OSX
#if (defined(VTKM_GCC) || defined(VTKM_CLANG)) && !defined(VTKM_PGI)
# pragma GCC diagnostic push
# pragma GCC diagnostic ignored "-Wdeprecated-declarations"
#endif
#if defined (__APPLE__)
# include <GLUT/glut.h>
#else
# include <GL/glut.h>
#endif
#include "../isosurface/quaternion.h"
typedef VTKM_DEFAULT_DEVICE_ADAPTER_TAG DeviceAdapter;
// Output data set shared with opengl
vtkm::worklet::StreamLineFilterUniformGrid<vtkm::Float32, DeviceAdapter> *streamLineFilter;
vtkm::cont::DataSet outDataSet;
// Input parameters
const vtkm::Id nSeeds = 25;
const vtkm::Id nSteps = 2000;
const vtkm::Float32 tStep = 0.5f;
const vtkm::Id direction = vtkm::worklet::internal::BOTH;
// Point location of vertices from a CastAndCall but needs a static cast eventually
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float32, 3> > vertexArray;
// OpenGL display variables
Quaternion qrot;
int lastx, lasty;
int mouse_state = 1;
//
// Functor to retrieve vertex locations from the CoordinateSystem
// Actually need a static cast to ArrayHandle from DynamicArrayHandleCoordinateSystem
// but haven't been able to figure out what that is
//
struct GetVertexArray
{
template <typename ArrayHandleType>
VTKM_CONT_EXPORT
void operator()(ArrayHandleType array) const
{
this->GetVertexPortal(array.GetPortalConstControl());
}
private:
template <typename PortalType>
VTKM_CONT_EXPORT
void GetVertexPortal(const PortalType &portal) const
{
for (vtkm::Id index = 0; index < portal.GetNumberOfValues(); index++)
{
vertexArray.GetPortalControl().Set(index, portal.Get(index));
}
}
};
//
// Initialize the OpenGL state
//
void initializeGL()
{
glClearColor(0.0f, 0.0f, 0.0f, 0.0f);
glEnable(GL_DEPTH_TEST);
glShadeModel(GL_SMOOTH);
float white[] = { 0.8f, 0.8f, 0.8f, 1.0f };
float black[] = { 0.0f, 0.0f, 0.0f, 1.0f };
float lightPos[] = { 10.0f, 10.0f, 10.5f, 1.0f };
glLightfv(GL_LIGHT0, GL_AMBIENT, white);
glLightfv(GL_LIGHT0, GL_DIFFUSE, white);
glLightfv(GL_LIGHT0, GL_SPECULAR, black);
glLightfv(GL_LIGHT0, GL_POSITION, lightPos);
glLightModeli(GL_LIGHT_MODEL_TWO_SIDE, 1);
glEnable(GL_LIGHTING);
glEnable(GL_LIGHT0);
glEnable(GL_NORMALIZE);
glColorMaterial(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE);
glEnable(GL_COLOR_MATERIAL);
}
//
// Render the output using simple OpenGL
//
void displayCall()
{
glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
glEnable(GL_DEPTH_TEST);
glMatrixMode(GL_PROJECTION);
glLoadIdentity();
gluPerspective( 60.0f, 1.0f, 1.0f, 100.0f);
glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
gluLookAt(0.0f, 0.0f, 100.0f, 0.0f, 0.0f, 0.0f, 0.0f, 1.0f, 0.0f);
glLineWidth(1.0f);
glPushMatrix();
float rotationMatrix[16];
qrot.getRotMat(rotationMatrix);
glMultMatrixf(rotationMatrix);
glTranslatef(-0.5f, -0.5f, -0.5f);
// Get the cell set, coordinate system and coordinate data
vtkm::cont::CellSetExplicit<> &cellSet =
outDataSet.GetCellSet(0).CastTo<vtkm::cont::CellSetExplicit<> >();
const vtkm::cont::DynamicArrayHandleCoordinateSystem &coordArray =
outDataSet.GetCoordinateSystem(0).GetData();
vtkm::Id numberOfCells = cellSet.GetNumberOfCells();
vtkm::Id numberOfPoints = coordArray.GetNumberOfValues();
// Need the actual vertex points from a static cast of the dynamic array but can't get it right
// So use cast and call on a functor that stores that dynamic array into static array we created
vertexArray.Allocate(numberOfPoints);
coordArray.CastAndCall(GetVertexArray());
// Each cell is a polyline
glColor3f(1.0f, 0.0f, 0.0f);
for (vtkm::Id polyline = 0; polyline < numberOfCells; polyline++)
{
vtkm::Vec<vtkm::Id, nSteps> polylineIndices;
vtkm::IdComponent numIndices = cellSet.GetNumberOfPointsInCell(polyline);
cellSet.GetIndices(polyline, polylineIndices);
glBegin(GL_LINE_STRIP);
for (vtkm::IdComponent i = 0; i < numIndices; i++)
{
vtkm::Vec<vtkm::Float32,3> pt = vertexArray.GetPortalConstControl().Get(polylineIndices[i]);
glVertex3f(pt[0], pt[1], pt[2]);
}
glEnd();
}
glPopMatrix();
glutSwapBuffers();
}
// Allow rotations of the view
void mouseMove(int x, int y)
{
vtkm::Float32 dx = static_cast<vtkm::Float32>(x - lastx);
vtkm::Float32 dy = static_cast<vtkm::Float32>(y - lasty);
if (mouse_state == 0)
{
vtkm::Float32 pi = static_cast<float>(vtkm::Pi());
Quaternion newRotX;
newRotX.setEulerAngles(-0.2f * dx * pi / 180.0f, 0.0f, 0.0f);
qrot.mul(newRotX);
Quaternion newRotY;
newRotY.setEulerAngles(0.0f, 0.0f, -0.2f * dy * pi / 180.0f);
qrot.mul(newRotY);
}
lastx = x;
lasty = y;
glutPostRedisplay();
}
// Respond to mouse button
void mouseCall(int button, int state, int x, int y)
{
if (button == 0) mouse_state = state;
if ((button == 0) && (state == 0)) { lastx = x; lasty = y; }
}
namespace {
template <typename T>
VTKM_EXEC_CONT_EXPORT
vtkm::Vec<T,3> Normalize(vtkm::Vec<T,3> v)
{
T magnitude = static_cast<T>(sqrt(vtkm::dot(v, v)));
T zero = static_cast<T>(0.0);
T one = static_cast<T>(1.0);
if (magnitude == zero)
return vtkm::make_Vec(zero, zero, zero);
else
return one / magnitude * v;
}
}
// Run streamlines on a uniform grid of vector data
int main(int argc, char* argv[])
{
std::cout << "StreamLineUniformGrid Example" << std::endl;
std::cout << "Parameters are fileName [numSeeds maxSteps timeStep direction]" << std::endl << std::endl;
std::cout << "Direction is FORWARD=0 BACKWARD=1 BOTH=2" << std::endl << std::endl;
std::cout << "File is expected to be binary with xdim ydim zdim as 32 bit integers " << std::endl;
std::cout << "followed by vector data per dimension point as 32 bit float" << std::endl;
// Read in the vector data for testing
FILE * pFile = fopen(argv[1], "rb");
if (pFile == NULL) perror ("Error opening file");
// Size of the dataset
int dims[3];
fread(dims, sizeof(int), 3, pFile);
const vtkm::Id3 vdims(dims[0], dims[1], dims[2]);
// Read vector data at each point of the uniform grid and store
vtkm::Id nElements = vdims[0] * vdims[1] * vdims[2] * 3;
float* data = new float[nElements];
fread(data, sizeof(float), nElements, pFile);
std::vector<vtkm::Vec<vtkm::Float32, 3> > field;
for (vtkm::Id i = 0; i < nElements; i++)
{
vtkm::Float32 x = data[i];
vtkm::Float32 y = data[++i];
vtkm::Float32 z = data[++i];
vtkm::Vec<vtkm::Float32, 3> vecData(x, y, z);
field.push_back(Normalize(vecData));
}
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float32, 3> > fieldArray;
fieldArray = vtkm::cont::make_ArrayHandle(&field[0], field.size());
// Construct the input dataset (uniform) to hold the input and set vector data
vtkm::cont::DataSet inDataSet;
vtkm::cont::ArrayHandleUniformPointCoordinates coordinates(vdims);
inDataSet.AddCoordinateSystem(vtkm::cont::CoordinateSystem("coordinates", 1, coordinates));
inDataSet.AddField(vtkm::cont::Field("vecData", 1, vtkm::cont::Field::ASSOC_POINTS, fieldArray));
vtkm::cont::CellSetStructured<3> inCellSet("cells");
inCellSet.SetPointDimensions(vtkm::make_Vec(vdims[0], vdims[1], vdims[2]));
inDataSet.AddCellSet(inCellSet);
// Create and run the filter
streamLineFilter = new vtkm::worklet::StreamLineFilterUniformGrid<vtkm::Float32, DeviceAdapter>();
outDataSet = streamLineFilter->Run(inDataSet,
direction,
nSeeds,
nSteps,
tStep);
// Render the output dataset of polylines
lastx = lasty = 0;
glutInit(&argc, argv);
glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE | GLUT_DEPTH);
glutInitWindowSize(1000, 1000);
glutCreateWindow("VTK-m Uniform 3D StreamLines");
initializeGL();
glutDisplayFunc(displayCall);
glutMotionFunc(mouseMove);
glutMouseFunc(mouseCall);
glutMainLoop();
return 0;
}
#if (defined(VTKM_GCC) || defined(VTKM_CLANG)) && !defined(VTKM_PGI)
# pragma GCC diagnostic pop
#endif

23
examples/streamline/StreamLineUniformGridTBB.cxx 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_TBB
#include "StreamLineUniformGrid.cxx"

BIN
examples/streamline/tornado.vec Normal file
View File

Binary file not shown.

1
vtkm/worklet/CMakeLists.txt
View File

@ -36,6 +36,7 @@ set(headers
ScatterCounting.h
ScatterIdentity.h
ScatterUniform.h
StreamLineUniformGrid.h
TetrahedralizeExplicitGrid.h
TetrahedralizeUniformGrid.h
Threshold.h

440
vtkm/worklet/StreamLineUniformGrid.h Normal file
View File

@ -0,0 +1,440 @@
//============================================================================
// 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_StreamLineUniformGrid_h
#define vtk_m_worklet_StreamLineUniformGrid_h
#include <vtkm/cont/DeviceAdapter.h>
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/ArrayHandleCounting.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/CellSetStructured.h>
#include <vtkm/cont/CellSetExplicit.h>
#include <vtkm/cont/Field.h>
#include <vtkm/worklet/DispatcherMapField.h>
#include <vtkm/worklet/ScatterUniform.h>
#include <vtkm/worklet/WorkletMapField.h>
#include <vtkm/exec/ExecutionWholeArray.h>
namespace vtkm {
// Take this out when defined in CellShape.h
const vtkm::UInt8 CELL_SHAPE_POLY_LINE = 4;
namespace worklet {
namespace internal {
enum StreamLineMode
{
FORWARD = 0,
BACKWARD = 1,
BOTH = 2
};
// Trilinear interpolation to calculate vector data at position
template <typename FieldType, typename PortalType>
VTKM_EXEC_EXPORT
vtkm::Vec<FieldType, 3> VecDataAtPos(
vtkm::Vec<FieldType, 3> pos,
const vtkm::Id3 &vdims,
const vtkm::Id &planesize,
const vtkm::Id &rowsize,
const PortalType &vecdata)
{
// Adjust initial position to be within bounding box of grid
for (vtkm::IdComponent d = 0; d < 3; d++)
{
if (pos[d] < 0.0f)
pos[d] = 0.0f;
if (pos[d] > static_cast<FieldType>(vdims[d] - 1))
pos[d] = static_cast<FieldType>(vdims[d] - 1);
}
// Set the eight corner indices with no wraparound
vtkm::Id3 idx000, idx001, idx010, idx011, idx100, idx101, idx110, idx111;
idx000[0] = static_cast<vtkm::Id>(floor(pos[0]));
idx000[1] = static_cast<vtkm::Id>(floor(pos[1]));
idx000[2] = static_cast<vtkm::Id>(floor(pos[2]));
idx001 = idx000; idx001[0] = (idx001[0] + 1) <= vdims[0] - 1 ? idx001[0] + 1 : vdims[0] - 1;
idx010 = idx000; idx010[1] = (idx010[1] + 1) <= vdims[1] - 1 ? idx010[1] + 1 : vdims[1] - 1;
idx011 = idx010; idx011[0] = (idx011[0] + 1) <= vdims[0] - 1 ? idx011[0] + 1 : vdims[0] - 1;
idx100 = idx000; idx100[2] = (idx100[2] + 1) <= vdims[2] - 1 ? idx100[2] + 1 : vdims[2] - 1;
idx101 = idx100; idx101[0] = (idx101[0] + 1) <= vdims[0] - 1 ? idx101[0] + 1 : vdims[0] - 1;
idx110 = idx100; idx110[1] = (idx110[1] + 1) <= vdims[1] - 1 ? idx110[1] + 1 : vdims[1] - 1;
idx111 = idx110; idx111[0] = (idx111[0] + 1) <= vdims[0] - 1 ? idx111[0] + 1 : vdims[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
vtkm::Vec<FieldType, 3> v;
a = pos[2] - static_cast<FieldType>(floor(pos[2]));
v[0] = (1.0f - a) * v0[0] + v1[0];
v[1] = (1.0f - a) * v0[1] + v1[1];
v[2] = (1.0f - a) * v0[2] + v1[2];
return v;
}
}
/// \brief Compute the streamline
template <typename FieldType, typename DeviceAdapter>
class StreamLineFilterUniformGrid
{
public:
struct IsUnity
{
template<typename T>
VTKM_EXEC_CONT_EXPORT bool operator()(const T &x) const
{
return x == T(1);
}
};
typedef vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3> > FieldHandle;
typedef typename FieldHandle::template ExecutionTypes<DeviceAdapter>::PortalConst FieldPortalConstType;
class MakeStreamLines : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(FieldIn<IdType> seedId,
FieldIn<> position,
ExecObject numIndices,
ExecObject validPoint,
ExecObject streamLines);
typedef void ExecutionSignature(_1, _2, _3, _4, _5, VisitIndex);
typedef _1 InputDomain;
typedef vtkm::worklet::ScatterUniform ScatterType;
VTKM_CONT_EXPORT
ScatterType GetScatter() const
{
return ScatterType(2);
}
FieldPortalConstType field;
const vtkm::Id3 vdims;
const vtkm::Id maxsteps;
const FieldType timestep;
const vtkm::Id planesize;
const vtkm::Id rowsize;
const vtkm::Id streammode;
VTKM_CONT_EXPORT
MakeStreamLines(const FieldType tStep,
const vtkm::Id sMode,
const vtkm::Id nSteps,
const vtkm::Id3 dims,
FieldPortalConstType fieldArray) :
timestep(tStep),
streammode(sMode),
maxsteps(nSteps),
vdims(dims),
planesize(dims[0] * dims[1]),
rowsize(dims[0]),
field(fieldArray)
{
}
VTKM_EXEC_EXPORT
void operator()(vtkm::Id &seedId,
vtkm::Vec<FieldType, 3> &seedPos,
vtkm::exec::ExecutionWholeArray<vtkm::IdComponent> &numIndices,
vtkm::exec::ExecutionWholeArray<vtkm::IdComponent> &validPoint,
vtkm::exec::ExecutionWholeArray<vtkm::Vec<FieldType, 3> > &slLists,
vtkm::IdComponent visitIndex) const
{
// Set initial offset into the output streams array
vtkm::Vec<FieldType, 3> pos = seedPos;
vtkm::Vec<FieldType, 3> pre_pos = seedPos;
// Forward tracing
if (visitIndex == 0 &&
(streammode == vtkm::worklet::internal::FORWARD ||
streammode == vtkm::worklet::internal::BOTH))
{
vtkm::Id index = (seedId * 2) * maxsteps;
bool done = false;
vtkm::Id step = 0;
validPoint.Set(index, 1);
slLists.Set(index++, pos);
while (done != true && step < maxsteps)
{
vtkm::Vec<FieldType, 3> vdata, adata, bdata, cdata, ddata;
vdata = internal::VecDataAtPos<FieldType, FieldPortalConstType>
(pos, vdims, planesize, rowsize, field);
for (vtkm::IdComponent d = 0; d < 3; d++)
{
adata[d] = timestep * vdata[d];
pos[d] += adata[d] / 2.0f;
}
vdata = internal::VecDataAtPos<FieldType, FieldPortalConstType>
(pos, vdims, planesize, rowsize, field);
for (vtkm::IdComponent d = 0; d < 3; d++)
{
bdata[d] = timestep * vdata[d];
pos[d] += bdata[d] / 2.0f;
}
vdata = internal::VecDataAtPos<FieldType, FieldPortalConstType>
(pos, vdims, planesize, rowsize, field);
for (vtkm::IdComponent d = 0; d < 3; d++)
{
cdata[d] = timestep * vdata[d];
pos[d] += cdata[d] / 2.0f;
}
vdata = internal::VecDataAtPos<FieldType, FieldPortalConstType>
(pos, vdims, planesize, rowsize, field);
for (vtkm::IdComponent d = 0; d < 3; d++)
{
ddata[d] = timestep * vdata[d];
pos[d] += (adata[d] + (2.0f * bdata[d]) + (2.0f * cdata[d]) + ddata[d]) / 6.0f;
}
if (pos[0] < 0.0f || pos[0] > vdims[0] ||
pos[1] < 0.0f || pos[1] > vdims[1] ||
pos[2] < 0.0f || pos[2] > vdims[2])
{
pos = pre_pos;
done = true;
} else {
validPoint.Set(index, 1);
slLists.Set(index++, pos);
pre_pos = pos;
}
step++;
}
numIndices.Set(seedId * 2, static_cast<vtkm::IdComponent>(step));
}
// Backward tracing
if (visitIndex == 1 &&
(streammode == vtkm::worklet::internal::BACKWARD ||
streammode == vtkm::worklet::internal::BOTH))
{
vtkm::Id index = (seedId * 2 + 1) * maxsteps;
bool done = false;
vtkm::Id step = 0;
validPoint.Set(index, 1);
slLists.Set(index++, pos);
while (done != true && step < maxsteps)
{
vtkm::Vec<FieldType, 3> vdata, adata, bdata, cdata, ddata;
vdata = internal::VecDataAtPos<FieldType, FieldPortalConstType>
(pos, vdims, planesize, rowsize, field);
for (vtkm::IdComponent d = 0; d < 3; d++)
{
adata[d] = timestep * (0.0f - vdata[d]);
pos[d] += adata[d] / 2.0f;
}
vdata = internal::VecDataAtPos<FieldType, FieldPortalConstType>
(pos, vdims, planesize, rowsize, field);
for (vtkm::IdComponent d = 0; d < 3; d++)
{
bdata[d] = timestep * (0.0f - vdata[d]);
pos[d] += bdata[d] / 2.0f;
}
vdata = internal::VecDataAtPos<FieldType, FieldPortalConstType>
(pos, vdims, planesize, rowsize, field);
for (vtkm::IdComponent d = 0; d < 3; d++)
{
cdata[d] = timestep * (0.0f - vdata[d]);
pos[d] += cdata[d] / 2.0f;
}
vdata = internal::VecDataAtPos<FieldType, FieldPortalConstType>
(pos, vdims, planesize, rowsize, field);
for (vtkm::IdComponent d = 0; d < 3; d++)
{
ddata[d] = timestep * (0.0f - vdata[d]);
pos[d] += (adata[d] + (2.0f * bdata[d]) + (2.0f * cdata[d]) + ddata[d]) / 6.0f;
}
if (pos[0] < 0.0f || pos[0] > vdims[0] ||
pos[1] < 0.0f || pos[1] > vdims[1] ||
pos[2] < 0.0f || pos[2] > vdims[2])
{
pos = pre_pos;
done = true;
} else {
validPoint.Set(index, 1);
slLists.Set(index++, pos);
pre_pos = pos;
}
step++;
}
numIndices.Set((seedId * 2) + 1, static_cast<vtkm::IdComponent>(step));
}
}
};
StreamLineFilterUniformGrid()
{
}
vtkm::cont::DataSet Run(const vtkm::cont::DataSet &InDataSet,
vtkm::Id streamMode,
vtkm::Id numSeeds,
vtkm::Id maxSteps,
FieldType timeStep)
{
typedef typename vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapter> DeviceAlgorithm;
// Get information from input dataset
vtkm::cont::CellSetStructured<3> &inCellSet =
InDataSet.GetCellSet(0).template CastTo<vtkm::cont::CellSetStructured<3> >();
vtkm::Id3 vdims= inCellSet.GetSchedulingRange(vtkm::TopologyElementTagPoint());
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3> > fieldArray =
InDataSet.GetField("vecData").GetData().
CastToArrayHandle<vtkm::Vec<FieldType, 3>, VTKM_DEFAULT_STORAGE_TAG>();
// Generate random seeds for starting streamlines
std::vector<vtkm::Vec<FieldType, 3> > seeds;
for (vtkm::Id i = 0; i < numSeeds; i++)
{
vtkm::Vec<FieldType, 3> seed;
seed[0] = static_cast<FieldType>(rand() % vdims[0]);
seed[1] = static_cast<FieldType>(rand() % vdims[1]);
seed[2] = static_cast<FieldType>(rand() % vdims[2]);
seeds.push_back(seed);
}
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3> > seedPosArray =
vtkm::cont::make_ArrayHandle(&seeds[0], numSeeds);
vtkm::cont::ArrayHandleCounting<vtkm::Id> seedIdArray(0, 1, numSeeds);
// Number of streams * number of steps * [forward, backward]
vtkm::Id numCells = numSeeds * 2;
vtkm::Id maxConnectivityLen = numCells * maxSteps;
// Stream array at max size will be filled with stream coordinates
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3> > streamArray;
streamArray.Allocate(maxConnectivityLen);
// NumIndices per polyline cell filled in by MakeStreamLines
vtkm::cont::ArrayHandle<vtkm::IdComponent> numIndices;
numIndices.Allocate(numCells);
// All cells are polylines
vtkm::cont::ArrayHandle<vtkm::UInt8> cellTypes;
cellTypes.Allocate(numCells);
vtkm::cont::ArrayHandleConstant<vtkm::UInt8> polyLineShape(vtkm::CELL_SHAPE_POLY_LINE, numCells);
DeviceAlgorithm::Copy(polyLineShape, cellTypes);
// Possible maxSteps points but if less use stencil
vtkm::cont::ArrayHandle<vtkm::IdComponent> validPoint;
vtkm::cont::ArrayHandleConstant<vtkm::Id> zeros(0, maxConnectivityLen);
validPoint.Allocate(maxConnectivityLen);
DeviceAlgorithm::Copy(zeros, validPoint);
// Worklet to make the streamlines
MakeStreamLines makeStreamLines(timeStep,
streamMode,
maxSteps,
vdims,
fieldArray.PrepareForInput(DeviceAdapter()));
typedef typename vtkm::worklet::DispatcherMapField<MakeStreamLines> MakeStreamLinesDispatcher;
MakeStreamLinesDispatcher makeStreamLinesDispatcher(makeStreamLines);
makeStreamLinesDispatcher.Invoke(
seedIdArray,
seedPosArray,
vtkm::exec::ExecutionWholeArray<vtkm::IdComponent>(numIndices, numCells),
vtkm::exec::ExecutionWholeArray<vtkm::IdComponent>(validPoint, maxConnectivityLen),
vtkm::exec::ExecutionWholeArray<vtkm::Vec<FieldType, 3> >(streamArray, maxConnectivityLen));
// Size of connectivity based on size of returned streamlines
vtkm::cont::ArrayHandle<vtkm::IdComponent> numIndicesOut;
vtkm::IdComponent connectivityLen = DeviceAlgorithm::ScanExclusive(numIndices, numIndicesOut);
// Connectivity is sequential
vtkm::cont::ArrayHandleCounting<vtkm::Id> connCount(0, 1, connectivityLen);
vtkm::cont::ArrayHandle<vtkm::Id> connectivity;
DeviceAlgorithm::Copy(connCount, connectivity);
// Compact the stream array so it only has valid points
vtkm::cont::ArrayHandle<vtkm::Vec<FieldType, 3> > coordinates;
DeviceAlgorithm::StreamCompact(streamArray,
validPoint,
coordinates,
IsUnity());
// Create the output data set
vtkm::cont::DataSet OutDataSet;
vtkm::cont::CellSetExplicit<> outCellSet;
outCellSet.Fill(cellTypes, numIndices, connectivity);
OutDataSet.AddCellSet(outCellSet);
OutDataSet.AddCoordinateSystem(vtkm::cont::CoordinateSystem("coordinates", 0, coordinates));
return OutDataSet;
}
};
}
}
#endif // vtk_m_worklet_StreamLineUniformGrid_h

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

@ -28,6 +28,7 @@ set(unit_tests
UnitTestPointElevation.cxx
UnitTestScatterCounting.cxx
UnitTestSplatKernels.cxx
UnitTestStreamLineUniformGrid.cxx
UnitTestTetrahedralizeExplicitGrid.cxx
UnitTestTetrahedralizeUniformGrid.cxx
UnitTestThreshold.cxx

241
vtkm/worklet/testing/UnitTestStreamLineUniformGrid.cxx Normal file
View File

@ -0,0 +1,241 @@
//============================================================================
// 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/worklet/StreamLineUniformGrid.h>
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/testing/Testing.h>
#include <fstream>
#include <vector>
#include <math.h>
namespace {
template <typename T>
VTKM_EXEC_CONT_EXPORT
vtkm::Vec<T,3> Normalize(vtkm::Vec<T,3> v)
{
T magnitude = static_cast<T>(sqrt(vtkm::dot(v, v)));
T zero = static_cast<T>(0.0);
T one = static_cast<T>(1.0);
if (magnitude == zero)
return vtkm::make_Vec(zero, zero, zero);
else
return one / magnitude * v;
}
float data[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 TestStreamLineUniformGrid()
{
std::cout << "Testing StreamLineUniformGrid Filter" << std::endl;
typedef VTKM_DEFAULT_DEVICE_ADAPTER_TAG DeviceAdapter;
// Parameters for streamlines
vtkm::Id numSeeds = 5;
vtkm::Id maxSteps = 50;
vtkm::Float32 timeStep = 0.5f;
// Size of the dataset
const vtkm::Id3 vdims(5, 5, 5);
// Read vector data at each point of the uniform grid and store
vtkm::Id nElements = vdims[0] * vdims[1] * vdims[2] * 3;
std::vector<vtkm::Vec<vtkm::Float32, 3> > field;
for (vtkm::Id i = 0; i < nElements; i++)
{
vtkm::Float32 x = data[i];
vtkm::Float32 y = data[++i];
vtkm::Float32 z = data[++i];
vtkm::Vec<vtkm::Float32, 3> vecData(x, y, z);
field.push_back(Normalize(vecData));
}
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float32, 3> > fieldArray;
fieldArray = vtkm::cont::make_ArrayHandle(&field[0], nElements);
// Construct the input dataset (uniform) to hold the input and set vector data
vtkm::cont::DataSet inDataSet;
vtkm::cont::ArrayHandleUniformPointCoordinates coordinates(vdims);
inDataSet.AddCoordinateSystem(vtkm::cont::CoordinateSystem("coordinates", 1, coordinates));
inDataSet.AddField(vtkm::cont::Field("vecData", 1, vtkm::cont::Field::ASSOC_POINTS, fieldArray));
vtkm::cont::CellSetStructured<3> inCellSet("cells");
inCellSet.SetPointDimensions(vtkm::make_Vec(vdims[0], vdims[1], vdims[2]));
inDataSet.AddCellSet(inCellSet);
// Create and run the filter
vtkm::worklet::StreamLineFilterUniformGrid<vtkm::Float32, DeviceAdapter> streamLineFilter;
vtkm::cont::DataSet outDataSet = streamLineFilter.Run(inDataSet,
vtkm::worklet::internal::BOTH,
numSeeds,
maxSteps,
timeStep);
// Check output
vtkm::cont::CellSetExplicit<> &outCellSet =
outDataSet.GetCellSet(0).CastTo<vtkm::cont::CellSetExplicit<> >();
const vtkm::cont::DynamicArrayHandleCoordinateSystem &coordArray =
outDataSet.GetCoordinateSystem(0).GetData();
vtkm::Id numberOfCells = outCellSet.GetNumberOfCells();
vtkm::Id numberOfPoints = coordArray.GetNumberOfValues();
std::cout << "Number of polylines " << numberOfCells << std::endl;
std::cout << "Number of coordinates " << numberOfPoints << std::endl;
VTKM_TEST_ASSERT(test_equal(numberOfCells, numSeeds * 2),
"Wrong number of cells for stream lines");
}
int UnitTestStreamLineUniformGrid(int, char *[])
{
return vtkm::cont::testing::Testing::Run(TestStreamLineUniformGrid);
}