Support for polylines.

This commit is contained in:
Dave Pugmire 2019-05-17 13:35:35 -04:00
parent be5b51fb01
commit 489995782f
12 changed files with 305 additions and 9 deletions

@ -26,7 +26,7 @@ enum CellShapeIdEnum
CELL_SHAPE_VERTEX = 1,
//CELL_SHAPE_POLY_VERTEX = 2,
CELL_SHAPE_LINE = 3,
//CELL_SHAPE_POLY_LINE = 4,
CELL_SHAPE_POLY_LINE = 4,
CELL_SHAPE_TRIANGLE = 5,
//CELL_SHAPE_TRIANGLE_STRIP = 6,
CELL_SHAPE_POLYGON = 7,
@ -110,7 +110,7 @@ VTKM_DEFINE_CELL_TAG(Empty, CELL_SHAPE_EMPTY);
VTKM_DEFINE_CELL_TAG(Vertex, CELL_SHAPE_VERTEX);
//VTKM_DEFINE_CELL_TAG(PolyVertex, CELL_SHAPE_POLY_VERTEX);
VTKM_DEFINE_CELL_TAG(Line, CELL_SHAPE_LINE);
//VTKM_DEFINE_CELL_TAG(PolyLine, CELL_SHAPE_POLY_LINE);
VTKM_DEFINE_CELL_TAG(PolyLine, CELL_SHAPE_POLY_LINE);
VTKM_DEFINE_CELL_TAG(Triangle, CELL_SHAPE_TRIANGLE);
//VTKM_DEFINE_CELL_TAG(TriangleStrip, CELL_SHAPE_TRIANGLE_STRIP);
VTKM_DEFINE_CELL_TAG(Polygon, CELL_SHAPE_POLYGON);
@ -182,6 +182,7 @@ struct CellShapeTagGeneric
vtkmGenericCellShapeMacroCase(CELL_SHAPE_EMPTY, call); \
vtkmGenericCellShapeMacroCase(CELL_SHAPE_VERTEX, call); \
vtkmGenericCellShapeMacroCase(CELL_SHAPE_LINE, call); \
vtkmGenericCellShapeMacroCase(CELL_SHAPE_POLY_LINE, call); \
vtkmGenericCellShapeMacroCase(CELL_SHAPE_TRIANGLE, call); \
vtkmGenericCellShapeMacroCase(CELL_SHAPE_POLYGON, call); \
vtkmGenericCellShapeMacroCase(CELL_SHAPE_QUAD, call); \

@ -104,7 +104,7 @@ VTKM_DEFINE_CELL_TRAITS(Empty, 0, 0);
VTKM_DEFINE_CELL_TRAITS(Vertex, 0, 1);
//VTKM_DEFINE_CELL_TRAITS_VARIABLE(PolyVertex, 0);
VTKM_DEFINE_CELL_TRAITS(Line, 1, 2);
//VTKM_DEFINE_CELL_TRAITS_VARIABLE(PolyLine, 1);
VTKM_DEFINE_CELL_TRAITS_VARIABLE(PolyLine, 1);
VTKM_DEFINE_CELL_TRAITS(Triangle, 2, 3);
//VTKM_DEFINE_CELL_TRAITS_VARIABLE(TriangleStrip, 2);
VTKM_DEFINE_CELL_TRAITS_VARIABLE(Polygon, 2);

