Implementation for tube filter.

This commit is contained in:
Dave Pugmire 2019-06-21 09:45:37 -04:00
parent 03e7b77fc9
commit 92348e28fd
2 changed files with 400 additions and 651 deletions

@ -13,10 +13,8 @@
#include <typeinfo>
#include <vtkm/VectorAnalysis.h>
#include <vtkm/cont/Algorithm.h>
#include <vtkm/cont/ArrayHandleCounting.h>
#include <vtkm/cont/CellSetPermutation.h>
#include <vtkm/cont/CellSetExplicit.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/rendering/raytracing/MeshConnectivityBuilder.h>
#include <vtkm/worklet/DispatcherMapTopology.h>
#include <vtkm/worklet/WorkletMapField.h>
#include <vtkm/worklet/WorkletMapTopology.h>
@ -31,40 +29,47 @@ class CountSegments : public vtkm::worklet::WorkletMapPointToCell
{
public:
VTKM_CONT
CountSegments(vtkm::Id& n)
CountSegments(const vtkm::Id& n)
: NumSides(n)
, NumVertsPerCell(4) //Quad
{
}
using ControlSignature = void(CellSetIn, FieldOut, FieldOut);
using ExecutionSignature = void(CellShape, PointCount, PointIndices, _2, _3);
using ControlSignature = void(CellSetIn,
FieldOut ptsPerPolyline,
FieldOut ptsPerTube,
FieldOut numTubeCells);
using ExecutionSignature = void(CellShape shapeType,
PointCount numPoints,
_2 ptsPerPolyline,
_3 ptsPerTube,
_4 numTubeCells);
using InputDomain = _1;
template <typename CellShapeTag, typename PointIndexType>
template <typename CellShapeTag>
VTKM_EXEC void operator()(const CellShapeTag& shapeType,
const vtkm::IdComponent& numPoints,
const PointIndexType& ptIndices,
vtkm::Id& ptsPerSegment,
vtkm::Id& totalPoints) const
vtkm::Id& ptsPerPolyline,
vtkm::Id& ptsPerTube,
vtkm::Id& numTubeCells) const
{
if (shapeType.Id == vtkm::CELL_SHAPE_POLY_LINE)
{
ptsPerSegment = numPoints;
totalPoints = this->NumSides * numPoints;
std::cout << "Pt indices: ";
for (int i = 0; i < numPoints; i++)
std::cout << ptIndices[i] << " ";
std::cout << std::endl;
ptsPerPolyline = numPoints;
ptsPerTube = this->NumSides * numPoints;
numTubeCells = this->NumVertsPerCell * (numPoints - 1) * this->NumSides;
}
else
{
ptsPerSegment = 0;
totalPoints = 0;
ptsPerPolyline = 0;
ptsPerTube = 0;
numTubeCells = 0;
}
}
private:
vtkm::Id NumSides;
vtkm::Id NumVertsPerCell;
};
class GenerateNormals : public vtkm::worklet::WorkletMapPointToCell
@ -80,16 +85,35 @@ public:
using ControlSignature = void(CellSetIn cellset,
WholeArrayIn pointCoords,
FieldInTo pointOffset,
FieldInTo polylineOffset,
WholeArrayOut newNormals);
using ExecutionSignature = void(CellShape shapeType,
PointCount numPoints,
PointIndices ptIndices,
_2 inPts,
_3 pointOffset,
_3 polylineOffset,
_4 outNormals);
using InputDomain = _1;
template <typename InPointsType, typename PointIndexType>
VTKM_EXEC vtkm::Id FindValidSegment(const InPointsType& inPts,
const PointIndexType& ptIndices,
const vtkm::Id& numPoints,
vtkm::Id start) const
{
auto ps = inPts.Get(ptIndices[start]);
vtkm::Id end = start + 1;
while (end < numPoints)
{
auto pe = inPts.Get(ptIndices[end]);
if (vtkm::Magnitude(pe - ps) > 0)
return end - 1;
}
return numPoints;
}
template <typename CellShapeTag,
typename PointIndexType,
typename InPointsType,
@ -98,7 +122,7 @@ public:
const vtkm::IdComponent& numPoints,
const PointIndexType& ptIndices,
const InPointsType& inPts,
const vtkm::Id& ptOffset,
const vtkm::Id& polylineOffset,
OutNormalType& outNormals) const
{
//Ignore non-polyline and 0 point polylines.
@ -113,33 +137,45 @@ public:
//The following follows the VTK implementation in:
//vtkPolyLine::GenerateSlidingNormals
vtkm::Vec<vtkm::FloatDefault, 3> sPrev, sNext, normal, p0, p1;
vtkm::IdComponent sNextId = FindValidSegment(inPts, ptIndices, numPoints, 0);
//Find the first normal.
// Find two non-parallel segments and compute the normal.
vtkm::Vec<vtkm::FloatDefault, 3> sPrev, sNext, normal;
vtkm::IdComponent sNextId = 0;
auto p0 = inPts.Get(ptIndices[sNextId]);
auto p1 = inPts.Get(ptIndices[sNextId + 1]);
std::cout << sNextId << ": P01= " << p0 << " " << p1 << std::endl;
if (sNextId != numPoints) // at least one valid segment
{
p0 = inPts.Get(ptIndices[sNextId]);
p1 = inPts.Get(ptIndices[sNextId + 1]);
sPrev = vtkm::Normal(p1 - p0);
p0 = p1;
sNextId = 1;
while (++sNextId < numPoints)
}
else // no valid segments. Set everything to the default normal.
{
p1 = inPts.Get(ptIndices[sNextId]);
std::cout << sNextId << ": P01= " << p0 << " " << p1 << std::endl;
sNext = vtkm::Normal(p1 - p0);
std::cout << " pn = " << sPrev << " " << sNext << std::endl;
auto n = vtkm::Cross(sPrev, sNext);
if (vtkm::Magnitude(n) > vecMagnitudeEps)
{
sPrev = sNext;
normal = n;
break;
for (vtkm::Id i = 0; i < numPoints; i++)
outNormals.Set(polylineOffset + i, this->DefaultNorm);
return;
}
p0 = p1;
// find the next valid, non-parallel segment
while (++sNextId < numPoints)
{
sNextId = FindValidSegment(inPts, ptIndices, numPoints, sNextId);
if (sNextId != numPoints)
{
p0 = inPts.Get(ptIndices[sNextId]);
p1 = inPts.Get(ptIndices[sNextId + 1]);
sNext = vtkm::Normal(p1 - p0);
// now the starting normal should simply be the cross product
// in the following if statement we check for the case where
// the two segments are parallel, in which case, continue searching
// for the next valid segment
auto n = vtkm::Cross(sPrev, sNext);
//std::cout<<" "<<sNextId<<"n= "<<n<<" "<<sPrev<<" x "<<sNext<<std::endl;
if (vtkm::Magnitude(n) > vecMagnitudeEps)
{
normal = n;
sPrev = sNext;
break;
}
}
}
//only one valid segment...
@ -153,33 +189,45 @@ public:
normal[j] = -sPrev[(j + 1) % 3] / sPrev[j];
break;
}
std::cout << "Only one valid segment.... " << normal << std::endl;
}
vtkm::Normalize(normal);
std::cout << "Normals: np= " << numPoints << " N= " << normal << std::endl;
std::cout << "ptIndices= {";
for (int i = 0; i < numPoints; i++)
std::cout << ptIndices[i] << " ";
std::cout << "}" << std::endl;
std::cout << "ptOffset= " << ptOffset << std::endl;
/*
std::cout<<"Normals: np= "<<numPoints<<" N= "<<normal<<std::endl;
std::cout<<"ptIndices= {";
for (int i = 0; i < numPoints; i++) std::cout<<ptIndices[i]<<" ";
std::cout<<"}"<<std::endl;
std::cout<<"polylineOffset= "<<polylineOffset<<std::endl;
*/
//std::cout<<"LAST while loop: sNextId= "<<sNextId<<" normal= "<<normal<<std::endl;
vtkm::Id lastNormalId = 0;
while (++sNextId < numPoints)
{
if (sNextId == numPoints - 1)
sNextId = FindValidSegment(inPts, ptIndices, numPoints, sNextId);
if (sNextId == numPoints)
break;
p0 = inPts.Get(ptIndices[sNextId]);
p1 = inPts.Get(ptIndices[sNextId + 1]);
sNext = vtkm::Normal(p1 - p0);
auto w = vtkm::Cross(sPrev, normal);
//std::cout<<std::endl;
//std::cout<<"while_"<<sNextId<<" p0= "<<p0<<" p1= "<<p1<<std::endl;
//std::cout<<"while_"<<sNextId<<" sPrev= "<<sPrev<<" "<<"sNext= "<<sNext<<" w= "<<w<<std::endl;
if (vtkm::Magnitude(w) == 0) //can't use this segment
continue;
vtkm::Normalize(w);
auto q = vtkm::Cross(sNext, sPrev);
//std::cout<<"q= "<<q<<" "<<sNext<<" x "<<sPrev<<std::endl;
if (vtkm::Magnitude(q) == 0) //can't use this segment
continue;
vtkm::Normalize(q);
//std::cout<<"normal= "<<normal<<" sprev= "<<sPrev<<" sNext= "<<sNext<<" w= "<<w<<" q= "<<q<<" || "<<vtkm::Magnitude(q)<<std::endl;
vtkm::FloatDefault f1 = vtkm::Dot(q, normal);
vtkm::FloatDefault f2 = 1 - (f1 * f1);
if (f2 > 0)
@ -188,29 +236,29 @@ public:
f2 = 0;
auto c = vtkm::Normal(sNext + sPrev);
//std::cout<<"c= "<<c<<" "<<sNext<<" x "<<sPrev<<std::endl;
w = vtkm::Cross(c, q);
c = vtkm::Cross(sPrev, q);
if ((vtkm::Dot(normal, c) * vtkm::Dot(w, c)) < 0)
f2 = -f2;
//std::cout<<"round0 update: "<<lastNormalId<<" "<<sNextId<<std::endl;
for (vtkm::Id i = lastNormalId; i < sNextId; i++)
{
std::cout << "round0: " << ptOffset + i << " " << normal << std::endl;
outNormals.Set(ptOffset + i, normal);
//std::cout<<"round0: "<<polylineOffset+i<<" "<<normal<<std::endl;
outNormals.Set(polylineOffset + i, normal);
}
lastNormalId = sNextId;
sPrev = sNext;
normal = (f1 * q) + (f2 * w);
}
for (vtkm::Id i = lastNormalId; i < numPoints; i++)
{
std::cout << "round1: " << ptOffset + i << " " << normal << std::endl;
outNormals.Set(ptOffset + i, normal);
//std::cout<<"stuff: "<<c<<" "<<w<<" "<<q<<" "<<f1<<" "<<f2<<" normal= "<<normal<<std::endl;
}
/*
for (int i = 0; i < numPoints; i++)
outNormals.Set(ptOffset+i, vtkm::Vec<vtkm::FloatDefault, 3>(ptOffset, i, 0));
*/
for (vtkm::Id i = lastNormalId; i < numPoints; i++)
{
//std::cout<<"round1: "<<polylineOffset+i<<" "<<normal<<std::endl;
outNormals.Set(polylineOffset + i, normal);
}
}
private:
@ -231,15 +279,17 @@ public:
using ControlSignature = void(CellSetIn cellset,
WholeArrayIn pointCoords,
WholeArrayIn normals,
FieldInTo pointOffset,
FieldInTo tubePointOffsets,
FieldInTo polylineOffset,
WholeArrayOut newPointCoords);
using ExecutionSignature = void(CellShape shapeType,
PointCount numPoints,
PointIndices ptIndices,
_2 inPts,
_3 inNormals,
_4 ptOffset,
_5 outPts);
_4 tubePointOffsets,
_5 polylineOffset,
_6 outPts);
using InputDomain = _1;
template <typename CellShapeTag,
@ -252,16 +302,17 @@ public:
const PointIndexType& ptIndices,
const InPointsType& inPts,
const InNormalsType& inNormals,
const vtkm::Id& ptOffset,
const vtkm::Id& tubePointOffsets,
const vtkm::Id& polylineOffset,
OutPointsType& outPts) const
{
// std::cout<<std::endl;
// std::cout<<std::endl;
// std::cout<<"GeneratePoints: offset="<<tubePointOffsets<<" n="<<numPoints<<std::endl;
if (shapeType.Id != vtkm::CELL_SHAPE_POLY_LINE)
return;
//std::cout<<"inPts.sz() "<<inPts.GetNumberOfValues()<<std::endl;
//std::cout<<"outPts.sz() "<<outPts.GetNumberOfValues()<<std::endl;
std::cout << "ptOffset= " << ptOffset << std::endl;
vtkm::Vec<vtkm::FloatDefault, 3> n, p, pNext, sNext, sPrev, startCapNorm, endCapNorm;
for (vtkm::IdComponent j = 0; j < numPoints; j++)
{
@ -288,53 +339,43 @@ public:
sPrev = sNext;
sNext = pNext - p;
}
n = inNormals.Get(ptOffset + j);
n = inNormals.Get(polylineOffset + j);
if (vtkm::Magnitude(sNext) == 0)
throw vtkm::cont::ErrorBadValue("Coincident points.");
vtkm::Normalize(sNext);
auto s = (sPrev + sNext) / 2.;
if (vtkm::Magnitude(s) == 0)
s = vtkm::Cross(sPrev, n);
vtkm::Normalize(s);
auto w = vtkm::Cross(s, n);
if (vtkm::Magnitude(w) == 0)
throw vtkm::cont::ErrorBadValue("Bad normal in Tube worklet.");
vtkm::Normalize(w);
//create orthogonal coordinate system.
auto nP = vtkm::Normal(vtkm::Cross(w, s));
auto nP = vtkm::Cross(w, s);
vtkm::Normalize(nP);
//std::cout<<"normal_"<<j<<" "<<n<<std::endl;
//std::cout<<"Vectors_"<<j<<" s= "<<s<<" w= "<<w<<" nP= "<<nP<<std::endl;
//this only implements the 'sides share vertices' line 476
vtkm::Vec<vtkm::FloatDefault, 3> normal;
vtkm::Id outIdx = tubePointOffsets + this->NumSides * j;
for (vtkm::IdComponent k = 0; k < this->NumSides; k++)
{
vtkm::FloatDefault cosValue = vtkm::Cos(static_cast<vtkm::FloatDefault>(k * this->Theta));
vtkm::FloatDefault sinValue = vtkm::Sin(static_cast<vtkm::FloatDefault>(k * this->Theta));
vtkm::FloatDefault angle = static_cast<vtkm::FloatDefault>(k) * this->Theta;
vtkm::FloatDefault cosValue = vtkm::Cos(angle);
vtkm::FloatDefault sinValue = vtkm::Sin(angle);
normal = w * cosValue + nP * sinValue;
auto newPt = p + this->Radius * normal;
//std::cout<<"outPts: "<<outIdx+k<<" = "<<newPt<<std::endl;
outPts.Set(outIdx + k, newPt);
//std::cout<<" outPt["<<outIdx+k<<"] = "<<newPt<<std::endl;
}
}
#if 0
for (vtkm::IdComponent i = 0; i < numPoints; i++)
{
vtkm::Id pidx = ptIndices[i];
std::cout<<"GeneratePoints: "<<pidx<<" --> ";
auto pt = inPts.Get(pidx);
vtkm::Id outIdx = ptOffset + this->NumSides*i;
for (vtkm::Id j = 0; j < this->NumSides; j++)
{
std::cout<<outIdx+j<<" ";
pt[0] = -1;
pt[1] = -1;
pt[2] = -1;
outPts.Set(outIdx+j, pt);
}
std::cout<<std::endl;
}
#endif
}
private:
@ -343,428 +384,65 @@ private:
vtkm::FloatDefault Theta;
};
#if 0
class GeneratePoints : public vtkm::worklet::WorkletMapPointToCell
{
public:
using ControlSignature = void(CellSetIn cellSet, FieldOut);
/*
WholeArrayIn pointCoords,
FieldInTo cellOffset,
FieldInTo ptOffset,
WholeArrayOut pts);
*/
// using ExecutionSignature = void(CellShape, PointCount, PointIndices); //, _2, _3, _4, _5);
using ExecutionSignature = void(CellShape, PointCount, PointIndices, _2);
using InputDomain = _1;
class GenerateCells : public vtkm::worklet::WorkletMapPointToCell
{
public:
VTKM_CONT
GeneratePoints(const vtkm::Id& n, const vtkm::FloatDefault& r)
GenerateCells(const vtkm::Id& n)
: NumSides(n)
, Radius(r)
{
}
// template <typename CellShapeTag, typename PointIndexType, typename PointPortal>
template <typename CellShapeTag, typename PointIndexType>
VTKM_EXEC
void operator()(const CellShapeTag &shapeType,
using ControlSignature = void(CellSetIn cellset,
FieldInTo tubePointOffsets,
FieldInTo tubeCellOffsets,
WholeArrayOut outConnectivity);
using ExecutionSignature = void(CellShape shapeType,
PointCount numPoints,
_2 tubePointOffset,
_3 tubeCellOffsets,
_4 outConn);
using InputDomain = _1;
template <typename CellShapeTag, typename OutConnType>
VTKM_EXEC void operator()(const CellShapeTag& shapeType,
const vtkm::IdComponent& numPoints,
const PointIndexType& ptIndices,
vtkm::Id &idx) const
// const vtkm::Id& cellOffset,
// const vtkm::Id& ptOffset,
// const PointPortal& inPointsPortal,
// PointPortal& outPointsPortal) const
const vtkm::Id& tubePointOffset,
const vtkm::Id& tubeCellOffset,
OutConnType& outConn) const
{
if (shapeType.Id != vtkm::CELL_SHAPE_POLY_LINE)
return;
// std::cout<<std::endl<<std::endl;
// std::cout<<"GenerateCells: ptOffset= "<<ptOffset<<" connOffset "<<connOffset<<std::endl;
// std::cout<<" pointIndices: ";
// for (int i = 0; i < numPoints; i++) std::cout<<ptIndices[i]<<" "; std::cout<<std::endl;
vtkm::Id outIdx = tubeCellOffset;
for (vtkm::IdComponent i = 0; i < numPoints - 1; i++)
{
for (vtkm::Id j = 0; j < this->NumSides; j++)
{
outConn.Set(outIdx + 0, tubePointOffset + i * this->NumSides + j);
outConn.Set(outIdx + 1, tubePointOffset + i * this->NumSides + (j + 1) % this->NumSides);
outConn.Set(outIdx + 2,
tubePointOffset + (i + 1) * this->NumSides + (j + 1) % this->NumSides);
outConn.Set(outIdx + 3, tubePointOffset + (i + 1) * this->NumSides + j);
// for (int k = 0; k < 4; k++)
// std::cout<<" outConn["<<outIdx+k<<"] = "<<outConn.Get(outIdx+k)<<std::endl;
outIdx += 4;
}
}
// std::cout<<"*********** gen cells done ********"<<std::endl;
}
private:
vtkm::Id NumSides;
vtkm::FloatDefault Radius;
};
#endif
#if 0
class MeowMeow : public vtkm::worklet::WorkletMapCellToPoint
{
public:
VTKM_CONT
MeowMeow() {}
using ControlSignature = void(CellSetIn, WholeArrayIn);
using ExecutionSignature = void(InputIndex pointIndex,
CellIndices incidentCells,
_1 cellSet,
_2 pointCoordsPortal);
using InputDomain = _1;
template <typename IncidentCellVecType, typename PointType>
VTKM_EXEC
void operator()(const vtkm::Id &idx,
const IncidentCellVecType& indices,
const PointType &pt) const
{
vtkm::IdComponent ncells = indices.GetNumberOfComponents();
}
};
class GeneratePoints : public vtkm::worklet::WorkletMapField
{
public:
using ControlSignature = void(FieldIn, FieldOut);
using ExecutionSignature = void(InputIndex, _1, _2);
using InputDomain = _1;
VTKM_CONT
GeneratePoints(vtkm::FloatDefault rad, vtkm::Id numSides)
: Radius(rad)
, NumSides(numSides)
{
}
VTKM_CONT
GeneratePoints()
{
}
template <typename T>
VTKM_EXEC
void operator()(const vtkm::Id& idx,
const vtkm::Vec<T,3>& pt,
vtkm::Vec<T,3>& newPt) const
{
newPt[0] = pt[0];
}
private:
vtkm::FloatDefault Radius;
vtkm::FloatDefault NumSides;
};
#endif
private:
vtkm::IdComponent NumSides;
};
}
#define SEG_PER_TRI 3
//CSS is CellSetStructured
#define TRI_PER_CSS 12
class Tube
{
public:
template <int DIM>
class SegmentedStructured : public vtkm::worklet::WorkletMapPointToCell
{
public:
typedef void ControlSignature(CellSetIn cellset, FieldInTo, WholeArrayOut);
typedef void ExecutionSignature(FromIndices, _2, _3);
//typedef _1 InputDomain;
VTKM_CONT
SegmentedStructured() {}
#if defined(VTKM_MSVC)
#pragma warning(push)
#pragma warning(disable : 4127) //conditional expression is constant
#endif
template <typename CellNodeVecType, typename OutIndicesPortal>
VTKM_EXEC void cell2seg(vtkm::Id3 idx,
vtkm::Vec<Id, 3>& segment,
const vtkm::Id offset,
const CellNodeVecType& cellIndices,
OutIndicesPortal& outputIndices) const
{
segment[1] = cellIndices[vtkm::IdComponent(idx[0])];
segment[2] = cellIndices[vtkm::IdComponent(idx[1])];
outputIndices.Set(offset, segment);
segment[1] = cellIndices[vtkm::IdComponent(idx[1])];
segment[2] = cellIndices[vtkm::IdComponent(idx[2])];
outputIndices.Set(offset + 1, segment);
segment[1] = cellIndices[vtkm::IdComponent(idx[2])];
segment[2] = cellIndices[vtkm::IdComponent(idx[0])];
outputIndices.Set(offset + 2, segment);
}
template <typename CellNodeVecType, typename OutIndicesPortal>
VTKM_EXEC void operator()(const CellNodeVecType& cellIndices,
const vtkm::Id& cellIndex,
OutIndicesPortal& outputIndices) const
{
if (DIM == 2)
{
// Do nothing mark says
}
else if (DIM == 3)
{
vtkm::Id offset = cellIndex * TRI_PER_CSS * SEG_PER_TRI;
vtkm::Vec<vtkm::Id, 3> segment;
segment[0] = cellIndex;
vtkm::Id3 idx;
idx[0] = 0;
idx[1] = 1;
idx[2] = 5;
cell2seg(idx, segment, offset, cellIndices, outputIndices);
idx[0] = 0;
idx[1] = 5;
idx[2] = 4;
offset += 3;
cell2seg(idx, segment, offset, cellIndices, outputIndices);
idx[0] = 1;
idx[1] = 2;
idx[2] = 6;
offset += 3;
cell2seg(idx, segment, offset, cellIndices, outputIndices);
idx[0] = 1;
idx[1] = 6;
idx[2] = 5;
offset += 3;
cell2seg(idx, segment, offset, cellIndices, outputIndices);
idx[0] = 3;
idx[1] = 7;
idx[2] = 6;
offset += 3;
cell2seg(idx, segment, offset, cellIndices, outputIndices);
idx[0] = 3;
idx[1] = 6;
idx[2] = 2;
offset += 3;
cell2seg(idx, segment, offset, cellIndices, outputIndices);
idx[0] = 0;
idx[1] = 4;
idx[2] = 7;
offset += 3;
cell2seg(idx, segment, offset, cellIndices, outputIndices);
idx[0] = 0;
idx[1] = 7;
idx[2] = 3;
offset += 3;
cell2seg(idx, segment, offset, cellIndices, outputIndices);
idx[0] = 0;
idx[1] = 3;
idx[2] = 2;
offset += 3;
cell2seg(idx, segment, offset, cellIndices, outputIndices);
idx[0] = 0;
idx[1] = 2;
idx[2] = 1;
offset += 3;
cell2seg(idx, segment, offset, cellIndices, outputIndices);
idx[0] = 4;
idx[1] = 5;
idx[2] = 6;
offset += 3;
cell2seg(idx, segment, offset, cellIndices, outputIndices);
idx[0] = 4;
idx[1] = 6;
idx[2] = 7;
offset += 3;
cell2seg(idx, segment, offset, cellIndices, outputIndices);
}
}
#if defined(VTKM_MSVC)
#pragma warning(pop)
#endif
};
class Cylinderize : public vtkm::worklet::WorkletMapPointToCell
{
public:
VTKM_CONT
Cylinderize() {}
typedef void ControlSignature(CellSetIn cellset, FieldInCell, WholeArrayOut);
typedef void ExecutionSignature(_2, CellShape, PointCount, PointIndices, WorkIndex, _3);
template <typename VecType, typename OutputPortal>
VTKM_EXEC void tri2seg(vtkm::Id& offset,
const VecType& cellIndices,
const vtkm::Id& cellId,
const vtkm::Id Id0,
const vtkm::Id Id1,
const vtkm::Id Id2,
OutputPortal& outputIndices) const
{
vtkm::Id3 segment;
segment[0] = cellId;
segment[1] = vtkm::Id(cellIndices[vtkm::IdComponent(Id0)]);
segment[2] = vtkm::Id(cellIndices[vtkm::IdComponent(Id1)]);
outputIndices.Set(offset++, segment);
segment[1] = vtkm::Id(cellIndices[vtkm::IdComponent(Id1)]);
segment[2] = vtkm::Id(cellIndices[vtkm::IdComponent(Id2)]);
outputIndices.Set(offset++, segment);
segment[1] = vtkm::Id(cellIndices[vtkm::IdComponent(Id2)]);
segment[2] = vtkm::Id(cellIndices[vtkm::IdComponent(Id0)]);
outputIndices.Set(offset++, segment);
}
template <typename VecType, typename OutputPortal>
VTKM_EXEC void operator()(const vtkm::Id& offset,
vtkm::CellShapeTagQuad shapeType,
const vtkm::IdComponent& numPoints,
const VecType& cellIndices,
const vtkm::Id& cellId,
OutputPortal& outputIndices) const
{
if (shapeType.Id == vtkm::CELL_SHAPE_QUAD)
{
vtkm::Id3 segment;
segment[0] = cellId;
segment[1] = cellIndices[0];
segment[2] = cellIndices[1];
outputIndices.Set(offset, segment);
segment[1] = cellIndices[1];
segment[2] = cellIndices[2];
outputIndices.Set(offset + 1, segment);
segment[1] = cellIndices[2];
segment[2] = cellIndices[3];
outputIndices.Set(offset + 2, segment);
segment[1] = cellIndices[3];
segment[2] = cellIndices[0];
outputIndices.Set(offset + 3, segment);
}
}
template <typename VecType, typename OutputPortal>
VTKM_EXEC void operator()(const vtkm::Id& pointOffset,
vtkm::CellShapeTagHexahedron vtkmNotUsed(shapeType),
const vtkm::IdComponent& vtkmNotUsed(numPoints),
const VecType& cellIndices,
const vtkm::Id& cellId,
OutputPortal& outputIndices) const
{
vtkm::Id offset = pointOffset;
tri2seg(offset, cellIndices, cellId, 0, 1, 5, outputIndices);
tri2seg(offset, cellIndices, cellId, 0, 5, 4, outputIndices);
tri2seg(offset, cellIndices, cellId, 1, 2, 6, outputIndices);
tri2seg(offset, cellIndices, cellId, 1, 6, 5, outputIndices);
tri2seg(offset, cellIndices, cellId, 3, 7, 6, outputIndices);
tri2seg(offset, cellIndices, cellId, 3, 6, 2, outputIndices);
tri2seg(offset, cellIndices, cellId, 0, 4, 7, outputIndices);
tri2seg(offset, cellIndices, cellId, 0, 7, 3, outputIndices);
tri2seg(offset, cellIndices, cellId, 0, 3, 2, outputIndices);
tri2seg(offset, cellIndices, cellId, 0, 2, 1, outputIndices);
tri2seg(offset, cellIndices, cellId, 4, 5, 6, outputIndices);
tri2seg(offset, cellIndices, cellId, 4, 6, 7, outputIndices);
}
template <typename VecType, typename OutputPortal>
VTKM_EXEC void operator()(const vtkm::Id& pointOffset,
vtkm::CellShapeTagGeneric shapeType,
const vtkm::IdComponent& numPoints,
const VecType& cellIndices,
const vtkm::Id& cellId,
OutputPortal& outputIndices) const
{
if (shapeType.Id == vtkm::CELL_SHAPE_LINE)
{
vtkm::Id3 segment;
segment[0] = cellId;
segment[1] = cellIndices[0];
segment[2] = cellIndices[1];
outputIndices.Set(pointOffset, segment);
}
else if (shapeType.Id == vtkm::CELL_SHAPE_POLY_LINE)
{
}
else if (shapeType.Id == vtkm::CELL_SHAPE_TRIANGLE)
{
vtkm::Id3 segment;
segment[0] = cellId;
segment[1] = cellIndices[0];
segment[2] = cellIndices[1];
outputIndices.Set(pointOffset, segment);
segment[1] = cellIndices[1];
segment[2] = cellIndices[2];
outputIndices.Set(pointOffset + 1, segment);
segment[1] = cellIndices[2];
segment[2] = cellIndices[0];
outputIndices.Set(pointOffset + 2, segment);
}
else if (shapeType.Id == vtkm::CELL_SHAPE_QUAD)
{
vtkm::Id3 segment;
segment[0] = cellId;
segment[1] = cellIndices[0];
segment[2] = cellIndices[1];
outputIndices.Set(pointOffset, segment);
segment[1] = cellIndices[1];
segment[2] = cellIndices[2];
outputIndices.Set(pointOffset + 1, segment);
segment[1] = cellIndices[2];
segment[2] = cellIndices[3];
outputIndices.Set(pointOffset + 2, segment);
segment[1] = cellIndices[3];
segment[2] = cellIndices[0];
outputIndices.Set(pointOffset + 3, segment);
}
else if (shapeType.Id == vtkm::CELL_SHAPE_TETRA)
{
vtkm::Id offset = pointOffset;
tri2seg(offset, cellIndices, cellId, 0, 3, 1, outputIndices);
tri2seg(offset, cellIndices, cellId, 1, 2, 3, outputIndices);
tri2seg(offset, cellIndices, cellId, 0, 2, 3, outputIndices);
tri2seg(offset, cellIndices, cellId, 0, 2, 1, outputIndices);
}
else if (shapeType.Id == vtkm::CELL_SHAPE_HEXAHEDRON)
{
vtkm::Id offset = pointOffset;
tri2seg(offset, cellIndices, cellId, 0, 1, 5, outputIndices);
tri2seg(offset, cellIndices, cellId, 0, 5, 4, outputIndices);
tri2seg(offset, cellIndices, cellId, 1, 2, 6, outputIndices);
tri2seg(offset, cellIndices, cellId, 1, 6, 5, outputIndices);
tri2seg(offset, cellIndices, cellId, 3, 7, 6, outputIndices);
tri2seg(offset, cellIndices, cellId, 3, 6, 2, outputIndices);
tri2seg(offset, cellIndices, cellId, 0, 4, 7, outputIndices);
tri2seg(offset, cellIndices, cellId, 0, 7, 3, outputIndices);
tri2seg(offset, cellIndices, cellId, 0, 3, 2, outputIndices);
tri2seg(offset, cellIndices, cellId, 0, 2, 1, outputIndices);
tri2seg(offset, cellIndices, cellId, 4, 5, 6, outputIndices);
tri2seg(offset, cellIndices, cellId, 4, 6, 7, outputIndices);
}
else if (shapeType.Id == vtkm::CELL_SHAPE_WEDGE)
{
vtkm::Id offset = pointOffset;
tri2seg(offset, cellIndices, cellId, 0, 1, 2, outputIndices);
tri2seg(offset, cellIndices, cellId, 3, 5, 4, outputIndices);
tri2seg(offset, cellIndices, cellId, 3, 0, 2, outputIndices);
tri2seg(offset, cellIndices, cellId, 3, 2, 5, outputIndices);
tri2seg(offset, cellIndices, cellId, 1, 4, 5, outputIndices);
tri2seg(offset, cellIndices, cellId, 1, 5, 2, outputIndices);
tri2seg(offset, cellIndices, cellId, 0, 3, 4, outputIndices);
tri2seg(offset, cellIndices, cellId, 0, 4, 1, outputIndices);
}
else if (shapeType.Id == vtkm::CELL_SHAPE_PYRAMID)
{
vtkm::Id offset = pointOffset;
tri2seg(offset, cellIndices, cellId, 0, 4, 1, outputIndices);
tri2seg(offset, cellIndices, cellId, 1, 2, 4, outputIndices);
tri2seg(offset, cellIndices, cellId, 2, 3, 4, outputIndices);
tri2seg(offset, cellIndices, cellId, 0, 4, 3, outputIndices);
tri2seg(offset, cellIndices, cellId, 3, 2, 1, outputIndices);
tri2seg(offset, cellIndices, cellId, 3, 1, 0, outputIndices);
}
}
}; //class cylinderize
public:
VTKM_CONT
Tube(const vtkm::Id& n, const vtkm::FloatDefault& r)
@ -775,74 +453,74 @@ public:
}
VTKM_CONT
void Run(const vtkm::cont::CoordinateSystem& coords, const vtkm::cont::DynamicCellSet& cellset)
void Run(const vtkm::cont::CoordinateSystem& coords,
const vtkm::cont::DynamicCellSet& cellset,
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::FloatDefault, 3>>& newPoints,
vtkm::cont::CellSetSingleType<>& newCells)
{
using ExplCoordsType = vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::FloatDefault, 3>>;
using NormalsType = vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::FloatDefault, 3>>;
if (coords.GetData().IsType<ExplCoordsType>() &&
if (!(coords.GetData().IsType<ExplCoordsType>() &&
(cellset.IsSameType(vtkm::cont::CellSetExplicit<>()) ||
cellset.IsSameType(vtkm::cont::CellSetSingleType<>())))
cellset.IsSameType(vtkm::cont::CellSetSingleType<>()))))
{
vtkm::cont::ArrayHandle<vtkm::Id> vertsPerCell, totalPointsPerCell;
throw vtkm::cont::ErrorBadValue("Tube filter only supported for polyline data.");
}
//Count number of polyline pts, tube pts and tube cells
vtkm::cont::ArrayHandle<vtkm::Id> ptsPerPolyline, ptsPerTube, numTubeCells;
detail::CountSegments countSegs(this->NumSides);
vtkm::worklet::DispatcherMapTopology<detail::CountSegments> countInvoker(countSegs);
countInvoker.Invoke(cellset, vertsPerCell, totalPointsPerCell);
countInvoker.Invoke(cellset, ptsPerPolyline, ptsPerTube, numTubeCells);
vtkm::Id totalVerts = vtkm::cont::Algorithm::Reduce(vertsPerCell, vtkm::Id(0));
if (totalVerts == 0)
vtkm::Id totalPolylinePts = vtkm::cont::Algorithm::Reduce(ptsPerPolyline, vtkm::Id(0));
if (totalPolylinePts == 0)
throw vtkm::cont::ErrorBadValue("Tube filter only supported for polyline data.");
vtkm::Id totalTubePts = vtkm::cont::Algorithm::Reduce(ptsPerTube, vtkm::Id(0));
vtkm::Id numPoints = vtkm::cont::Algorithm::Reduce(totalPointsPerCell, vtkm::Id(0));
ExplCoordsType inCoords = coords.GetData().Cast<ExplCoordsType>();
ExplCoordsType newPts;
newPts.Allocate(numPoints);
vtkm::cont::ArrayHandle<vtkm::Id> cellOffsets, pointOffsets;
vtkm::cont::Algorithm::ScanExclusive(vertsPerCell, cellOffsets);
vtkm::cont::Algorithm::ScanExclusive(totalPointsPerCell, pointOffsets);
NormalsType normals;
normals.Allocate(totalVerts);
vtkm::worklet::DispatcherMapTopology<detail::GenerateNormals> genNormalsDisp;
std::cout << std::endl << std::endl;
std::cout << "------------ GenerateNormals --------------------" << std::endl;
genNormalsDisp.Invoke(cellset, inCoords, cellOffsets, normals);
detail::GeneratePoints genPts(this->NumSides, this->Radius);
vtkm::worklet::DispatcherMapTopology<detail::GeneratePoints> genPtsDisp(genPts);
genPtsDisp.Invoke(cellset, inCoords, normals, pointOffsets, newPts);
vtkm::cont::ArrayHandle<vtkm::Id> polylineOffset, tubePointOffsets, tubeCellOffsets;
vtkm::cont::Algorithm::ScanExclusive(ptsPerPolyline, polylineOffset);
vtkm::cont::Algorithm::ScanExclusive(ptsPerTube, tubePointOffsets);
vtkm::cont::Algorithm::ScanExclusive(numTubeCells, tubeCellOffsets);
/*
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 3>> outputIndices;
outputIndices.Allocate(total);
std::cout<<"ptsPerPolyline: ";
vtkm::cont::printSummary_ArrayHandle(ptsPerPolyline, std::cout, true);
std::cout<<"totalPtsPerCell: ";
vtkm::cont::printSummary_ArrayHandle(ptsPerTube, std::cout, true);
std::cout<<"polylineOffset: ";
vtkm::cont::printSummary_ArrayHandle(polylineOffset, std::cout, true);
std::cout<<"tubePointOffsets: ";
vtkm::cont::printSummary_ArrayHandle(tubePointOffsets, std::cout, true);
std::cout<<"numTubeCells: ";
vtkm::cont::printSummary_ArrayHandle(numTubeCells, std::cout, true);
std::cout<<"tubeCellOffsets: ";
vtkm::cont::printSummary_ArrayHandle(tubeCellOffsets, std::cout, true);
std::cout<<"totalPolylinePts= "<<totalPolylinePts<<std::endl;
*/
vtkm::cont::printSummary_ArrayHandle(inCoords, std::cout, true);
vtkm::cont::printSummary_ArrayHandle(newPts, std::cout, true);
std::cout << "Normals: ";
vtkm::cont::printSummary_ArrayHandle(normals, std::cout, true);
vtkm::cont::printSummary_ArrayHandle(vertsPerCell, std::cout, true);
std::cout << "totalPtsPerCell: ";
vtkm::cont::printSummary_ArrayHandle(totalPointsPerCell, std::cout, true);
std::cout << "cellOffset: ";
vtkm::cont::printSummary_ArrayHandle(cellOffsets, std::cout, true);
std::cout << "pointOffsets: ";
vtkm::cont::printSummary_ArrayHandle(pointOffsets, std::cout, true);
std::cout << "totalVerts= " << totalVerts << std::endl;
//Generate normals at each point on all polylines
ExplCoordsType inCoords = coords.GetData().Cast<ExplCoordsType>();
NormalsType normals;
normals.Allocate(totalPolylinePts);
vtkm::worklet::DispatcherMapTopology<detail::GenerateNormals> genNormalsDisp;
genNormalsDisp.Invoke(cellset, inCoords, polylineOffset, normals);
//Generate the tube points
newPoints.Allocate(totalTubePts);
detail::GeneratePoints genPts(this->NumSides, this->Radius);
vtkm::worklet::DispatcherMapTopology<detail::GeneratePoints> genPtsDisp(genPts);
genPtsDisp.Invoke(cellset, inCoords, normals, tubePointOffsets, polylineOffset, newPoints);
//Generate tube cells
vtkm::cont::ArrayHandle<vtkm::Id> newConnectivity;
newConnectivity.Allocate(this->NumSides * (totalPolylinePts - 1) * 4);
detail::GenerateCells genCells(this->NumSides);
vtkm::worklet::DispatcherMapTopology<detail::GenerateCells> genCellsDisp(genCells);
genCellsDisp.Invoke(cellset, tubePointOffsets, tubeCellOffsets, newConnectivity);
std::cout << std::endl << std::endl << std::endl;
std::cout << "CellOffsets: " << std::endl;
vtkm::cont::printSummary_ArrayHandle(cellOffsets, std::cout, true);
std::cout << "VertsPerCell: " << std::endl;
vtkm::cont::printSummary_ArrayHandle(vertsPerCell, std::cout, true);
}
else
throw vtkm::cont::ErrorBadValue("Tube filter only supported for polyline data.");
newCells.Fill(totalTubePts, vtkm::CELL_SHAPE_QUAD, 4, newConnectivity);
}
private:

