diff --git a/vtkm/CellShape.h b/vtkm/CellShape.h index 7dfe1c96a..e65afce05 100644 --- a/vtkm/CellShape.h +++ b/vtkm/CellShape.h @@ -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); \ diff --git a/vtkm/CellTraits.h b/vtkm/CellTraits.h index d7a90cb24..178515912 100644 --- a/vtkm/CellTraits.h +++ b/vtkm/CellTraits.h @@ -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); diff --git a/vtkm/exec/CellDerivative.h b/vtkm/exec/CellDerivative.h index ba3a872fd..575ecb484 100644 --- a/vtkm/exec/CellDerivative.h +++ b/vtkm/exec/CellDerivative.h @@ -550,6 +550,79 @@ VTKM_EXEC vtkm::Vec CellDerivative( static_cast((field[1] - field[0]) / wCoords.GetSpacing()[0]), T(0), T(0)); } +template +VTKM_EXEC vtkm::Vec CellDerivative( + const FieldVecType& field, + const WorldCoordType& wCoords, + const vtkm::Vec& pcoords, + vtkm::CellShapeTagPolyLine, + const vtkm::exec::FunctorBase& worklet) +{ + vtkm::IdComponent numPoints = field.GetNumberOfComponents(); + VTKM_ASSERT(numPoints >= 1); + VTKM_ASSERT(numPoints == wCoords.GetNumberOfComponents()); + + switch (numPoints) + { + case 1: + return CellDerivative(field, wCoords, pcoords, vtkm::CellShapeTagVertex(), worklet); + case 2: + return CellDerivative(field, wCoords, pcoords, vtkm::CellShapeTagLine(), worklet); + } + + using FieldType = typename FieldVecType::ComponentType; + using BaseComponentType = typename BaseComponent::Type; + + ParametricCoordType dt; + dt = static_cast(1) / static_cast(numPoints - 1); + vtkm::IdComponent idx = static_cast(vtkm::Ceil(pcoords[0] / dt)); + if (idx == 0) + idx = 1; + if (idx > numPoints - 1) + idx = numPoints - 1; + + FieldType deltaField(field[idx] - field[idx - 1]); + vtkm::Vec vec(wCoords[idx] - wCoords[idx - 1]); + + return detail::CellDerivativeLineImpl(deltaField, + vec, + vtkm::MagnitudeSquared(vec), + typename vtkm::TypeTraits::DimensionalityTag()); +} + +template +VTKM_EXEC vtkm::Vec CellDerivative( + const FieldVecType& field, + const vtkm::VecAxisAlignedPointCoordinates<1>& wCoords, + const vtkm::Vec& 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); + } + + ParametricCoordType dt; + dt = static_cast(1) / static_cast(numPoints - 1); + vtkm::IdComponent idx = static_cast(vtkm::Ceil(pcoords[0] / dt)); + if (idx == 0) + idx = 1; + if (idx > numPoints - 1) + idx = numPoints - 1; + + using T = typename FieldVecType::ComponentType; + + return vtkm::Vec( + static_cast((field[idx] - field[idx - 1]) / wCoords.GetSpacing()[0]), T(0), T(0)); +} + //----------------------------------------------------------------------------- namespace detail { diff --git a/vtkm/exec/CellEdge.h b/vtkm/exec/CellEdge.h index 1fec8c970..624dafd78 100644 --- a/vtkm/exec/CellEdge.h +++ b/vtkm/exec/CellEdge.h @@ -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 @@ -133,6 +133,15 @@ static inline VTKM_EXEC vtkm::IdComponent CellEdgeNumberOfEdges(vtkm::IdComponen return numPoints; } +static inline VTKM_EXEC vtkm::IdComponent CellEdgeNumberOfEdges(vtkm::IdComponent numPoints, + vtkm::CellShapeTagPolyLine, + const vtkm::exec::FunctorBase&) +{ + (void)numPoints; // Silence compiler warnings. + VTKM_ASSUME(numPoints > 0); + return detail::CellEdgeTables{}.NumEdges(vtkm::CELL_SHAPE_POLY_LINE); +} + static inline VTKM_EXEC vtkm::IdComponent CellEdgeNumberOfEdges( vtkm::IdComponent numPoints, vtkm::CellShapeTagGeneric shape, @@ -142,6 +151,10 @@ static inline VTKM_EXEC vtkm::IdComponent CellEdgeNumberOfEdges( { return CellEdgeNumberOfEdges(numPoints, vtkm::CellShapeTagPolygon(), worklet); } + else if (shape.Id == vtkm::CELL_SHAPE_POLY_LINE) + { + return CellEdgeNumberOfEdges(numPoints, vtkm::CellShapeTagPolyLine(), worklet); + } else { return detail::CellEdgeTables{}.NumEdges(shape.Id); diff --git a/vtkm/exec/CellFace.h b/vtkm/exec/CellFace.h index d5442e24a..7b1ccd95f 100644 --- a/vtkm/exec/CellFace.h +++ b/vtkm/exec/CellFace.h @@ -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 diff --git a/vtkm/exec/CellInside.h b/vtkm/exec/CellInside.h index 621315e24..99bc777de 100644 --- a/vtkm/exec/CellInside.h +++ b/vtkm/exec/CellInside.h @@ -36,6 +36,12 @@ static inline VTKM_EXEC bool CellInside(const vtkm::Vec& pcoords, vtkm::Ce return pcoords[0] >= T(0) && pcoords[0] <= T(1); } +template +static inline VTKM_EXEC bool CellInside(const vtkm::Vec& pcoords, vtkm::CellShapeTagPolyLine) +{ + return pcoords[0] >= T(0) && pcoords[0] <= T(1); +} + template static inline VTKM_EXEC bool CellInside(const vtkm::Vec& pcoords, vtkm::CellShapeTagTriangle) { diff --git a/vtkm/exec/CellInterpolate.h b/vtkm/exec/CellInterpolate.h index 1e6b0371c..45eb5f58d 100644 --- a/vtkm/exec/CellInterpolate.h +++ b/vtkm/exec/CellInterpolate.h @@ -197,6 +197,63 @@ VTKM_EXEC vtkm::Vec CellInterpolate( origin[0] + static_cast(pcoords[0]) * spacing[0], origin[1], origin[2]); } +//----------------------------------------------------------------------------- +template +VTKM_EXEC typename FieldVecType::ComponentType CellInterpolate( + const FieldVecType& field, + const vtkm::Vec& 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 = ParametricCoordType; + + T dt = 1 / static_cast(numPoints - 1); + vtkm::IdComponent idx = static_cast(pcoords[0] / dt); + if (idx == numPoints - 1) + return field[numPoints - 1]; + + T t = (pcoords[0] - static_cast(idx) * dt) / dt; + + return vtkm::Lerp(field[idx], field[idx + 1], t); +} + +template +VTKM_EXEC vtkm::Vec CellInterpolate( + const vtkm::VecAxisAlignedPointCoordinates<1>& field, + const vtkm::Vec& 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; + const T& origin = field.GetOrigin(); + const T& spacing = field.GetSpacing(); + + return T( + origin[0] + static_cast(pcoords[0]) * spacing[0], origin[1], origin[2]); +} + //----------------------------------------------------------------------------- template VTKM_EXEC typename FieldVecType::ComponentType CellInterpolate( diff --git a/vtkm/exec/ParametricCoordinates.h b/vtkm/exec/ParametricCoordinates.h index c9e6bbe6a..8a58c3d85 100644 --- a/vtkm/exec/ParametricCoordinates.h +++ b/vtkm/exec/ParametricCoordinates.h @@ -66,6 +66,26 @@ static inline VTKM_EXEC void ParametricCoordinatesCenter(vtkm::IdComponent numPo pcoords[2] = 0; } +template +static inline VTKM_EXEC void ParametricCoordinatesCenter(vtkm::IdComponent numPoints, + vtkm::Vec& pcoords, + vtkm::CellShapeTagPolyLine, + const vtkm::exec::FunctorBase& worklet) +{ + switch (numPoints) + { + case 1: + ParametricCoordinatesCenter(numPoints, pcoords, vtkm::CellShapeTagVertex(), worklet); + return; + case 2: + ParametricCoordinatesCenter(numPoints, pcoords, vtkm::CellShapeTagLine(), worklet); + return; + } + pcoords[0] = 0.5; + pcoords[1] = 0; + pcoords[2] = 0; +} + template static inline VTKM_EXEC void ParametricCoordinatesCenter(vtkm::IdComponent numPoints, vtkm::Vec& pcoords, @@ -249,6 +269,29 @@ static inline VTKM_EXEC void ParametricCoordinatesPoint(vtkm::IdComponent numPoi pcoords[2] = 0; } +template +static inline VTKM_EXEC void ParametricCoordinatesPoint(vtkm::IdComponent numPoints, + vtkm::IdComponent pointIndex, + vtkm::Vec& pcoords, + vtkm::CellShapeTagPolyLine, + const vtkm::exec::FunctorBase& functor) +{ + switch (numPoints) + { + case 1: + ParametricCoordinatesPoint( + numPoints, pointIndex, pcoords, vtkm::CellShapeTagVertex(), functor); + return; + case 2: + ParametricCoordinatesPoint(numPoints, pointIndex, pcoords, vtkm::CellShapeTagLine(), functor); + return; + } + pcoords[0] = + static_cast(pointIndex) / static_cast(numPoints - 1); + pcoords[1] = 0; + pcoords[2] = 0; +} + template static inline VTKM_EXEC void ParametricCoordinatesPoint(vtkm::IdComponent numPoints, vtkm::IdComponent pointIndex, @@ -634,6 +677,75 @@ WorldCoordinatesToParametricCoordinates(const WorldCoordVector& pointWCoords, return Vector3(numerator / denominator, 0, 0); } +template +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); + } + + using Vector3 = typename WorldCoordVector::ComponentType; + using T = typename Vector3::ComponentType; + + //Find the closest vertex to the point. + vtkm::IdComponent idx = 0; + Vector3 vec = pointWCoords[0] - wcoords; + T minDistSq = vtkm::Dot(vec, vec); + for (vtkm::IdComponent i = 1; i < numPoints; i++) + { + vec = pointWCoords[i] - wcoords; + T d = vtkm::Dot(vec, vec); + + if (d < minDistSq) + { + idx = i; + minDistSq = d; + } + } + + //Find the right segment, and the parameterization along that segment. + //Closest to 0, so segment is (0,1) + if (idx == 0) + idx = 1; + + //Find the pt projection onto the line segment at points idx and idx-1. + vec = pointWCoords[idx] - pointWCoords[idx - 1]; + T numerator = vtkm::Dot(vec, wcoords - pointWCoords[idx - 1]); + T denominator = vtkm::MagnitudeSquared(vec); + T segmentParam = numerator / denominator; + + //The point is on the OTHER side of idx. If there is a next segment reparam onto it. + if (segmentParam > 1 && idx < numPoints - 1) + { + idx = idx + 1; + vec = pointWCoords[idx] - pointWCoords[idx - 1]; + numerator = vtkm::Dot(vec, wcoords - pointWCoords[idx - 1]); + denominator = vtkm::MagnitudeSquared(vec); + segmentParam = numerator / denominator; + } + + //Segment param is [0,1] on that segment. + //Map that onto the param for the entire segment. + T dParam = static_cast(1) / static_cast(numPoints - 1); + T polyLineParam = static_cast(idx - 1) * dParam + segmentParam * dParam; + + return Vector3(polyLineParam, 0, 0); +} + template static inline VTKM_EXEC typename WorldCoordVector::ComponentType WorldCoordinatesToParametricCoordinates(const WorldCoordVector& pointWCoords, diff --git a/vtkm/exec/testing/UnitTestCellEdgeFace.cxx b/vtkm/exec/testing/UnitTestCellEdgeFace.cxx index 29a2e3468..928eb2199 100644 --- a/vtkm/exec/testing/UnitTestCellEdgeFace.cxx +++ b/vtkm/exec/testing/UnitTestCellEdgeFace.cxx @@ -215,6 +215,14 @@ struct TestCellFacesFunctor this->TryShapeWithNumPoints(vtkm::CellTraits::NUM_POINTS, CellShapeTag()); } + void operator()(vtkm::CellShapeTagPolyLine) const + { + for (vtkm::IdComponent numPoints = 3; numPoints < 7; numPoints++) + { + this->TryShapeWithNumPoints(numPoints, vtkm::CellShapeTagPolyLine()); + } + } + void operator()(vtkm::CellShapeTagPolygon) const { for (vtkm::IdComponent numPoints = 3; numPoints < 7; numPoints++) diff --git a/vtkm/exec/testing/UnitTestCellInterpolate.cxx b/vtkm/exec/testing/UnitTestCellInterpolate.cxx index a0eee7828..39852270a 100644 --- a/vtkm/exec/testing/UnitTestCellInterpolate.cxx +++ b/vtkm/exec/testing/UnitTestCellInterpolate.cxx @@ -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."); } diff --git a/vtkm/worklet/particleadvection/ParticleAdvectionWorklets.h b/vtkm/worklet/particleadvection/ParticleAdvectionWorklets.h index 15552c327..f63970f07 100644 --- a/vtkm/worklet/particleadvection/ParticleAdvectionWorklets.h +++ b/vtkm/worklet/particleadvection/ParticleAdvectionWorklets.h @@ -216,7 +216,8 @@ private: vtkm::cont::ArrayHandle cellTypes; cellTypes.Allocate(numSeeds); - vtkm::cont::ArrayHandleConstant polyLineShape(vtkm::CELL_SHAPE_LINE, numSeeds); + vtkm::cont::ArrayHandleConstant polyLineShape(vtkm::CELL_SHAPE_POLY_LINE, + numSeeds); vtkm::cont::ArrayCopy(polyLineShape, cellTypes); vtkm::cont::ArrayHandle cellCounts; diff --git a/vtkm/worklet/testing/UnitTestParticleAdvection.cxx b/vtkm/worklet/testing/UnitTestParticleAdvection.cxx index c77f371a6..924a83e15 100644 --- a/vtkm/worklet/testing/UnitTestParticleAdvection.cxx +++ b/vtkm/worklet/testing/UnitTestParticleAdvection.cxx @@ -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()