@ -550,6 +550,82 @@ VTKM_EXEC vtkm::Vec<typename FieldVecType::ComponentType, 3> CellDerivative(
static_cast<T>((field[1] - field[0]) / wCoords.GetSpacing()[0]), T(0), T(0));
}
template <typename FieldVecType, typename WorldCoordType, typename ParametricCoordType>
VTKM_EXEC vtkm::Vec<typename FieldVecType::ComponentType, 3> CellDerivative(
const FieldVecType& field,
const WorldCoordType& wCoords,
const vtkm::Vec<ParametricCoordType, 3>& pcoords,
vtkm::CellShapeTagPolyLine,
const vtkm::exec::FunctorBase& worklet)
{
vtkm::IdComponent numPoints = field.GetNumberOfComponents();
VTKM_ASSERT(numPoints >= 1);
VTKM_ASSERT(numPoints == wCoords.GetNumberOfComponents());
/*
std::cout<<"CellDeriv_wWorld:"<<std::endl;
std::cout<<" pcoords= "<<pcoords<<std::endl;
std::cout<<" wCoords= (";
for (int i = 0; i < numPoints; i++) std::cout<<wCoords[i]<<" ";
std::cout<<")"<<std::endl;
std::cout<<" field= (";
for (int i = 0; i < numPoints; i++) std::cout<<field[i]<<" ";
std::cout<<")"<<std::endl;
*/
switch (numPoints)
{
case 1:
return CellDerivative(field, wCoords, pcoords, vtkm::CellShapeTagVertex(), worklet);
case 2:
return CellDerivative(field, wCoords, pcoords, vtkm::CellShapeTagLine(), worklet);
}
vtkm::FloatDefault dt = 1 / static_cast<vtkm::FloatDefault>(numPoints - 1);
vtkm::IdComponent idx = static_cast<vtkm::IdComponent>(pcoords[0] / dt);
if (idx == 0)
idx = 1;
using FieldType = typename FieldVecType::ComponentType;
using BaseComponentType = typename BaseComponent<FieldType>::Type;
FieldType deltaField(field[idx] - field[idx - 1]);
vtkm::Vec<BaseComponentType, 3> vec(wCoords[idx] - wCoords[idx - 1]);
return detail::CellDerivativeLineImpl(deltaField,
vec,
vtkm::MagnitudeSquared(vec),
typename vtkm::TypeTraits<FieldType>::DimensionalityTag());
}
template <typename FieldVecType, typename ParametricCoordType>
VTKM_EXEC vtkm::Vec<typename FieldVecType::ComponentType, 3> CellDerivative(
const FieldVecType& field,
const vtkm::VecAxisAlignedPointCoordinates<1>& wCoords,
const vtkm::Vec<ParametricCoordType, 3>& pcoords,
vtkm::CellShapeTagPolyLine,
const vtkm::exec::FunctorBase& worklet)
{
vtkm::IdComponent numPoints = field.GetNumberOfComponents();
VTKM_ASSERT(numPoints >= 1);
switch (numPoints)
{
case 1:
return CellDerivative(field, wCoords, pcoords, vtkm::CellShapeTagVertex(), worklet);
case 2:
return CellDerivative(field, wCoords, pcoords, vtkm::CellShapeTagLine(), worklet);
}
vtkm::FloatDefault dt = 1 / static_cast<vtkm::FloatDefault>(numPoints - 1);
vtkm::IdComponent idx = static_cast<vtkm::IdComponent>(pcoords[0] / dt);
if (idx == 0)
idx = 1;
using T = typename FieldVecType::ComponentType;
return vtkm::Vec<T, 3>(
static_cast<T>((field[idx] - field[idx - 1]) / wCoords.GetSpacing()[0]), T(0), T(0));
}
//-----------------------------------------------------------------------------
namespace detail
{

@ -39,7 +39,7 @@ public:
0, // 1: CELL_SHAPE_VERTEX
0, // 2: Unused
0, // 3: CELL_SHAPE_LINE
0, // 4: Unused
0, // 4: CELL_SHAPE_POLY_LINE
3, // 5: CELL_SHAPE_TRIANGLE
0, // 6: Unused
-1, // 7: CELL_SHAPE_POLYGON ---special case---
@ -73,7 +73,7 @@ public:
// 3: CELL_SHAPE_LINE
{ { -1, -1 }, { -1, -1 }, { -1, -1 }, { -1, -1 }, { -1, -1 }, { -1, -1 },
{ -1, -1 }, { -1, -1 }, { -1, -1 }, { -1, -1 }, { -1, -1 }, { -1, -1 } },
// 4: Unused
// 4: CELL_SHAPE_POLY_LINE
{ { -1, -1 }, { -1, -1 }, { -1, -1 }, { -1, -1 }, { -1, -1 }, { -1, -1 },
{ -1, -1 }, { -1, -1 }, { -1, -1 }, { -1, -1 }, { -1, -1 }, { -1, -1 } },
// 5: CELL_SHAPE_TRIANGLE

@ -38,7 +38,7 @@ public:
0, // 1: CELL_SHAPE_VERTEX
0, // 2: Unused
0, // 3: CELL_SHAPE_LINE
0, // 4: Unused
0, // 4: CELL_SHAPE_POLY_LINE
0, // 5: CELL_SHAPE_TRIANGLE
0, // 6: Unused
0, // 7: CELL_SHAPE_POLYGON
@ -62,7 +62,7 @@ public:
{ -1, -1, -1, -1, -1, -1 }, // 1: CELL_SHAPE_VERTEX
{ -1, -1, -1, -1, -1, -1 }, // 2: Unused
{ -1, -1, -1, -1, -1, -1 }, // 3: CELL_SHAPE_LINE
{ -1, -1, -1, -1, -1, -1 }, // 4: Unused
{ -1, -1, -1, -1, -1, -1 }, // 4: CELL_SHAPE_POLY_LINE
{ -1, -1, -1, -1, -1, -1 }, // 5: CELL_SHAPE_TRIANGLE
{ -1, -1, -1, -1, -1, -1 }, // 6: Unused
{ -1, -1, -1, -1, -1, -1 }, // 7: CELL_SHAPE_POLYGON
@ -98,7 +98,7 @@ public:
// 3: CELL_SHAPE_LINE
{ { -1, -1, -1, -1 }, { -1, -1, -1, -1 }, { -1, -1, -1, -1 },
{ -1, -1, -1, -1 }, { -1, -1, -1, -1 }, { -1, -1, -1, -1 } },
// 4: Unused
// 4: CELL_SHAPE_POLY_LINE
{ { -1, -1, -1, -1 }, { -1, -1, -1, -1 }, { -1, -1, -1, -1 },
{ -1, -1, -1, -1 }, { -1, -1, -1, -1 }, { -1, -1, -1, -1 } },
// 5: CELL_SHAPE_TRIANGLE

@ -36,6 +36,12 @@ static inline VTKM_EXEC bool CellInside(const vtkm::Vec<T, 3>& pcoords, vtkm::Ce
return pcoords[0] >= T(0) && pcoords[0] <= T(1);
}
template <typename T>
static inline VTKM_EXEC bool CellInside(const vtkm::Vec<T, 3>& pcoords, vtkm::CellShapeTagPolyLine)
{
return pcoords[0] >= T(0) && pcoords[0] <= T(1);
}
template <typename T>
static inline VTKM_EXEC bool CellInside(const vtkm::Vec<T, 3>& pcoords, vtkm::CellShapeTagTriangle)
{

@ -197,6 +197,83 @@ VTKM_EXEC vtkm::Vec<vtkm::FloatDefault, 3> CellInterpolate(
origin[0] + static_cast<vtkm::FloatDefault>(pcoords[0]) * spacing[0], origin[1], origin[2]);
}
//-----------------------------------------------------------------------------
template <typename FieldVecType, typename ParametricCoordType>
VTKM_EXEC typename FieldVecType::ComponentType CellInterpolate(
const FieldVecType& field,
const vtkm::Vec<ParametricCoordType, 3>& pcoords,
vtkm::CellShapeTagPolyLine,
const vtkm::exec::FunctorBase& worklet)
{
const vtkm::IdComponent numPoints = field.GetNumberOfComponents();
VTKM_ASSERT(numPoints >= 1);
/*
std::cout<<"PolyLine: numPts= "<<numPoints<<std::endl;
std::cout<<" field= [";
for (int i = 0; i < numPoints; i++)
std::cout<<field[i]<<" ";
std::cout<<"]"<<std::endl;
std::cout<<" nodeT= [";
if (numPoints > 1)
{
float dn = 1.0f/(float)(numPoints-1);
float t = 0;
for (int i = 0; i < numPoints; i++, t+=dn)
std::cout<<t<<" ";
std::cout<<"]"<<std::endl;
}
std::cout<<" pcoords= "<<pcoords<<std::endl;
*/
switch (numPoints)
{
case 1:
return CellInterpolate(field, pcoords, vtkm::CellShapeTagVertex(), worklet);
case 2:
return CellInterpolate(field, pcoords, vtkm::CellShapeTagLine(), worklet);
}
vtkm::FloatDefault dt = 1 / static_cast<vtkm::FloatDefault>(numPoints - 1);
vtkm::IdComponent idx = static_cast<vtkm::IdComponent>(pcoords[0] / dt);
if (idx == numPoints - 1)
return field[numPoints - 1];
vtkm::FloatDefault t = pcoords[0] - static_cast<vtkm::FloatDefault>(idx) * dt;
t = t / dt;
/*
std::cout<<" dt= "<<dt<<" idx= "<<idx<<" T= "<<t<<std::endl;
std::cout<<" Lerp("<<field[idx]<<", "<<field[idx+1]<<", "<<t<<")"<<std::endl;
*/
return vtkm::Lerp(field[idx], field[idx + 1], t);
}
template <typename ParametricCoordType>
VTKM_EXEC vtkm::Vec<vtkm::FloatDefault, 3> CellInterpolate(
const vtkm::VecAxisAlignedPointCoordinates<1>& field,
const vtkm::Vec<ParametricCoordType, 3>& pcoords,
vtkm::CellShapeTagPolyLine,
const vtkm::exec::FunctorBase& worklet)
{
const vtkm::IdComponent numPoints = field.GetNumberOfComponents();
VTKM_ASSERT(numPoints >= 1);
switch (numPoints)
{
case 1:
return CellInterpolate(field, pcoords, vtkm::CellShapeTagVertex(), worklet);
case 2:
return CellInterpolate(field, pcoords, vtkm::CellShapeTagLine(), worklet);
}
using T = vtkm::Vec<vtkm::FloatDefault, 3>;
const T& origin = field.GetOrigin();
const T& spacing = field.GetSpacing();
return T(
origin[0] + static_cast<vtkm::FloatDefault>(pcoords[0]) * spacing[0], origin[1], origin[2]);
}
//-----------------------------------------------------------------------------
template <typename FieldVecType, typename ParametricCoordType>
VTKM_EXEC typename FieldVecType::ComponentType CellInterpolate(

@ -66,6 +66,24 @@ static inline VTKM_EXEC void ParametricCoordinatesCenter(vtkm::IdComponent numPo
pcoords[2] = 0;
}
template <typename ParametricCoordType>
static inline VTKM_EXEC void ParametricCoordinatesCenter(vtkm::IdComponent numPoints,
vtkm::Vec<ParametricCoordType, 3>& pcoords,
vtkm::CellShapeTagPolyLine,
const vtkm::exec::FunctorBase& worklet)
{
switch (numPoints)
{
case 1:
return ParametricCoordinatesCenter(numPoints, pcoords, vtkm::CellShapeTagVertex(), worklet);
case 2:
return ParametricCoordinatesCenter(numPoints, pcoords, vtkm::CellShapeTagLine(), worklet);
}
pcoords[0] = 0.5;
pcoords[1] = 0;
pcoords[2] = 0;
}
template <typename ParametricCoordType>
static inline VTKM_EXEC void ParametricCoordinatesCenter(vtkm::IdComponent numPoints,
vtkm::Vec<ParametricCoordType, 3>& pcoords,
@ -249,6 +267,28 @@ static inline VTKM_EXEC void ParametricCoordinatesPoint(vtkm::IdComponent numPoi
pcoords[2] = 0;
}
template <typename ParametricCoordType>
static inline VTKM_EXEC void ParametricCoordinatesPoint(vtkm::IdComponent numPoints,
vtkm::IdComponent pointIndex,
vtkm::Vec<ParametricCoordType, 3>& pcoords,
vtkm::CellShapeTagPolyLine,
const vtkm::exec::FunctorBase& functor)
{
switch (numPoints)
{
case 1:
return ParametricCoordinatesPoint(
numPoints, pointIndex, pcoords, vtkm::CellShapeTagVertex(), functor);
case 2:
return ParametricCoordinatesPoint(
numPoints, pointIndex, pcoords, vtkm::CellShapeTagLine(), functor);
}
pcoords[0] =
static_cast<ParametricCoordType>(pointIndex) / static_cast<ParametricCoordType>(numPoints - 1);
pcoords[1] = 0;
pcoords[2] = 0;
}
template <typename ParametricCoordType>
static inline VTKM_EXEC void ParametricCoordinatesPoint(vtkm::IdComponent numPoints,
vtkm::IdComponent pointIndex,
@ -634,6 +674,90 @@ WorldCoordinatesToParametricCoordinates(const WorldCoordVector& pointWCoords,
return Vector3(numerator / denominator, 0, 0);
}
template <typename WorldCoordVector>
static inline VTKM_EXEC typename WorldCoordVector::ComponentType
WorldCoordinatesToParametricCoordinates(const WorldCoordVector& pointWCoords,
const typename WorldCoordVector::ComponentType& wcoords,
vtkm::CellShapeTagPolyLine,
bool& success,
const vtkm::exec::FunctorBase& worklet)
{
vtkm::IdComponent numPoints = pointWCoords.GetNumberOfComponents();
VTKM_ASSERT(pointWCoords.GetNumberOfComponents() >= 1);
switch (numPoints)
{
case 1:
return WorldCoordinatesToParametricCoordinates(
pointWCoords, wcoords, vtkm::CellShapeTagVertex(), success, worklet);
case 2:
return WorldCoordinatesToParametricCoordinates(
pointWCoords, wcoords, vtkm::CellShapeTagLine(), success, worklet);
}
vtkm::IdComponent idx = 0;
vtkm::FloatDefault minDistSq = vtkm::Dot(pointWCoords[0], wcoords);
for (vtkm::IdComponent i = 1; i < numPoints; i++)
{
vtkm::FloatDefault d = vtkm::Dot(pointWCoords[i], wcoords);
if (d < minDistSq)
{
idx = i;
minDistSq = d;
}
}
if (idx == 0)
idx = 1;
using Vector3 = typename WorldCoordVector::ComponentType;
using T = typename Vector3::ComponentType;
Vector3 vec = pointWCoords[idx] - pointWCoords[idx - 1];
T numerator = vtkm::Dot(vec, wcoords - pointWCoords[idx - 1]);
T denominator = vtkm::MagnitudeSquared(vec);
std::cout << "worldCoordsToParam:" << std::endl;
std::cout << " PointWcoords= (";
for (int i = 0; i < numPoints; i++)
std::cout << pointWCoords[i] << " ";
std::cout << ")" << std::endl;
std::cout << "wCoords= " << wcoords << std::endl;
std::cout << "idx= " << idx << std::endl;
std::cout << "Param= " << Vector3(numerator / denominator, 0, 0) << std::endl;
return Vector3(numerator / denominator, 0, 0);
#if 0
std::cout<<"worldCoordsToParam:"<<std::endl;
std::cout<<" PointWcoords= (";
for (int i = 0; i < numPoints; i++) std::cout<<pointWCoords[i]<<" ";
std::cout<<")"<<std::endl;
std::cout<<"wCoords= "<<wcoords<<std::endl;
success = true;
using Vector3 = typename WorldCoordVector::ComponentType;
return Vector3(0,0,0);
#endif
/*
vtkm::FloatDefault dt = 1/static_cast<vtkm::FloatDefault>(numPoints-1);
vtkm::IdComponent idx = 1; //static_cast<vtkm::IdComponent>(pcoords[0]/dt);
if (idx == 0)
idx = 1;
using Vector3 = typename WorldCoordVector::ComponentType;
using T = typename Vector3::ComponentType;
Vector3 vec = pointWCoords[idx] - pointWCoords[idx-1];
T numerator = vtkm::Dot(vec, wcoords - pointWCoords[idx-1]);
T denominator = vtkm::MagnitudeSquared(vec);
return Vector3(numerator / denominator, 0, 0);
*/
}
template <typename WorldCoordVector>
static inline VTKM_EXEC typename WorldCoordVector::ComponentType
WorldCoordinatesToParametricCoordinates(const WorldCoordVector& pointWCoords,

@ -215,6 +215,8 @@ struct TestCellFacesFunctor
this->TryShapeWithNumPoints(vtkm::CellTraits<CellShapeTag>::NUM_POINTS, CellShapeTag());
}
void operator()(vtkm::CellShapeTagPolyLine) const { throw "POLYLINE not supported yet"; }
void operator()(vtkm::CellShapeTagPolygon) const
{
for (vtkm::IdComponent numPoints = 3; numPoints < 7; numPoints++)

@ -94,6 +94,7 @@ struct TestInterpolateFunctor
vtkm::exec::CellInterpolate(fieldValues, pcoord, shape, workletProxy);
VTKM_TEST_ASSERT(!errorMessage.IsErrorRaised(), messageBuffer);
std::cout << "AVG= " << averageValue << " interp= " << interpolatedValue << std::endl;
VTKM_TEST_ASSERT(test_equal(averageValue, interpolatedValue),
"Interpolation at center not average value.");
}

@ -216,7 +216,8 @@ private:
vtkm::cont::ArrayHandle<vtkm::UInt8> cellTypes;
cellTypes.Allocate(numSeeds);
vtkm::cont::ArrayHandleConstant<vtkm::UInt8> polyLineShape(vtkm::CELL_SHAPE_LINE, numSeeds);
vtkm::cont::ArrayHandleConstant<vtkm::UInt8> polyLineShape(vtkm::CELL_SHAPE_POLY_LINE,
numSeeds);
vtkm::cont::ArrayCopy(polyLineShape, cellTypes);
vtkm::cont::ArrayHandle<vtkm::IdComponent> cellCounts;

@ -453,6 +453,14 @@ void ValidateStreamlineResult(const vtkm::worklet::StreamlineResult& res,
for (vtkm::Id i = 0; i < nSeeds; i++)
VTKM_TEST_ASSERT(res.stepsTaken.GetPortalConstControl().Get(i) <= maxSteps,
"Too many steps taken in streamline");
vtkm::cont::DataSet ds;
ds.AddCoordinateSystem(vtkm::cont::CoordinateSystem("coords", res.positions));
ds.AddCellSet(res.polyLines);
ds.PrintSummary(std::cout);
vtkm::io::writer::VTKDataSetWriter writer1("ds.vtk");
writer1.WriteDataSet(ds);
exit(0);
}
void TestParticleWorklets()