@ -10,6 +10,7 @@
#include <vtkm/cont/DataSetBuilderExplicit.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/io/writer/VTKDataSetWriter.h>
#include <vtkm/worklet/DispatcherMapField.h>
#include <vtkm/worklet/Tube.h>
@ -41,102 +42,172 @@ void appendPts(vtkm::cont::DataSetBuilderExplicitIterative& dsb,
ids.push_back(pid);
}
void TestTubeWorklets()
#if 0
void TestTubeWorklets2()
{
std::cout << "Testing Tube Worklet" << std::endl;
vtkm::cont::DataSetBuilderExplicitIterative dsb;
std::vector<vtkm::Id> ids;
ids.clear();
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault,3>(0,0,4), ids);
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault,3>(1,0,3), ids);
dsb.AddCell(vtkm::CELL_SHAPE_LINE, ids);
ids.clear();
vtkm::FloatDefault x0 = 0, x1 = 6.28, dx = 0.05;
//x1 = 6.28 / 6.;
for (vtkm::FloatDefault x = x0; x < x1; x += dx)
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault,3>(x,vtkm::Cos(x),.5*vtkm::Sin(x)), ids);
dsb.AddCell(vtkm::CELL_SHAPE_POLY_LINE, ids);
ids.clear();
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault,3>(0,0,4), ids);
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault,3>(1,0,3), ids);
dsb.AddCell(vtkm::CELL_SHAPE_POLY_LINE, ids);
ids.clear();
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault,3>(0,0,4), ids);
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault,3>(1,0,3), ids);
dsb.AddCell(vtkm::CELL_SHAPE_LINE, ids);
ids.clear();
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault,3>(0,0,4), ids);
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault,3>(1,0,3), ids);
dsb.AddCell(vtkm::CELL_SHAPE_LINE, ids);
ids.clear();
for (vtkm::FloatDefault x = x0; x < x1; x += dx)
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault,3>(x,2+vtkm::Cos(x),.5*vtkm::Sin(x)), ids);
dsb.AddCell(vtkm::CELL_SHAPE_POLY_LINE, ids);
vtkm::cont::DataSet ds = dsb.Create();
ds.PrintSummary(std::cout);
vtkm::worklet::Tube tubeWorklet(13, 0.05f);
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::FloatDefault,3>> newPoints;
vtkm::cont::CellSetSingleType<> newCells;
tubeWorklet.Run(ds.GetCoordinateSystem(0), ds.GetCellSet(0), newPoints, newCells);
vtkm::cont::DataSet ds2;
ds2.AddCoordinateSystem(vtkm::cont::CoordinateSystem("coords", newPoints));
ds2.AddCellSet(newCells);
ds2.PrintSummary(std::cout);
vtkm::io::writer::VTKDataSetWriter writer("tube.vtk");
writer.WriteDataSet(ds2);
if (1)
{
vtkm::io::writer::VTKDataSetWriter writer("poly.vtk");
writer.WriteDataSet(ds);
}
if (1)
{
dsb = vtkm::cont::DataSetBuilderExplicitIterative();
int nPts = newPoints.GetNumberOfValues();
ids.clear();
auto portal = newPoints.GetPortalControl();
for (int i = 0; i < nPts; i++)
{
vtkm::Id id = dsb.AddPoint(portal.Get(i));
dsb.AddCell(vtkm::CELL_SHAPE_VERTEX, {id});
}
ds = dsb.Create();
vtkm::io::writer::VTKDataSetWriter writer("pts.vtk");
writer.WriteDataSet(ds);
}
}
#endif
// Test with 3 polylines
void TestTube1()
{
vtkm::cont::DataSetBuilderExplicitIterative dsb;
std::vector<vtkm::Id> ids;
ids.clear();
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault, 3>(0, 0, 0), ids);
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault, 3>(1, 0, 0), ids);
dsb.AddCell(vtkm::CELL_SHAPE_POLY_LINE, ids);
ids.clear();
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault, 3>(0, 0, 0), ids);
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault, 3>(1, 0, 0), ids);
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault, 3>(2, 0, 0), ids);
dsb.AddCell(vtkm::CELL_SHAPE_POLY_LINE, ids);
#if 0
/*
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault,3>(3,1,0), ids);
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault,3>(4,0,0), ids);
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault,3>(5,0,0), ids);
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault,3>(6,0,0), ids);
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault,3>(7,0,0), ids);
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault,3>(8,0,0), ids);
*/
ids.clear();
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault, 3>(0, 0, 0), ids);
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault, 3>(1, 0, 0), ids);
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault, 3>(2, 1, 0), ids);
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault, 3>(3, 0, 0), ids);
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault, 3>(4, 0, 0), ids);
dsb.AddCell(vtkm::CELL_SHAPE_POLY_LINE, ids);
std::cout<<"PolyLine: ids= "<<ids<<std::endl;
ids.clear();
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault,3>(0,0,0), ids);
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault,3>(1,0,0), ids);
dsb.AddCell(vtkm::CELL_SHAPE_LINE, ids);
ids.clear();
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault,3>(0,0,0), ids);
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault,3>(1,0,0), ids);
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault,3>(2,0,0), ids);
dsb.AddCell(vtkm::CELL_SHAPE_POLY_LINE, ids);
ids.clear();
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault,3>(0,0,0), ids);
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault,3>(1,0,0), ids);
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault,3>(2,0,0), ids);
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault,3>(3,0,0), ids);
dsb.AddCell(vtkm::CELL_SHAPE_POLY_LINE, ids);
ids.clear();
#endif
#if 0
//0-5
dsb.AddPoint(vtkm::Vec<vtkm::FloatDefault,3>(0,0,0));
dsb.AddPoint(vtkm::Vec<vtkm::FloatDefault,3>(1,0,0));
dsb.AddPoint(vtkm::Vec<vtkm::FloatDefault,3>(2,1,0));
dsb.AddPoint(vtkm::Vec<vtkm::FloatDefault,3>(3,0,0));
dsb.AddPoint(vtkm::Vec<vtkm::FloatDefault,3>(4,1,0));
dsb.AddPoint(vtkm::Vec<vtkm::FloatDefault,3>(5,0,0));
dsb.AddCell(vtkm::CELL_SHAPE_POLY_LINE, {0,1,2,3,4,5});
//6-7
dsb.AddPoint(vtkm::Vec<vtkm::FloatDefault,3>(0,0,1));
dsb.AddPoint(vtkm::Vec<vtkm::FloatDefault,3>(0,0,2));
dsb.AddCell(vtkm::CELL_SHAPE_POLY_LINE, {6,7});
//8-9
dsb.AddPoint(vtkm::Vec<vtkm::FloatDefault,3>(0,1,1));
dsb.AddPoint(vtkm::Vec<vtkm::FloatDefault,3>(0,1,2));
dsb.AddCell(vtkm::CELL_SHAPE_LINE, {8,9});
//10-12
dsb.AddPoint(vtkm::Vec<vtkm::FloatDefault,3>(0,0,1));
dsb.AddPoint(vtkm::Vec<vtkm::FloatDefault,3>(1,0,1));
dsb.AddPoint(vtkm::Vec<vtkm::FloatDefault,3>(1,1,1));
dsb.AddCell(vtkm::CELL_SHAPE_TRIANGLE, {10,11,12});
//13-20
dsb.AddPoint(vtkm::Vec<vtkm::FloatDefault,3>(0,0,1));
dsb.AddPoint(vtkm::Vec<vtkm::FloatDefault,3>(1,0,1));
dsb.AddPoint(vtkm::Vec<vtkm::FloatDefault,3>(2,1,1));
dsb.AddPoint(vtkm::Vec<vtkm::FloatDefault,3>(3,0,1));
dsb.AddPoint(vtkm::Vec<vtkm::FloatDefault,3>(4,1,1));
dsb.AddPoint(vtkm::Vec<vtkm::FloatDefault,3>(5,0,1));
dsb.AddPoint(vtkm::Vec<vtkm::FloatDefault,3>(6,1,1));
int x = dsb.AddPoint(vtkm::Vec<vtkm::FloatDefault,3>(7,0,1));
std::cout<<"X= "<<x<<std::endl;
dsb.AddCell(vtkm::CELL_SHAPE_POLY_LINE, {13,14,15,16,17,18,19,20});
#endif
vtkm::cont::DataSet ds = dsb.Create();
ds.PrintSummary(std::cout);
vtkm::cont::CoordinateSystem pts;
vtkm::cont::CellSetSingleType<> cells("polyLines");
vtkm::worklet::Tube tubeWorklet(13, 0.05f);
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::FloatDefault, 3>> newPoints;
vtkm::cont::CellSetSingleType<> newCells;
tubeWorklet.Run(ds.GetCoordinateSystem(0), ds.GetCellSet(0), newPoints, newCells);
vtkm::worklet::Tube tubeWorklet(2, 0.1f);
tubeWorklet.Run(ds.GetCoordinateSystem(0), ds.GetCellSet(0));
VTKM_TEST_ASSERT(newPoints.GetNumberOfValues() == 130, "Wrong number of points in Tube worklet");
VTKM_TEST_ASSERT(newCells.GetCellShape(0) == vtkm::CELL_SHAPE_QUAD,
"Wrong cell shape in Tube worklet");
VTKM_TEST_ASSERT(newCells.GetNumberOfCells() == 117, "Wrong cell shape in Tube worklet");
}
// Test with 2 polylines and 1 triangle (which should be skipped).
void TestTube2()
{
vtkm::cont::DataSetBuilderExplicitIterative dsb;
std::vector<vtkm::Id> ids;
ids.clear();
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault, 3>(0, 0, 0), ids);
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault, 3>(1, 0, 0), ids);
dsb.AddCell(vtkm::CELL_SHAPE_POLY_LINE, ids);
ids.clear();
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault, 3>(0, 0, 0), ids);
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault, 3>(1, 0, 0), ids);
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault, 3>(1, 1, 0), ids);
dsb.AddCell(vtkm::CELL_SHAPE_TRIANGLE, ids);
ids.clear();
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault, 3>(0, 0, 0), ids);
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault, 3>(1, 0, 0), ids);
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault, 3>(2, 1, 0), ids);
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault, 3>(3, 0, 0), ids);
appendPts(dsb, vtkm::Vec<vtkm::FloatDefault, 3>(4, 0, 0), ids);
dsb.AddCell(vtkm::CELL_SHAPE_POLY_LINE, ids);
vtkm::cont::DataSet ds = dsb.Create();
vtkm::worklet::Tube tubeWorklet(13, 0.05f);
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::FloatDefault, 3>> newPoints;
vtkm::cont::CellSetSingleType<> newCells;
tubeWorklet.Run(ds.GetCoordinateSystem(0), ds.GetCellSet(0), newPoints, newCells);
VTKM_TEST_ASSERT(newPoints.GetNumberOfValues() == 91, "Wrong number of points in Tube worklet");
VTKM_TEST_ASSERT(newCells.GetCellShape(0) == vtkm::CELL_SHAPE_QUAD,
"Wrong cell shape in Tube worklet");
VTKM_TEST_ASSERT(newCells.GetNumberOfCells() == 78, "Wrong cell shape in Tube worklet");
}
void TestTubeWorklets()
{
std::cout << "Testing Tube Worklet" << std::endl;
TestTube1();
TestTube2();
}
}