diff --git a/vtkm/CMakeLists.txt b/vtkm/CMakeLists.txt index 5e4c62dd1..e69821509 100644 --- a/vtkm/CMakeLists.txt +++ b/vtkm/CMakeLists.txt @@ -76,6 +76,8 @@ if(VTKm_ENABLE_LOGGING) endif() add_subdirectory(thirdparty/optionparser) add_subdirectory(thirdparty/taotuple) +add_subdirectory(thirdparty/vtkc) + add_subdirectory(testing) add_subdirectory(internal) diff --git a/vtkm/CellShape.h b/vtkm/CellShape.h index e65afce05..e055ea660 100644 --- a/vtkm/CellShape.h +++ b/vtkm/CellShape.h @@ -13,6 +13,17 @@ #include #include +#include +#include + +// Vtk-c does not have tags for Empty and PolyLine. Define dummy tags here to +// avoid compilation errors. These tags are not used anywhere. +namespace vtkc +{ +struct Empty; +struct PolyLine; +} + namespace vtkm { @@ -22,21 +33,21 @@ namespace vtkm enum CellShapeIdEnum { // Linear cells - CELL_SHAPE_EMPTY = 0, - CELL_SHAPE_VERTEX = 1, + CELL_SHAPE_EMPTY = vtkc::ShapeId::EMPTY, + CELL_SHAPE_VERTEX = vtkc::ShapeId::VERTEX, //CELL_SHAPE_POLY_VERTEX = 2, - CELL_SHAPE_LINE = 3, + CELL_SHAPE_LINE = vtkc::ShapeId::LINE, CELL_SHAPE_POLY_LINE = 4, - CELL_SHAPE_TRIANGLE = 5, + CELL_SHAPE_TRIANGLE = vtkc::ShapeId::TRIANGLE, //CELL_SHAPE_TRIANGLE_STRIP = 6, - CELL_SHAPE_POLYGON = 7, + CELL_SHAPE_POLYGON = vtkc::ShapeId::POLYGON, //CELL_SHAPE_PIXEL = 8, - CELL_SHAPE_QUAD = 9, - CELL_SHAPE_TETRA = 10, + CELL_SHAPE_QUAD = vtkc::ShapeId::QUAD, + CELL_SHAPE_TETRA = vtkc::ShapeId::TETRA, //CELL_SHAPE_VOXEL = 11, - CELL_SHAPE_HEXAHEDRON = 12, - CELL_SHAPE_WEDGE = 13, - CELL_SHAPE_PYRAMID = 14, + CELL_SHAPE_HEXAHEDRON = vtkc::ShapeId::HEXAHEDRON, + CELL_SHAPE_WEDGE = vtkc::ShapeId::WEDGE, + CELL_SHAPE_PYRAMID = vtkc::ShapeId::PYRAMID, NUMBER_OF_CELL_SHAPES }; @@ -58,6 +69,10 @@ struct CellShapeTagCheck : std::false_type { }; +/// Convert VTK-m tag to VTK-c tag +template +struct CellShapeTagVtkmToVtkc; + } // namespace internal /// Checks that the argument is a proper cell shape tag. This is a handy @@ -94,6 +109,11 @@ struct CellShapeIdToTag struct CellShapeTagCheck : std::true_type \ { \ }; \ + template <> \ + struct CellShapeTagVtkmToVtkc \ + { \ + using Type = vtkc::name; \ + }; \ } \ static inline VTKM_EXEC_CONT const char* GetCellShapeName(vtkm::CellShapeTag##name) \ { \ @@ -140,6 +160,35 @@ struct CellShapeTagGeneric vtkm::UInt8 Id; }; +namespace internal +{ + +template +VTKM_EXEC_CONT inline typename CellShapeTagVtkmToVtkc::Type make_VtkcCellShapeTag( + const VtkmCellShapeTag&, + vtkm::IdComponent numPoints = 0) +{ + using VtkcCellShapeTag = typename CellShapeTagVtkmToVtkc::Type; + static_cast(numPoints); // unused + return VtkcCellShapeTag{}; +} + +VTKM_EXEC_CONT +inline vtkc::Polygon make_VtkcCellShapeTag(const vtkm::CellShapeTagPolygon&, + vtkm::IdComponent numPoints = 0) +{ + return vtkc::Polygon(numPoints); +} + +VTKM_EXEC_CONT +inline vtkc::Cell make_VtkcCellShapeTag(const vtkm::CellShapeTagGeneric& tag, + vtkm::IdComponent numPoints = 0) +{ + return vtkc::Cell(static_cast(tag.Id), numPoints); +} + +} // namespace internal + #define vtkmGenericCellShapeMacroCase(cellShapeId, call) \ case vtkm::cellShapeId: \ { \ diff --git a/vtkm/cont/CMakeLists.txt b/vtkm/cont/CMakeLists.txt index 35b2a2deb..2c9ea9068 100644 --- a/vtkm/cont/CMakeLists.txt +++ b/vtkm/cont/CMakeLists.txt @@ -211,7 +211,7 @@ if(TARGET vtkm::openmp) endif() target_link_libraries(vtkm_cont PUBLIC vtkm_compiler_flags ${backends}) -target_link_libraries(vtkm_cont PUBLIC vtkm_taotuple vtkm_optionparser vtkm_diy) +target_link_libraries(vtkm_cont PUBLIC vtkm_taotuple vtkm_optionparser vtkm_diy vtkm_vtkc) if(TARGET vtkm_loguru) target_link_libraries(vtkm_cont PRIVATE vtkm_loguru) endif() diff --git a/vtkm/exec/CMakeLists.txt b/vtkm/exec/CMakeLists.txt index a8f13b578..a38369d36 100644 --- a/vtkm/exec/CMakeLists.txt +++ b/vtkm/exec/CMakeLists.txt @@ -29,7 +29,6 @@ set(headers ExecutionWholeArray.h FieldNeighborhood.h FunctorBase.h - Jacobian.h ParametricCoordinates.h PointLocator.h PointLocatorUniformGrid.h diff --git a/vtkm/exec/CellDerivative.h b/vtkm/exec/CellDerivative.h index ed37aed2e..43a143388 100644 --- a/vtkm/exec/CellDerivative.h +++ b/vtkm/exec/CellDerivative.h @@ -12,392 +12,19 @@ #include #include -#include -#include #include #include -#include #include #include -#include + +#include namespace vtkm { namespace exec { -// The derivative for a 2D polygon in 3D space is underdetermined since there -// is no information in the direction perpendicular to the polygon. To compute -// derivatives for general polygons, we build a 2D space for the polygon's -// plane and solve the derivative there. - -namespace -{ -#define VTKM_ACCUM_PARAMETRIC_DERIVATIVE_3D(pointIndex, weight0, weight1, weight2) \ - parametricDerivative[0] += field[pointIndex] * (weight0); \ - parametricDerivative[1] += field[pointIndex] * (weight1); \ - parametricDerivative[2] += field[pointIndex] * (weight2) - -// Find the derivative of a field in parametric space. That is, find the -// vector [ds/du, ds/dv, ds/dw]. -template -VTKM_EXEC vtkm::Vec ParametricDerivative( - const FieldVecType& field, - const vtkm::Vec& pcoords, - vtkm::CellShapeTagHexahedron) -{ - using FieldType = typename FieldVecType::ComponentType; - using GradientType = vtkm::Vec; - - GradientType pc(pcoords); - GradientType rc = GradientType(FieldType(1)) - pc; - - GradientType parametricDerivative(FieldType(0)); - VTKM_DERIVATIVE_WEIGHTS_HEXAHEDRON(pc, rc, VTKM_ACCUM_PARAMETRIC_DERIVATIVE_3D); - return parametricDerivative; -} - -template -VTKM_EXEC vtkm::Vec ParametricDerivative( - const FieldVecType& field, - const vtkm::Vec& pcoords, - vtkm::CellShapeTagWedge) -{ - using FieldType = typename FieldVecType::ComponentType; - using GradientType = vtkm::Vec; - - GradientType pc(pcoords); - GradientType rc = GradientType(FieldType(1)) - pc; - - GradientType parametricDerivative(FieldType(0)); - VTKM_DERIVATIVE_WEIGHTS_WEDGE(pc, rc, VTKM_ACCUM_PARAMETRIC_DERIVATIVE_3D); - - return parametricDerivative; -} - -template -VTKM_EXEC vtkm::Vec ParametricDerivative( - const FieldVecType& field, - const vtkm::Vec& pcoords, - vtkm::CellShapeTagPyramid) -{ - using FieldType = typename FieldVecType::ComponentType; - using GradientType = vtkm::Vec; - - GradientType pc(pcoords); - GradientType rc = GradientType(FieldType(1)) - pc; - - GradientType parametricDerivative(FieldType(0)); - VTKM_DERIVATIVE_WEIGHTS_PYRAMID(pc, rc, VTKM_ACCUM_PARAMETRIC_DERIVATIVE_3D); - - return parametricDerivative; -} - -#undef VTKM_ACCUM_PARAMETRIC_DERIVATIVE_3D - -#define VTKM_ACCUM_PARAMETRIC_DERIVATIVE_2D(pointIndex, weight0, weight1) \ - parametricDerivative[0] += field[pointIndex] * (weight0); \ - parametricDerivative[1] += field[pointIndex] * (weight1) - -template -VTKM_EXEC vtkm::Vec ParametricDerivative( - const FieldVecType& field, - const vtkm::Vec& pcoords, - vtkm::CellShapeTagQuad) -{ - using FieldType = typename FieldVecType::ComponentType; - using GradientType = vtkm::Vec; - - GradientType pc(static_cast(pcoords[0]), static_cast(pcoords[1])); - GradientType rc = GradientType(FieldType(1)) - pc; - - GradientType parametricDerivative(FieldType(0)); - VTKM_DERIVATIVE_WEIGHTS_QUAD(pc, rc, VTKM_ACCUM_PARAMETRIC_DERIVATIVE_2D); - - return parametricDerivative; -} - -#if 0 -// This code doesn't work, so I'm bailing on it. Instead, I'm just grabbing a -// triangle and finding the derivative of that. If you can do better, please -// implement it. -template -VTKM_EXEC -vtkm::Vec -ParametricDerivative(const FieldVecType &field, - const vtkm::Vec &pcoords, - vtkm::CellShapeTagPolygon) -{ - using FieldType = typename FieldVecType::ComponentType; - using GradientType = vtkm::Vec; - - const vtkm::IdComponent numPoints = field.GetNumberOfComponents(); - FieldType deltaAngle = static_cast(2*vtkm::Pi()/numPoints); - - GradientType pc(pcoords[0], pcoords[1]); - - GradientType parametricDerivative(0); - for (vtkm::IdComponent pointIndex = 0; pointIndex < numPoints; pointIndex++) - { - FieldType angle = pointIndex*deltaAngle; - vtkm::Vec nodePCoords(0.5f*(vtkm::Cos(angle)+1), - 0.5f*(vtkm::Sin(angle)+1)); - - // This is the vector pointing from the user provided parametric coordinate - // to the node at pointIndex in parametric space. - vtkm::Vec pvec = nodePCoords - pc; - - // The weight (the derivative of the interpolation factor) happens to be - // pvec scaled by the cube root of pvec's magnitude. - FieldType magSqr = vtkm::MagnitudeSquared(pvec); - FieldType invMag = vtkm::RSqrt(magSqr); - FieldType scale = invMag*invMag*invMag; - vtkm::Vec weight = scale*pvec; - - parametricDerivative[0] += field[pointIndex] * weight[0]; - parametricDerivative[1] += field[pointIndex] * weight[1]; - } - - return parametricDerivative; -} -#endif - -#undef VTKM_ACCUM_PARAMETRIC_DERIVATIVE_2D - -} // namespace unnamed - -namespace detail -{ - -template -VTKM_EXEC vtkm::Vec CellDerivativeFor3DCell( - const vtkm::Vec& field, - const WorldCoordType& wCoords, - const vtkm::Vec& pcoords, - CellShapeTag tag) -{ - using GradientType = vtkm::Vec; - - // For reasons that should become apparent in a moment, we actually want - // the transpose of the Jacobian. - vtkm::Matrix jacobianTranspose; - vtkm::exec::JacobianFor3DCell(wCoords, pcoords, jacobianTranspose, tag); - jacobianTranspose = vtkm::MatrixTranspose(jacobianTranspose); - - GradientType parametricDerivative = ParametricDerivative(field, pcoords, tag); - - // If we write out the matrices below, it should become clear that the - // Jacobian transpose times the field derivative in world space equals - // the field derivative in parametric space. - // - // | | | | | | - // | dx/du dy/du dz/du | | ds/dx | | ds/du | - // | | | | | | - // | dx/dv dy/dv dz/dv | | ds/dy | = | ds/dv | - // | | | | | | - // | dx/dw dy/dw dz/dw | | ds/dz | | ds/dw | - // | | | | | | - // - // Now we just need to solve this linear system to find the derivative in - // world space. - - bool valid; // Ignored. - return vtkm::SolveLinearSystem(jacobianTranspose, parametricDerivative, valid); -} - -template -VTKM_EXEC vtkm::Vec, 3> CellDerivativeFor3DCell( - const vtkm::Vec, NumPoints>& field, - const WorldCoordType& wCoords, - const vtkm::Vec& pcoords, - CellShapeTag tag) -{ - //We have been given a vector field so we need to solve for each - //component of the vector. For explanation of the logic used see the - //scalar version of CellDerivativeFor3DCell. - vtkm::Matrix perComponentJacobianTranspose; - vtkm::exec::JacobianFor3DCell(wCoords, pcoords, perComponentJacobianTranspose, tag); - perComponentJacobianTranspose = vtkm::MatrixTranspose(perComponentJacobianTranspose); - - bool valid; // Ignored. - vtkm::Vec, 3> result(vtkm::Vec(0.0f)); - vtkm::Vec perPoint; - for (vtkm::IdComponent i = 0; i < N; ++i) - { - for (vtkm::IdComponent c = 0; c < NumPoints; ++c) - { - perPoint[c] = field[c][i]; - } - vtkm::Vec p = ParametricDerivative(perPoint, pcoords, tag); - vtkm::Vec grad = vtkm::SolveLinearSystem(perComponentJacobianTranspose, p, valid); - - result[0][i] += grad[0]; - result[1][i] += grad[1]; - result[2][i] += grad[2]; - } - return result; -} - -template -VTKM_EXEC vtkm::Vec CellDerivativeFor2DCellFinish( - const vtkm::Vec& field, - const vtkm::Matrix& LUFactorization, - const vtkm::IdComponent2& LUPermutation, - const vtkm::Vec& pcoords, - const vtkm::exec::internal::Space2D& space, - CellShapeTag, - vtkm::TypeTraitsScalarTag) -{ - // Finish solving linear equation. See CellDerivativeFor2DCell implementation - // for more detail. - vtkm::Vec parametricDerivative = - ParametricDerivative(field, pcoords, CellShapeTag()); - - vtkm::Vec gradient2D = - vtkm::detail::MatrixLUPSolve(LUFactorization, LUPermutation, parametricDerivative); - - return space.ConvertVecFromSpace(gradient2D); -} - -template -VTKM_EXEC vtkm::Vec CellDerivativeFor2DCellFinish( - const vtkm::Vec& field, - const vtkm::Matrix& LUFactorization, - const vtkm::IdComponent2& LUPermutation, - const vtkm::Vec& pcoords, - const vtkm::exec::internal::Space2D& space, - CellShapeTag, - vtkm::TypeTraitsVectorTag) -{ - using FieldTraits = vtkm::VecTraits; - using FieldComponentType = typename FieldTraits::ComponentType; - - vtkm::Vec gradient(vtkm::TypeTraits::ZeroInitialization()); - - for (vtkm::IdComponent fieldComponent = 0; - fieldComponent < FieldTraits::GetNumberOfComponents(field[0]); - fieldComponent++) - { - vtkm::Vec subField(FieldTraits::GetComponent(field[0], fieldComponent), - FieldTraits::GetComponent(field[1], fieldComponent), - FieldTraits::GetComponent(field[2], fieldComponent), - FieldTraits::GetComponent(field[3], fieldComponent)); - - vtkm::Vec subGradient = CellDerivativeFor2DCellFinish( - subField, - LUFactorization, - LUPermutation, - pcoords, - space, - CellShapeTag(), - typename vtkm::TypeTraits::DimensionalityTag()); - - FieldTraits::SetComponent(gradient[0], fieldComponent, subGradient[0]); - FieldTraits::SetComponent(gradient[1], fieldComponent, subGradient[1]); - FieldTraits::SetComponent(gradient[2], fieldComponent, subGradient[2]); - } - - return gradient; -} - -template -VTKM_EXEC vtkm::Vec CellDerivativeFor2DCellFinish( - const vtkm::Vec& field, - const vtkm::Matrix& LUFactorization, - const vtkm::IdComponent2& LUPermutation, - const vtkm::Vec& pcoords, - const vtkm::exec::internal::Space2D& space, - CellShapeTag, - vtkm::TypeTraitsMatrixTag) -{ - return CellDerivativeFor2DCellFinish(field, - LUFactorization, - LUPermutation, - pcoords, - space, - CellShapeTag(), - vtkm::TypeTraitsVectorTag()); -} - -template -VTKM_EXEC vtkm::Vec CellDerivativeFor2DCell( - const FieldVecType& field, - const WorldCoordType& wCoords, - const vtkm::Vec& pcoords, - CellShapeTag) -{ - using FieldType = typename FieldVecType::ComponentType; - using BaseFieldType = typename vtkm::VecTraits::BaseComponentType; - - // We have an underdetermined system in 3D, so create a 2D space in the - // plane that the polygon sits. - vtkm::exec::internal::Space2D space( - wCoords[0], wCoords[1], wCoords[wCoords.GetNumberOfComponents() - 1]); - - // For reasons that should become apparent in a moment, we actually want - // the transpose of the Jacobian. - vtkm::Matrix jacobianTranspose; - vtkm::exec::JacobianFor2DCell(wCoords, pcoords, space, jacobianTranspose, CellShapeTag()); - jacobianTranspose = vtkm::MatrixTranspose(jacobianTranspose); - - // Find the derivative of the field in parametric coordinate space. That is, - // find the vector [ds/du, ds/dv]. - // Commented because this is actually done in CellDerivativeFor2DCellFinish - // to handle vector fields. - // vtkm::Vec parametricDerivative = - // ParametricDerivative(field,pcoords,CellShapeTag()); - - // If we write out the matrices below, it should become clear that the - // Jacobian transpose times the field derivative in world space equals - // the field derivative in parametric space. - // - // | | | | | | - // | db0/du db1/du | | ds/db0 | | ds/du | - // | | | | = | | - // | db0/dv db1/dv | | ds/db1 | | ds/dv | - // | | | | | | - // - // Now we just need to solve this linear system to find the derivative in - // world space. - - bool valid; // Ignored. - // If you look at the implementation of vtkm::SolveLinearSystem, you will see - // that it is done in two parts. First, it does an LU factorization and then - // uses that result to complete the solve. The factorization part talkes the - // longest amount of time, and if we are performing the gradient on a vector - // field, the factorization can be reused for each component of the vector - // field. Thus, we are going to call the internals of SolveLinearSystem - // ourselves to do the factorization and then apply it to all components. - vtkm::IdComponent2 permutation; - BaseFieldType inversionParity; // Unused - vtkm::detail::MatrixLUPFactor(jacobianTranspose, permutation, inversionParity, valid); - // MatrixLUPFactor does in place factorization. jacobianTranspose is now the - // LU factorization. - return CellDerivativeFor2DCellFinish(field, - jacobianTranspose, - permutation, - pcoords, - space, - CellShapeTag(), - typename vtkm::TypeTraits::DimensionalityTag()); -} - -} // namespace detail - //----------------------------------------------------------------------------- /// \brief Take the derivative (get the gradient) of a point field in a cell. /// @@ -428,6 +55,60 @@ VTKM_EXEC vtkm::Vec CellDerivative( } //----------------------------------------------------------------------------- +namespace internal +{ + +template +VTKM_EXEC vtkm::Vec CellDerivativeImpl( + VtkcCellShapeTag tag, + const FieldVecType& field, + const WorldCoordType& wCoords, + const ParametricCoordType& pcoords, + const vtkm::exec::FunctorBase& worklet) +{ + VTKM_ASSERT(field.GetNumberOfComponents() == tag.numberOfPoints()); + VTKM_ASSERT(wCoords.GetNumberOfComponents() == tag.numberOfPoints()); + + using FieldType = typename FieldVecType::ComponentType; + + auto fieldNumComponents = vtkm::VecTraits::GetNumberOfComponents(field[0]); + vtkm::Vec derivs; + auto status = vtkc::derivative(tag, + vtkc::makeFieldAccessorNestedSOA(wCoords, 3), + vtkc::makeFieldAccessorNestedSOA(field, fieldNumComponents), + pcoords, + derivs[0], + derivs[1], + derivs[2]); + if (status != vtkc::ErrorCode::SUCCESS) + { + worklet.RaiseError(vtkc::errorString(status)); + derivs = vtkm::TypeTraits>::ZeroInitialization(); + } + + return derivs; +} + +} // namespace internal + +template +VTKM_EXEC vtkm::Vec CellDerivative( + const FieldVecType& field, + const WorldCoordType& wCoords, + const vtkm::Vec& pcoords, + CellShapeTag shape, + const vtkm::exec::FunctorBase& worklet) +{ + return internal::CellDerivativeImpl( + vtkm::internal::make_VtkcCellShapeTag(shape), field, wCoords, pcoords, worklet); +} + template VTKM_EXEC vtkm::Vec CellDerivative( const FieldVecType&, @@ -440,116 +121,6 @@ VTKM_EXEC vtkm::Vec CellDerivative( return vtkm::Vec(); } -//----------------------------------------------------------------------------- -template -VTKM_EXEC vtkm::Vec CellDerivative( - const FieldVecType& field, - const WorldCoordType& wCoords, - const vtkm::Vec&, - vtkm::CellShapeTagVertex, - const vtkm::exec::FunctorBase& vtkmNotUsed(worklet)) -{ - (void)field; - (void)wCoords; - VTKM_ASSERT(field.GetNumberOfComponents() == 1); - VTKM_ASSERT(wCoords.GetNumberOfComponents() == 1); - - using GradientType = vtkm::Vec; - return vtkm::TypeTraits::ZeroInitialization(); -} - -//----------------------------------------------------------------------------- -namespace detail -{ - -template -VTKM_EXEC vtkm::Vec CellDerivativeLineImpl( - const FieldType& deltaField, - const WorldCoordType& vec, // direction of line - const typename WorldCoordType::ComponentType& vecMagSqr, - vtkm::TypeTraitsScalarTag) -{ - using GradientType = vtkm::Vec; - - // The derivative of a line is in the direction of the line. Its length is - // equal to the difference of the scalar field divided by the length of the - // line segment. Thus, the derivative is characterized by - // (deltaField*vec)/mag(vec)^2. - - return (deltaField / static_cast(vecMagSqr)) * GradientType(vec); -} - -template -VTKM_EXEC vtkm::Vec CellDerivativeLineImpl( - const FieldType& deltaField, - const WorldCoordType& vec, // direction of line - const typename WorldCoordType::ComponentType& vecMagSqr, - VectorTypeTraitsTag) -{ - using FieldTraits = vtkm::VecTraits; - using FieldComponentType = typename FieldTraits::ComponentType; - using GradientType = vtkm::Vec; - - GradientType gradient(vtkm::TypeTraits::ZeroInitialization()); - for (vtkm::IdComponent fieldComponent = 0; - fieldComponent < FieldTraits::GetNumberOfComponents(deltaField); - fieldComponent++) - { - using SubGradientType = vtkm::Vec; - SubGradientType subGradient = - CellDerivativeLineImpl(FieldTraits::GetComponent(deltaField, fieldComponent), - vec, - vecMagSqr, - typename vtkm::TypeTraits::DimensionalityTag()); - FieldTraits::SetComponent(gradient[0], fieldComponent, subGradient[0]); - FieldTraits::SetComponent(gradient[1], fieldComponent, subGradient[1]); - FieldTraits::SetComponent(gradient[2], fieldComponent, subGradient[2]); - } - - return gradient; -} - -} // namespace detail - -template -VTKM_EXEC vtkm::Vec CellDerivative( - const FieldVecType& field, - const WorldCoordType& wCoords, - const vtkm::Vec& vtkmNotUsed(pcoords), - vtkm::CellShapeTagLine, - const vtkm::exec::FunctorBase& vtkmNotUsed(worklet)) -{ - VTKM_ASSERT(field.GetNumberOfComponents() == 2); - VTKM_ASSERT(wCoords.GetNumberOfComponents() == 2); - - using FieldType = typename FieldVecType::ComponentType; - using BaseComponentType = typename vtkm::VecTraits::BaseComponentType; - - FieldType deltaField(field[1] - field[0]); - vtkm::Vec vec(wCoords[1] - wCoords[0]); - - 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& vtkmNotUsed(pcoords), - vtkm::CellShapeTagLine, - const vtkm::exec::FunctorBase& vtkmNotUsed(worklet)) -{ - VTKM_ASSERT(field.GetNumberOfComponents() == 2); - - using T = typename FieldVecType::ComponentType; - - return vtkm::Vec( - static_cast((field[1] - field[0]) / wCoords.GetSpacing()[0]), T(0), T(0)); -} - template VTKM_EXEC vtkm::Vec CellDerivative( const FieldVecType& field, @@ -570,286 +141,23 @@ VTKM_EXEC vtkm::Vec CellDerivative( return CellDerivative(field, wCoords, pcoords, vtkm::CellShapeTagLine(), worklet); } - using FieldType = typename FieldVecType::ComponentType; - using BaseComponentType = typename vtkm::VecTraits::BaseComponentType; - - ParametricCoordType dt; - dt = static_cast(1) / static_cast(numPoints - 1); - vtkm::IdComponent idx = static_cast(vtkm::Ceil(pcoords[0] / dt)); + auto dt = static_cast(1) / static_cast(numPoints - 1); + auto 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)); + auto lineField = vtkm::make_Vec(field[idx - 1], field[idx]); + auto lineWCoords = vtkm::make_Vec(wCoords[idx - 1], wCoords[idx]); + auto pc = (pcoords[0] - static_cast(idx) * dt) / dt; + return internal::CellDerivativeImpl(vtkc::Line{}, lineField, lineWCoords, &pc, worklet); } -//----------------------------------------------------------------------------- -namespace detail -{ - -template -VTKM_EXEC vtkm::Vec TriangleDerivativeFinish( - const vtkm::Vec& field, - const vtkm::Matrix& LUFactorization, - const vtkm::IdComponent3& LUPermutation, - vtkm::TypeTraitsScalarTag) -{ - // Finish solving linear equation. See TriangleDerivative implementation - // for more detail. - vtkm::Vec b(field[1] - field[0], field[2] - field[0], 0); - - return vtkm::detail::MatrixLUPSolve(LUFactorization, LUPermutation, b); -} - -template -VTKM_EXEC vtkm::Vec TriangleDerivativeFinish( - const vtkm::Vec& field, - const vtkm::Matrix& LUFactorization, - const vtkm::IdComponent3& LUPermutation, - vtkm::TypeTraitsVectorTag) -{ - using FieldTraits = vtkm::VecTraits; - using FieldComponentType = typename FieldTraits::ComponentType; - - vtkm::Vec gradient(vtkm::TypeTraits::ZeroInitialization()); - - for (vtkm::IdComponent fieldComponent = 0; - fieldComponent < FieldTraits::GetNumberOfComponents(field[0]); - fieldComponent++) - { - vtkm::Vec subField(FieldTraits::GetComponent(field[0], fieldComponent), - FieldTraits::GetComponent(field[1], fieldComponent), - FieldTraits::GetComponent(field[2], fieldComponent)); - - vtkm::Vec subGradient = - TriangleDerivativeFinish(subField, - LUFactorization, - LUPermutation, - typename vtkm::TypeTraits::DimensionalityTag()); - - FieldTraits::SetComponent(gradient[0], fieldComponent, subGradient[0]); - FieldTraits::SetComponent(gradient[1], fieldComponent, subGradient[1]); - FieldTraits::SetComponent(gradient[2], fieldComponent, subGradient[2]); - } - - return gradient; -} - -template -VTKM_EXEC vtkm::Vec TriangleDerivativeFinish( - const vtkm::Vec& field, - const vtkm::Matrix& LUFactorization, - const vtkm::IdComponent3& LUPermutation, - vtkm::TypeTraitsMatrixTag) -{ - return TriangleDerivativeFinish( - field, LUFactorization, LUPermutation, vtkm::TypeTraitsVectorTag()); -} - -template -VTKM_EXEC vtkm::Vec TriangleDerivative(const vtkm::Vec& field, - const vtkm::Vec& wCoords) -{ - using BaseComponentType = typename vtkm::VecTraits::BaseComponentType; - - // The scalar values of the three points in a triangle completely specify a - // linear field (with constant gradient) assuming the field is constant in - // the normal direction to the triangle. The field, defined by the 3-vector - // gradient g and scalar value s_origin, can be found with this set of 4 - // equations and 4 unknowns. - // - // Dot(p0, g) + s_origin = s0 - // Dot(p1, g) + s_origin = s1 - // Dot(p2, g) + s_origin = s2 - // Dot(n, g) = 0 - // - // Where the p's are point coordinates and n is the normal vector. But we - // don't really care about s_origin. We just want to find the gradient g. - // With some simple elimination we, we can get rid of s_origin and be left - // with 3 equations and 3 unknowns. - // - // Dot(p1-p0, g) = s1 - s0 - // Dot(p2-p0, g) = s2 - s0 - // Dot(n, g) = 0 - // - // We'll solve this by putting this in matrix form Ax = b where the rows of A - // are the differences in points and normal, b has the scalar differences, - // and x is really the gradient g. - - vtkm::Vec v0 = wCoords[1] - wCoords[0]; - vtkm::Vec v1 = wCoords[2] - wCoords[0]; - vtkm::Vec n = vtkm::Cross(v0, v1); - - vtkm::Matrix A; - vtkm::MatrixSetRow(A, 0, v0); - vtkm::MatrixSetRow(A, 1, v1); - vtkm::MatrixSetRow(A, 2, n); - - // If the triangle is degenerate, then valid will be false. For now we are - // ignoring it. We could detect it if we determine we need to although I have - // seen singular matrices missed due to floating point error. - // - bool valid; - - // If you look at the implementation of vtkm::SolveLinearSystem, you will see - // that it is done in two parts. First, it does an LU factorization and then - // uses that result to complete the solve. The factorization part talkes the - // longest amount of time, and if we are performing the gradient on a vector - // field, the factorization can be reused for each component of the vector - // field. Thus, we are going to call the internals of SolveLinearSystem - // ourselves to do the factorization and then apply it to all components. - vtkm::IdComponent3 permutation; - BaseComponentType inversionParity; // Unused - vtkm::detail::MatrixLUPFactor(A, permutation, inversionParity, valid); - // MatrixLUPFactor does in place factorization. A is now the LU factorization. - return TriangleDerivativeFinish( - field, A, permutation, typename vtkm::TypeTraits::DimensionalityTag()); -} - -} // namespace detail - -//----------------------------------------------------------------------------- -template -VTKM_EXEC vtkm::Vec CellDerivative( - const FieldVecType& inputField, - const WorldCoordType& wCoords, - const vtkm::Vec& vtkmNotUsed(pcoords), - vtkm::CellShapeTagTriangle, - const vtkm::exec::FunctorBase& vtkmNotUsed(worklet)) -{ - VTKM_ASSERT(inputField.GetNumberOfComponents() == 3); - VTKM_ASSERT(wCoords.GetNumberOfComponents() == 3); - - using ValueType = typename FieldVecType::ComponentType; - using WCoordType = typename WorldCoordType::ComponentType; - vtkm::Vec field; - inputField.CopyInto(field); - vtkm::Vec wpoints; - wCoords.CopyInto(wpoints); - return detail::TriangleDerivative(field, wpoints); -} - -namespace detail -{ - -//----------------------------------------------------------------------------- -template -VTKM_EXEC void PolygonSubTriangle(const FieldVecType& inField, - const WorldCoordType& inWCoords, - const vtkm::Vec& pcoords, - vtkm::Vec& outField, - vtkm::Vec& outWCoords, - const vtkm::exec::FunctorBase& worklet) -{ - // To find the gradient in a polygon (of 5 or more points), we will extract a small triangle near - // the desired parameteric coordinates (pcoords). We return the field values (outField) and world - // coordinates (outWCoords) for this triangle, which is all that is needed to find the gradient - // in a triangle. - // - // The trangle will be "pointing" away from the center of the polygon, and pcoords will be placed - // at the apex of the triangle. This is because if pcoords is at or near the edge of the polygon, - // we do not want to push any of the points over the edge, and it is not trivial to determine - // exactly where the edge of the polygon is. - - // First point is right at pcoords - outField[0] = vtkm::exec::CellInterpolate(inField, pcoords, vtkm::CellShapeTagPolygon{}, worklet); - outWCoords[0] = - vtkm::exec::CellInterpolate(inWCoords, pcoords, vtkm::CellShapeTagPolygon{}, worklet); - - // Find the unit vector pointing from the center of the polygon to pcoords - vtkm::Vec radialVector = { pcoords[0] - 0.5f, pcoords[1] - 0.5f }; - ParametricCoordType magnitudeSquared = vtkm::MagnitudeSquared(radialVector); - if (magnitudeSquared > 8 * vtkm::Epsilon()) - { - radialVector = vtkm::RSqrt(magnitudeSquared) * radialVector; - } - else - { - // pcoords is in the center of the polygon. Just point in an arbitrary direction - radialVector = { 1.0, 0.0 }; - } - - // We want the two points away from pcoords to be back toward the center but moved at 45 degrees - // off the radius. Simple geometry shows us that the (not quite unit) vectors of those two - // directions are (-radialVector[1] - radialVector[0], radialVector[0] - radialVector[1]) and - // (radialVector[1] - radialVector[0], -radialVector[0] - radialVector[1]). - // - // *\ (-radialVector[1], radialVector[0]) // - // | \ // - // | \ (-radialVector[1] - radialVector[0], radialVector[0] - radialVector[1]) // - // | \ // - // +-------* radialVector // - // | / // - // | / (radialVector[1] - radialVector[0], -radialVector[0] - radialVector[1]) // - // | / // - // */ (radialVector[1], -radialVector[0]) // - - // This scaling value is somewhat arbitrary. It is small enough to be "close" to the selected - // point and small enough to be guaranteed to be inside the polygon, but large enough to to - // get an accurate gradient. - static constexpr ParametricCoordType scale = 0.05f; - - vtkm::Vec backPcoord = { - pcoords[0] + scale * (-radialVector[1] - radialVector[0]), - pcoords[1] + scale * (radialVector[0] - radialVector[1]), - 0 - }; - outField[1] = - vtkm::exec::CellInterpolate(inField, backPcoord, vtkm::CellShapeTagPolygon{}, worklet); - outWCoords[1] = - vtkm::exec::CellInterpolate(inWCoords, backPcoord, vtkm::CellShapeTagPolygon{}, worklet); - - backPcoord = { pcoords[0] + scale * (radialVector[1] - radialVector[0]), - pcoords[1] + scale * (-radialVector[0] - radialVector[1]), - 0 }; - outField[2] = - vtkm::exec::CellInterpolate(inField, backPcoord, vtkm::CellShapeTagPolygon{}, worklet); - outWCoords[2] = - vtkm::exec::CellInterpolate(inWCoords, backPcoord, vtkm::CellShapeTagPolygon{}, worklet); -} - -} // namespace detail - //----------------------------------------------------------------------------- template VTKM_EXEC vtkm::Vec CellDerivative( @@ -870,232 +178,22 @@ VTKM_EXEC vtkm::Vec CellDerivative( return CellDerivative(field, wCoords, pcoords, vtkm::CellShapeTagVertex(), worklet); case 2: return CellDerivative(field, wCoords, pcoords, vtkm::CellShapeTagLine(), worklet); - case 3: - return CellDerivative(field, wCoords, pcoords, vtkm::CellShapeTagTriangle(), worklet); - case 4: - return CellDerivative(field, wCoords, pcoords, vtkm::CellShapeTagQuad(), worklet); + default: + return internal::CellDerivativeImpl( + vtkc::Polygon(numPoints), field, wCoords, pcoords, worklet); } - - vtkm::Vec triangleField; - vtkm::Vec triangleWCoords; - detail::PolygonSubTriangle(field, wCoords, pcoords, triangleField, triangleWCoords, worklet); - return detail::TriangleDerivative(triangleField, triangleWCoords); } //----------------------------------------------------------------------------- -template -VTKM_EXEC vtkm::Vec CellDerivative( - const FieldVecType& inputField, - const WorldCoordType& wCoords, - const vtkm::Vec& pcoords, - vtkm::CellShapeTagQuad, - const vtkm::exec::FunctorBase& vtkmNotUsed(worklet)) -{ - VTKM_ASSERT(inputField.GetNumberOfComponents() == 4); - VTKM_ASSERT(wCoords.GetNumberOfComponents() == 4); - - using ValueType = typename FieldVecType::ComponentType; - using WCoordType = typename WorldCoordType::ComponentType; - vtkm::Vec field; - inputField.CopyInto(field); - vtkm::Vec wpoints; - wCoords.CopyInto(wpoints); - - return detail::CellDerivativeFor2DCell(field, wpoints, pcoords, vtkm::CellShapeTagQuad()); -} - template VTKM_EXEC vtkm::Vec CellDerivative( const FieldVecType& field, const vtkm::VecAxisAlignedPointCoordinates<2>& wCoords, const vtkm::Vec& pcoords, vtkm::CellShapeTagQuad, - const vtkm::exec::FunctorBase& vtkmNotUsed(worklet)) + const vtkm::exec::FunctorBase& worklet) { - VTKM_ASSERT(field.GetNumberOfComponents() == 4); - - using T = typename FieldVecType::ComponentType; - using VecT = vtkm::Vec; - - VecT pc(static_cast::ComponentType>(pcoords[0]), - static_cast::ComponentType>(pcoords[1])); - VecT rc = VecT(T(1)) - pc; - - VecT sum = field[0] * VecT(-rc[1], -rc[0]); - sum = sum + field[1] * VecT(rc[1], -pc[0]); - sum = sum + field[2] * VecT(pc[1], pc[0]); - sum = sum + field[3] * VecT(-pc[1], rc[0]); - - return vtkm::Vec(static_cast(sum[0] / wCoords.GetSpacing()[0]), - static_cast(sum[1] / wCoords.GetSpacing()[1]), - T(0)); -} - -//----------------------------------------------------------------------------- -namespace detail -{ - -template -VTKM_EXEC vtkm::Vec TetraDerivativeFinish( - const vtkm::Vec& field, - const vtkm::Matrix& LUFactorization, - const vtkm::IdComponent3& LUPermutation, - vtkm::TypeTraitsScalarTag) -{ - // Finish solving linear equation. See TriangleDerivative implementation - // for more detail. - vtkm::Vec b(field[1] - field[0], field[2] - field[0], field[3] - field[0]); - - return vtkm::detail::MatrixLUPSolve(LUFactorization, LUPermutation, b); -} - -template -VTKM_EXEC vtkm::Vec TetraDerivativeFinish( - const vtkm::Vec& field, - const vtkm::Matrix& LUFactorization, - const vtkm::IdComponent3& LUPermutation, - vtkm::TypeTraitsVectorTag) -{ - using FieldTraits = vtkm::VecTraits; - using FieldComponentType = typename FieldTraits::ComponentType; - - vtkm::Vec gradient(vtkm::TypeTraits::ZeroInitialization()); - - for (vtkm::IdComponent fieldComponent = 0; - fieldComponent < FieldTraits::GetNumberOfComponents(field[0]); - fieldComponent++) - { - vtkm::Vec subField(FieldTraits::GetComponent(field[0], fieldComponent), - FieldTraits::GetComponent(field[1], fieldComponent), - FieldTraits::GetComponent(field[2], fieldComponent), - FieldTraits::GetComponent(field[3], fieldComponent)); - - vtkm::Vec subGradient = - TetraDerivativeFinish(subField, - LUFactorization, - LUPermutation, - typename vtkm::TypeTraits::DimensionalityTag()); - - FieldTraits::SetComponent(gradient[0], fieldComponent, subGradient[0]); - FieldTraits::SetComponent(gradient[1], fieldComponent, subGradient[1]); - FieldTraits::SetComponent(gradient[2], fieldComponent, subGradient[2]); - } - - return gradient; -} - -template -VTKM_EXEC vtkm::Vec TetraDerivativeFinish( - const vtkm::Vec& field, - const vtkm::Matrix& LUFactorization, - const vtkm::IdComponent3& LUPermutation, - vtkm::TypeTraitsMatrixTag) -{ - return TetraDerivativeFinish(field, LUFactorization, LUPermutation, vtkm::TypeTraitsVectorTag()); -} - -template -VTKM_EXEC vtkm::Vec TetraDerivative(const vtkm::Vec& field, - const vtkm::Vec& wCoords) -{ - using BaseComponentType = typename vtkm::VecTraits::BaseComponentType; - - // The scalar values of the four points in a tetrahedron completely specify a - // linear field (with constant gradient). The field, defined by the 3-vector - // gradient g and scalar value s_origin, can be found with this set of 4 - // equations and 4 unknowns. - // - // Dot(p0, g) + s_origin = s0 - // Dot(p1, g) + s_origin = s1 - // Dot(p2, g) + s_origin = s2 - // Dot(p3, g) + s_origin = s3 - // - // Where the p's are point coordinates. But we don't really care about - // s_origin. We just want to find the gradient g. With some simple - // elimination we, we can get rid of s_origin and be left with 3 equations - // and 3 unknowns. - // - // Dot(p1-p0, g) = s1 - s0 - // Dot(p2-p0, g) = s2 - s0 - // Dot(p3-p0, g) = s3 - s0 - // - // We'll solve this by putting this in matrix form Ax = b where the rows of A - // are the differences in points and normal, b has the scalar differences, - // and x is really the gradient g. - - vtkm::Vec v0 = wCoords[1] - wCoords[0]; - vtkm::Vec v1 = wCoords[2] - wCoords[0]; - vtkm::Vec v2 = wCoords[3] - wCoords[0]; - - vtkm::Matrix A; - vtkm::MatrixSetRow(A, 0, v0); - vtkm::MatrixSetRow(A, 1, v1); - vtkm::MatrixSetRow(A, 2, v2); - - // If the tetrahedron is degenerate, then valid will be false. For now we are - // ignoring it. We could detect it if we determine we need to although I have - // seen singular matrices missed due to floating point error. - // - bool valid; - - // If you look at the implementation of vtkm::SolveLinearSystem, you will see - // that it is done in two parts. First, it does an LU factorization and then - // uses that result to complete the solve. The factorization part talkes the - // longest amount of time, and if we are performing the gradient on a vector - // field, the factorization can be reused for each component of the vector - // field. Thus, we are going to call the internals of SolveLinearSystem - // ourselves to do the factorization and then apply it to all components. - vtkm::IdComponent3 permutation; - BaseComponentType inversionParity; // Unused - vtkm::detail::MatrixLUPFactor(A, permutation, inversionParity, valid); - // MatrixLUPFactor does in place factorization. A is now the LU factorization. - return TetraDerivativeFinish( - field, A, permutation, typename vtkm::TypeTraits::DimensionalityTag()); -} - -} // namespace detail - -//----------------------------------------------------------------------------- -template -VTKM_EXEC vtkm::Vec CellDerivative( - const FieldVecType& inputField, - const WorldCoordType& wCoords, - const vtkm::Vec& vtkmNotUsed(pcoords), - vtkm::CellShapeTagTetra, - const vtkm::exec::FunctorBase& vtkmNotUsed(worklet)) -{ - VTKM_ASSERT(inputField.GetNumberOfComponents() == 4); - VTKM_ASSERT(wCoords.GetNumberOfComponents() == 4); - - using ValueType = typename FieldVecType::ComponentType; - using WCoordType = typename WorldCoordType::ComponentType; - vtkm::Vec field; - inputField.CopyInto(field); - vtkm::Vec wpoints; - wCoords.CopyInto(wpoints); - return detail::TetraDerivative(field, wpoints); -} - -//----------------------------------------------------------------------------- -template -VTKM_EXEC vtkm::Vec CellDerivative( - const FieldVecType& inputField, - const WorldCoordType& wCoords, - const vtkm::Vec& pcoords, - vtkm::CellShapeTagHexahedron, - const vtkm::exec::FunctorBase& vtkmNotUsed(worklet)) -{ - VTKM_ASSERT(inputField.GetNumberOfComponents() == 8); - VTKM_ASSERT(wCoords.GetNumberOfComponents() == 8); - - using ValueType = typename FieldVecType::ComponentType; - using WCoordType = typename WorldCoordType::ComponentType; - vtkm::Vec field; - inputField.CopyInto(field); - vtkm::Vec wpoints; - wCoords.CopyInto(wpoints); - - return detail::CellDerivativeFor3DCell(field, wpoints, pcoords, vtkm::CellShapeTagHexahedron()); + return internal::CellDerivativeImpl(vtkc::Pixel{}, field, wCoords, pcoords, worklet); } template @@ -1104,93 +202,9 @@ VTKM_EXEC vtkm::Vec CellDerivative( const vtkm::VecAxisAlignedPointCoordinates<3>& wCoords, const vtkm::Vec& pcoords, vtkm::CellShapeTagHexahedron, - const vtkm::exec::FunctorBase& vtkmNotUsed(worklet)) + const vtkm::exec::FunctorBase& worklet) { - VTKM_ASSERT(field.GetNumberOfComponents() == 8); - - using T = typename FieldVecType::ComponentType; - using VecT = vtkm::Vec; - - VecT pc(static_cast::ComponentType>(pcoords[0]), - static_cast::ComponentType>(pcoords[1]), - static_cast::ComponentType>(pcoords[2])); - VecT rc = VecT(T(1)) - pc; - - VecT sum = field[0] * VecT(-rc[1] * rc[2], -rc[0] * rc[2], -rc[0] * rc[1]); - sum = sum + field[1] * VecT(rc[1] * rc[2], -pc[0] * rc[2], -pc[0] * rc[1]); - sum = sum + field[2] * VecT(pc[1] * rc[2], pc[0] * rc[2], -pc[0] * pc[1]); - sum = sum + field[3] * VecT(-pc[1] * rc[2], rc[0] * rc[2], -rc[0] * pc[1]); - sum = sum + field[4] * VecT(-rc[1] * pc[2], -rc[0] * pc[2], rc[0] * rc[1]); - sum = sum + field[5] * VecT(rc[1] * pc[2], -pc[0] * pc[2], pc[0] * rc[1]); - sum = sum + field[6] * VecT(pc[1] * pc[2], pc[0] * pc[2], pc[0] * pc[1]); - sum = sum + field[7] * VecT(-pc[1] * pc[2], rc[0] * pc[2], rc[0] * pc[1]); - - return sum / wCoords.GetSpacing(); -} - -//----------------------------------------------------------------------------- -template -VTKM_EXEC vtkm::Vec CellDerivative( - const FieldVecType& inputField, - const WorldCoordType& wCoords, - const vtkm::Vec& pcoords, - vtkm::CellShapeTagWedge, - const vtkm::exec::FunctorBase& vtkmNotUsed(worklet)) -{ - VTKM_ASSERT(inputField.GetNumberOfComponents() == 6); - VTKM_ASSERT(wCoords.GetNumberOfComponents() == 6); - - using ValueType = typename FieldVecType::ComponentType; - using WCoordType = typename WorldCoordType::ComponentType; - vtkm::Vec field; - inputField.CopyInto(field); - vtkm::Vec wpoints; - wCoords.CopyInto(wpoints); - - return detail::CellDerivativeFor3DCell(field, wpoints, pcoords, vtkm::CellShapeTagWedge()); -} - -//----------------------------------------------------------------------------- -template -VTKM_EXEC vtkm::Vec CellDerivative( - const FieldVecType& inputField, - const WorldCoordType& wCoords, - const vtkm::Vec& pcoords, - vtkm::CellShapeTagPyramid, - const vtkm::exec::FunctorBase& vtkmNotUsed(worklet)) -{ - VTKM_ASSERT(inputField.GetNumberOfComponents() == 5); - VTKM_ASSERT(wCoords.GetNumberOfComponents() == 5); - - using ValueType = typename FieldVecType::ComponentType; - using WCoordType = typename WorldCoordType::ComponentType; - vtkm::Vec field; - inputField.CopyInto(field); - vtkm::Vec wpoints; - wCoords.CopyInto(wpoints); - - if (pcoords[2] > static_cast(0.999)) - { - // If we are at the apex of the pyramid we need to do something special. - // As we approach the apex, the derivatives of the parametric shape - // functions in x and y go to 0 while the inverse of the Jacobian - // also goes to 0. This results in 0/0 but using l'Hopital's rule - // we could actually compute the value of the limit, if we had a - // functional expression to compute the gradient. We're on a computer - // so we don't but we can cheat and do a linear extrapolation of the - // derivatives which really ends up as the same thing. - using PCoordType = vtkm::Vec; - using ReturnType = vtkm::Vec; - PCoordType pcoords1(0.5f, 0.5f, (2 * 0.998f) - pcoords[2]); - ReturnType derivative1 = - detail::CellDerivativeFor3DCell(field, wpoints, pcoords1, vtkm::CellShapeTagPyramid{}); - PCoordType pcoords2(0.5f, 0.5f, 0.998f); - ReturnType derivative2 = - detail::CellDerivativeFor3DCell(field, wpoints, pcoords2, vtkm::CellShapeTagPyramid{}); - return (ValueType(2) * derivative2) - derivative1; - } - - return detail::CellDerivativeFor3DCell(field, wpoints, pcoords, vtkm::CellShapeTagPyramid()); + return internal::CellDerivativeImpl(vtkc::Voxel{}, field, wCoords, pcoords, worklet); } } } // namespace vtkm::exec diff --git a/vtkm/exec/CellInside.h b/vtkm/exec/CellInside.h index 99bc777de..b3c7a6499 100644 --- a/vtkm/exec/CellInside.h +++ b/vtkm/exec/CellInside.h @@ -13,84 +13,32 @@ #include #include +#include + namespace vtkm { namespace exec { +template +static inline VTKM_EXEC bool CellInside(const vtkm::Vec& pcoords, CellShapeTag) +{ + using VtkcTagType = typename vtkm::internal::CellShapeTagVtkmToVtkc::Type; + return vtkc::cellInside(VtkcTagType{}, pcoords); +} + template static inline VTKM_EXEC bool CellInside(const vtkm::Vec&, vtkm::CellShapeTagEmpty) { return false; } -template -static inline VTKM_EXEC bool CellInside(const vtkm::Vec& pcoords, vtkm::CellShapeTagVertex) -{ - return pcoords[0] == T(0) && pcoords[1] == T(0) && pcoords[2] == T(0); -} - -template -static inline VTKM_EXEC bool CellInside(const vtkm::Vec& pcoords, vtkm::CellShapeTagLine) -{ - 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) -{ - return pcoords[0] >= T(0) && pcoords[1] >= T(0) && (pcoords[0] + pcoords[1] <= T(1)); -} - -template -static inline VTKM_EXEC bool CellInside(const vtkm::Vec& pcoords, vtkm::CellShapeTagPolygon) -{ - return ((pcoords[0] - T(0.5)) * (pcoords[0] - T(0.5))) + - ((pcoords[1] - T(0.5)) * (pcoords[1] - T(0.5))) <= - T(0.25); -} - -template -static inline VTKM_EXEC bool CellInside(const vtkm::Vec& pcoords, vtkm::CellShapeTagQuad) -{ - return pcoords[0] >= T(0) && pcoords[0] <= T(1) && pcoords[1] >= T(0) && pcoords[1] <= T(1); -} - -template -static inline VTKM_EXEC bool CellInside(const vtkm::Vec& pcoords, vtkm::CellShapeTagTetra) -{ - return pcoords[0] >= T(0) && pcoords[1] >= T(0) && pcoords[2] >= T(0) && - (pcoords[0] + pcoords[1] + pcoords[2] <= T(1)); -} - -template -static inline VTKM_EXEC bool CellInside(const vtkm::Vec& pcoords, - vtkm::CellShapeTagHexahedron) -{ - return pcoords[0] >= T(0) && pcoords[0] <= T(1) && pcoords[1] >= T(0) && pcoords[1] <= T(1) && - pcoords[2] >= T(0) && pcoords[2] <= T(1); -} - -template -static inline VTKM_EXEC bool CellInside(const vtkm::Vec& pcoords, vtkm::CellShapeTagWedge) -{ - return pcoords[0] >= T(0) && pcoords[1] >= T(0) && pcoords[2] >= T(0) && pcoords[2] <= T(1) && - (pcoords[0] + pcoords[1] <= T(1)); -} - -template -static inline VTKM_EXEC bool CellInside(const vtkm::Vec& pcoords, vtkm::CellShapeTagPyramid) -{ - return pcoords[0] >= T(0) && pcoords[0] <= T(1) && pcoords[1] >= T(0) && pcoords[1] <= T(1) && - pcoords[2] >= T(0) && pcoords[2] <= T(1); -} - /// Checks if the parametric coordinates `pcoords` are on the inside for the /// specified cell type. /// diff --git a/vtkm/exec/CellInterpolate.h b/vtkm/exec/CellInterpolate.h index 35bfa4bf9..0c66e2c8b 100644 --- a/vtkm/exec/CellInterpolate.h +++ b/vtkm/exec/CellInterpolate.h @@ -12,11 +12,11 @@ #include #include -#include #include -#include #include +#include + #if (defined(VTKM_GCC) || defined(VTKM_CLANG)) #pragma GCC diagnostic push #pragma GCC diagnostic ignored "-Wconversion" @@ -30,97 +30,30 @@ namespace exec namespace internal { -// This is really the WorldCoorindatesToParametericCoordinates function, but -// moved to this header file because it is required to interpolate in a -// polygon, which is divided into triangles. -template -VTKM_EXEC typename WorldCoordVector::ComponentType ReverseInterpolateTriangle( - const WorldCoordVector& pointWCoords, - const typename WorldCoordVector::ComponentType& wcoords) +template +VTKM_EXEC typename FieldVecType::ComponentType CellInterpolateImpl( + VtkcCellShapeTag tag, + const FieldVecType& field, + const ParametricCoordType& pcoords, + const vtkm::exec::FunctorBase& worklet) { - VTKM_ASSERT(pointWCoords.GetNumberOfComponents() == 3); + VTKM_ASSERT(tag.numberOfPoints() == field.GetNumberOfComponents()); - // We will solve the world to parametric coordinates problem geometrically. - // Consider the parallelogram formed by wcoords and p0 of the triangle and - // the two adjacent edges. This parallelogram is equivalent to the - // axis-aligned rectangle anchored at the origin of parametric space. - // - // p2 |\ (1,0) |\ // - // | \ | \ // - // | \ | \ // - // | \ | \ // - // | \ | \ // - // | \ | (u,v) \ // - // | --- \ |-------* \ // - // | ---*wcoords | | \ // - // | | \ | | \ // - // p0 *--- | \ (0,0) *------------------\ (1,0) // - // ---| \ // - // x-- \ // - // --- \ // - // ---\ p1 // - // - // In this diagram, the distance between p0 and the point marked x divided by - // the length of the edge it is on is equal, by proportionality, to the u - // parametric coordinate. (The v coordinate follows the other edge - // accordingly.) Thus, if we can find the intersection at x (or more - // specifically the distance between p0 and x), then we can find that - // parametric coordinate. - // - // Because the triangle is in 3-space, we are actually going to intersect the - // edge with a plane that is parallel to the opposite edge of p0 and - // perpendicular to the triangle. This is partially because it is easy to - // find the intersection between a plane and a line and partially because the - // computation will work for points not on the plane. (The result is - // equivalent to a point projected on the plane.) - // - // First, we define an implicit plane as: - // - // Dot((p - wcoords), planeNormal) = 0 - // - // where planeNormal is the normal to the plane (easily computed from the - // triangle), and p is any point in the plane. Next, we define the parametric - // form of the line: - // - // p(d) = (p1 - p0)d + p0 - // - // Where d is the fraction of distance from p0 toward p1. Note that d is - // actually equal to the parametric coordinate we are trying to find. Once we - // compute it, we are done. We can skip the part about finding the actual - // coordinates of the intersection. - // - // Solving for the intersection is as simple as substituting the line's - // definition of p(d) into p for the plane equation. With some basic algebra - // you get: - // - // d = Dot((wcoords - p0), planeNormal)/Dot((p1-p0), planeNormal) - // - // From here, the u coordinate is simply d. The v coordinate follows - // similarly. - // - - using Vector3 = typename WorldCoordVector::ComponentType; - using T = typename Vector3::ComponentType; - - Vector3 pcoords(T(0)); - Vector3 triangleNormal = vtkm::TriangleNormal(pointWCoords[0], pointWCoords[1], pointWCoords[2]); - - for (vtkm::IdComponent dimension = 0; dimension < 2; dimension++) + using FieldValueType = typename FieldVecType::ComponentType; + IdComponent numComponents = vtkm::VecTraits::GetNumberOfComponents(field[0]); + FieldValueType result(0); + auto status = + vtkc::interpolate(tag, vtkc::makeFieldAccessorNestedSOA(field, numComponents), pcoords, result); + if (status != vtkc::ErrorCode::SUCCESS) { - Vector3 p0 = pointWCoords[0]; - Vector3 p1 = pointWCoords[dimension + 1]; - Vector3 p2 = pointWCoords[2 - dimension]; - Vector3 planeNormal = vtkm::Cross(triangleNormal, p2 - p0); - - T d = vtkm::Dot(wcoords - p0, planeNormal) / vtkm::Dot(p1 - p0, planeNormal); - - pcoords[dimension] = d; + worklet.RaiseError(vtkc::errorString(status)); } - - return pcoords; -} + return result; } +} // namespace internal + +//----------------------------------------------------------------------------- /// \brief Interpolate a point field in a cell. /// /// Given the point field values for each node and the parametric coordinates @@ -145,6 +78,19 @@ VTKM_EXEC typename FieldVecType::ComponentType CellInterpolate( return result; } +//----------------------------------------------------------------------------- +template +VTKM_EXEC typename FieldVecType::ComponentType CellInterpolate( + const FieldVecType& pointFieldValues, + const vtkm::Vec& pcoords, + CellShapeTag tag, + const vtkm::exec::FunctorBase& worklet) +{ + auto vtkcTag = + vtkm::internal::make_VtkcCellShapeTag(tag, pointFieldValues.GetNumberOfComponents()); + return internal::CellInterpolateImpl(vtkcTag, pointFieldValues, pcoords, worklet); +} + //----------------------------------------------------------------------------- template VTKM_EXEC typename FieldVecType::ComponentType CellInterpolate( @@ -157,45 +103,6 @@ VTKM_EXEC typename FieldVecType::ComponentType CellInterpolate( return typename FieldVecType::ComponentType(); } -//----------------------------------------------------------------------------- -template -VTKM_EXEC typename FieldVecType::ComponentType CellInterpolate( - const FieldVecType& pointFieldValues, - const vtkm::Vec, - vtkm::CellShapeTagVertex, - const vtkm::exec::FunctorBase& vtkmNotUsed(worklet)) -{ - VTKM_ASSERT(pointFieldValues.GetNumberOfComponents() == 1); - return pointFieldValues[0]; -} - -//----------------------------------------------------------------------------- -template -VTKM_EXEC typename FieldVecType::ComponentType CellInterpolate( - const FieldVecType& pointFieldValues, - const vtkm::Vec& parametricCoords, - vtkm::CellShapeTagLine, - const vtkm::exec::FunctorBase& vtkmNotUsed(worklet)) -{ - VTKM_ASSERT(pointFieldValues.GetNumberOfComponents() == 2); - return vtkm::Lerp(pointFieldValues[0], pointFieldValues[1], parametricCoords[0]); -} - -template -VTKM_EXEC vtkm::Vec3f CellInterpolate(const vtkm::VecAxisAlignedPointCoordinates<1>& field, - const vtkm::Vec& pcoords, - vtkm::CellShapeTagLine, - const vtkm::exec::FunctorBase&) -{ - using T = vtkm::Vec3f; - - 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( @@ -207,12 +114,9 @@ VTKM_EXEC typename FieldVecType::ComponentType CellInterpolate( const vtkm::IdComponent numPoints = field.GetNumberOfComponents(); VTKM_ASSERT(numPoints >= 1); - switch (numPoints) + if (numPoints == 1) { - case 1: - return CellInterpolate(field, pcoords, vtkm::CellShapeTagVertex(), worklet); - case 2: - return CellInterpolate(field, pcoords, vtkm::CellShapeTagLine(), worklet); + return CellInterpolate(field, pcoords, vtkm::CellShapeTagVertex(), worklet); } using T = ParametricCoordType; @@ -220,50 +124,13 @@ VTKM_EXEC typename FieldVecType::ComponentType CellInterpolate( 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::Vec3f 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); + return field[numPoints - 1]; } - using T = vtkm::Vec3f; - 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( - const FieldVecType& field, - const vtkm::Vec& pcoords, - vtkm::CellShapeTagTriangle, - const vtkm::exec::FunctorBase& vtkmNotUsed(worklet)) -{ - VTKM_ASSERT(field.GetNumberOfComponents() == 3); - using T = typename FieldVecType::ComponentType; - return static_cast((field[0] * (1 - pcoords[0] - pcoords[1])) + (field[1] * pcoords[0]) + - (field[2] * pcoords[1])); + T pc = (pcoords[0] - static_cast(idx) * dt) / dt; + return internal::CellInterpolateImpl( + vtkc::Line{}, vtkm::make_Vec(field[idx], field[idx + 1]), &pc, worklet); } //----------------------------------------------------------------------------- @@ -282,208 +149,29 @@ VTKM_EXEC typename FieldVecType::ComponentType CellInterpolate( return CellInterpolate(field, pcoords, vtkm::CellShapeTagVertex(), worklet); case 2: return CellInterpolate(field, pcoords, vtkm::CellShapeTagLine(), worklet); - case 3: - return CellInterpolate(field, pcoords, vtkm::CellShapeTagTriangle(), worklet); - case 4: - return CellInterpolate(field, pcoords, vtkm::CellShapeTagQuad(), worklet); + default: + return internal::CellInterpolateImpl(vtkc::Polygon(numPoints), field, pcoords, worklet); } - - // If we are here, then there are 5 or more points on this polygon. - - // Arrange the points such that they are on the circle circumscribed in the - // unit square from 0 to 1. That is, the point are on the circle centered at - // coordinate 0.5,0.5 with radius 0.5. The polygon is divided into regions - // defined by they triangle fan formed by the points around the center. This - // is C0 continuous but not necessarily C1 continuous. It is also possible to - // have a non 1 to 1 mapping between parametric coordinates world coordinates - // if the polygon is not planar or convex. - - using FieldType = typename FieldVecType::ComponentType; - - // Find the interpolation for the center point. - FieldType fieldCenter = field[0]; - for (vtkm::IdComponent pointIndex = 1; pointIndex < numPoints; pointIndex++) - { - fieldCenter = fieldCenter + field[pointIndex]; - } - fieldCenter = fieldCenter * FieldType(1.0f / static_cast(numPoints)); - - if ((vtkm::Abs(pcoords[0] - 0.5f) < 4 * vtkm::Epsilon()) && - (vtkm::Abs(pcoords[1] - 0.5f) < 4 * vtkm::Epsilon())) - { - return fieldCenter; - } - - ParametricCoordType angle = vtkm::ATan2(pcoords[1] - 0.5f, pcoords[0] - 0.5f); - if (angle < 0) - { - angle += 2 * vtkm::Pi(); - } - const ParametricCoordType deltaAngle = 2 * vtkm::Pi() / numPoints; - vtkm::IdComponent firstPointIndex = - static_cast(vtkm::Floor(angle / deltaAngle)); - vtkm::IdComponent secondPointIndex = firstPointIndex + 1; - if (secondPointIndex == numPoints) - { - secondPointIndex = 0; - } - - // Transform pcoords for polygon into pcoords for triangle. - vtkm::Vec, 3> polygonCoords; - polygonCoords[0][0] = 0.5f; - polygonCoords[0][1] = 0.5f; - polygonCoords[0][2] = 0; - - polygonCoords[1][0] = - 0.5f * (vtkm::Cos(deltaAngle * static_cast(firstPointIndex)) + 1); - polygonCoords[1][1] = - 0.5f * (vtkm::Sin(deltaAngle * static_cast(firstPointIndex)) + 1); - polygonCoords[1][2] = 0.0f; - - polygonCoords[2][0] = - 0.5f * (vtkm::Cos(deltaAngle * static_cast(secondPointIndex)) + 1); - polygonCoords[2][1] = - 0.5f * (vtkm::Sin(deltaAngle * static_cast(secondPointIndex)) + 1); - polygonCoords[2][2] = 0.0f; - - vtkm::Vec trianglePCoords = - vtkm::exec::internal::ReverseInterpolateTriangle(polygonCoords, pcoords); - - // Set up parameters for triangle that pcoords is in - vtkm::Vec triangleField; - triangleField[0] = fieldCenter; - triangleField[1] = field[firstPointIndex]; - triangleField[2] = field[secondPointIndex]; - - // Now use the triangle interpolate - return vtkm::exec::CellInterpolate( - triangleField, trianglePCoords, vtkm::CellShapeTagTriangle(), worklet); } //----------------------------------------------------------------------------- -template -VTKM_EXEC typename FieldVecType::ComponentType CellInterpolate( - const FieldVecType& field, - const vtkm::Vec& pcoords, - vtkm::CellShapeTagQuad, - const vtkm::exec::FunctorBase& vtkmNotUsed(worklet)) -{ - VTKM_ASSERT(field.GetNumberOfComponents() == 4); - - using T = typename FieldVecType::ComponentType; - - T bottomInterp = vtkm::Lerp(field[0], field[1], pcoords[0]); - T topInterp = vtkm::Lerp(field[3], field[2], pcoords[0]); - - return vtkm::Lerp(bottomInterp, topInterp, pcoords[1]); -} - template VTKM_EXEC vtkm::Vec3f CellInterpolate(const vtkm::VecAxisAlignedPointCoordinates<2>& field, const vtkm::Vec& pcoords, vtkm::CellShapeTagQuad, - const vtkm::exec::FunctorBase&) + const vtkm::exec::FunctorBase& worklet) { - using T = vtkm::Vec3f; - - const T& origin = field.GetOrigin(); - const T& spacing = field.GetSpacing(); - - return T(origin[0] + static_cast(pcoords[0]) * spacing[0], - origin[1] + static_cast(pcoords[1]) * spacing[1], - origin[2]); + return internal::CellInterpolateImpl(vtkc::Pixel{}, field, pcoords, worklet); } //----------------------------------------------------------------------------- -template -VTKM_EXEC typename FieldVecType::ComponentType CellInterpolate( - const FieldVecType& field, - const vtkm::Vec& pcoords, - vtkm::CellShapeTagTetra, - const vtkm::exec::FunctorBase& vtkmNotUsed(worklet)) -{ - VTKM_ASSERT(field.GetNumberOfComponents() == 4); - using T = typename FieldVecType::ComponentType; - return static_cast((field[0] * (1 - pcoords[0] - pcoords[1] - pcoords[2])) + - (field[1] * pcoords[0]) + (field[2] * pcoords[1]) + - (field[3] * pcoords[2])); -} - -//----------------------------------------------------------------------------- -template -VTKM_EXEC typename FieldVecType::ComponentType CellInterpolate( - const FieldVecType& field, - const vtkm::Vec& pcoords, - vtkm::CellShapeTagHexahedron, - const vtkm::exec::FunctorBase& vtkmNotUsed(worklet)) -{ - VTKM_ASSERT(field.GetNumberOfComponents() == 8); - - using T = typename FieldVecType::ComponentType; - - T bottomFrontInterp = vtkm::Lerp(field[0], field[1], pcoords[0]); - T bottomBackInterp = vtkm::Lerp(field[3], field[2], pcoords[0]); - T topFrontInterp = vtkm::Lerp(field[4], field[5], pcoords[0]); - T topBackInterp = vtkm::Lerp(field[7], field[6], pcoords[0]); - - T bottomInterp = vtkm::Lerp(bottomFrontInterp, bottomBackInterp, pcoords[1]); - T topInterp = vtkm::Lerp(topFrontInterp, topBackInterp, pcoords[1]); - - return vtkm::Lerp(bottomInterp, topInterp, pcoords[2]); -} - template VTKM_EXEC vtkm::Vec3f CellInterpolate(const vtkm::VecAxisAlignedPointCoordinates<3>& field, const vtkm::Vec& pcoords, vtkm::CellShapeTagHexahedron, - const vtkm::exec::FunctorBase&) + const vtkm::exec::FunctorBase& worklet) { - vtkm::Vec3f pcoordsCast(static_cast(pcoords[0]), - static_cast(pcoords[1]), - static_cast(pcoords[2])); - - return field.GetOrigin() + pcoordsCast * field.GetSpacing(); -} - -//----------------------------------------------------------------------------- -template -VTKM_EXEC typename FieldVecType::ComponentType CellInterpolate( - const FieldVecType& field, - const vtkm::Vec& pcoords, - vtkm::CellShapeTagWedge, - const vtkm::exec::FunctorBase& vtkmNotUsed(worklet)) -{ - VTKM_ASSERT(field.GetNumberOfComponents() == 6); - - using T = typename FieldVecType::ComponentType; - - T bottomInterp = static_cast((field[0] * (1 - pcoords[0] - pcoords[1])) + - (field[1] * pcoords[1]) + (field[2] * pcoords[0])); - - T topInterp = static_cast((field[3] * (1 - pcoords[0] - pcoords[1])) + - (field[4] * pcoords[1]) + (field[5] * pcoords[0])); - - return vtkm::Lerp(bottomInterp, topInterp, pcoords[2]); -} - -//----------------------------------------------------------------------------- -template -VTKM_EXEC typename FieldVecType::ComponentType CellInterpolate( - const FieldVecType& field, - const vtkm::Vec& pcoords, - vtkm::CellShapeTagPyramid, - const vtkm::exec::FunctorBase& vtkmNotUsed(worklet)) -{ - VTKM_ASSERT(field.GetNumberOfComponents() == 5); - - using T = typename FieldVecType::ComponentType; - - T frontInterp = vtkm::Lerp(field[0], field[1], pcoords[0]); - T backInterp = vtkm::Lerp(field[3], field[2], pcoords[0]); - - T baseInterp = vtkm::Lerp(frontInterp, backInterp, pcoords[1]); - - return vtkm::Lerp(baseInterp, field[4], pcoords[2]); + return internal::CellInterpolateImpl(vtkc::Voxel{}, field, pcoords, worklet); } } } // namespace vtkm::exec diff --git a/vtkm/exec/Jacobian.h b/vtkm/exec/Jacobian.h deleted file mode 100644 index d64aad78e..000000000 --- a/vtkm/exec/Jacobian.h +++ /dev/null @@ -1,239 +0,0 @@ -//============================================================================ -// 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. -//============================================================================ -#ifndef vtk_m_exec_Jacobian_h -#define vtk_m_exec_Jacobian_h - -#include -#include -#include -#include -#include - -namespace vtkm -{ -namespace exec -{ -namespace internal -{ - -template -struct Space2D -{ - using Vec3 = vtkm::Vec; - using Vec2 = vtkm::Vec; - - Vec3 Origin; - Vec3 Basis0; - Vec3 Basis1; - - VTKM_EXEC - Space2D(const Vec3& origin, const Vec3& pointFirst, const Vec3& pointLast) - { - this->Origin = origin; - - this->Basis0 = vtkm::Normal(pointFirst - this->Origin); - - Vec3 n = vtkm::Cross(this->Basis0, pointLast - this->Origin); - this->Basis1 = vtkm::Normal(vtkm::Cross(this->Basis0, n)); - } - - VTKM_EXEC - Vec2 ConvertCoordToSpace(const Vec3 coord) const - { - Vec3 vec = coord - this->Origin; - return Vec2(vtkm::Dot(vec, this->Basis0), vtkm::Dot(vec, this->Basis1)); - } - - template - VTKM_EXEC vtkm::Vec ConvertVecFromSpace(const vtkm::Vec vec) const - { - return vec[0] * this->Basis0 + vec[1] * this->Basis1; - } -}; - -} //namespace internal - -#define VTKM_DERIVATIVE_WEIGHTS_HEXAHEDRON(pc, rc, call) \ - call(0, (-rc[1] * rc[2]), (-rc[0] * rc[2]), (-rc[0] * rc[1])); \ - call(1, (rc[1] * rc[2]), (-pc[0] * rc[2]), (-pc[0] * rc[1])); \ - call(2, (pc[1] * rc[2]), (pc[0] * rc[2]), (-pc[0] * pc[1])); \ - call(3, (-pc[1] * rc[2]), (rc[0] * rc[2]), (-rc[0] * pc[1])); \ - call(4, (-rc[1] * pc[2]), (-rc[0] * pc[2]), (rc[0] * rc[1])); \ - call(5, (rc[1] * pc[2]), (-pc[0] * pc[2]), (pc[0] * rc[1])); \ - call(6, (pc[1] * pc[2]), (pc[0] * pc[2]), (pc[0] * pc[1])); \ - call(7, (-pc[1] * pc[2]), (rc[0] * pc[2]), (rc[0] * pc[1])) - -#define VTKM_DERIVATIVE_WEIGHTS_WEDGE(pc, rc, call) \ - call(0, (-rc[2]), (-rc[2]), (pc[0] - rc[1])); \ - call(1, (0.0f), (rc[2]), (-pc[1])); \ - call(2, (rc[2]), (0.0f), (-pc[0])); \ - call(3, (-pc[2]), (-pc[2]), (rc[0] - pc[1])); \ - call(4, (0.0f), (pc[2]), (pc[1])); \ - call(5, (pc[2]), (0.0f), (pc[0])) - -#define VTKM_DERIVATIVE_WEIGHTS_PYRAMID(pc, rc, call) \ - call(0, (-rc[1] * rc[2]), (-rc[0] * rc[2]), (-rc[0] * rc[1])); \ - call(1, (rc[1] * rc[2]), (-pc[0] * rc[2]), (-pc[0] * rc[1])); \ - call(2, (pc[1] * rc[2]), (pc[0] * rc[2]), (-pc[0] * pc[1])); \ - call(3, (-pc[1] * rc[2]), (rc[0] * rc[2]), (-rc[0] * pc[1])); \ - call(4, (0.0f), (0.0f), (1.0f)) - -#define VTKM_DERIVATIVE_WEIGHTS_QUAD(pc, rc, call) \ - call(0, (-rc[1]), (-rc[0])); \ - call(1, (rc[1]), (-pc[0])); \ - call(2, (pc[1]), (pc[0])); \ - call(3, (-pc[1]), (rc[0])) - -//----------------------------------------------------------------------------- -// This returns the Jacobian of a hexahedron's (or other 3D cell's) coordinates -// with respect to parametric coordinates. Explicitly, this is (d is partial -// derivative): -// -// | | -// | dx/du dx/dv dx/dw | -// | | -// | dy/du dy/dv dy/dw | -// | | -// | dz/du dz/dv dz/dw | -// | | -// - -#define VTKM_ACCUM_JACOBIAN_3D(pointIndex, weight0, weight1, weight2) \ - jacobian(0, 0) += static_cast(wCoords[pointIndex][0] * (weight0)); \ - jacobian(1, 0) += static_cast(wCoords[pointIndex][1] * (weight0)); \ - jacobian(2, 0) += static_cast(wCoords[pointIndex][2] * (weight0)); \ - jacobian(0, 1) += static_cast(wCoords[pointIndex][0] * (weight1)); \ - jacobian(1, 1) += static_cast(wCoords[pointIndex][1] * (weight1)); \ - jacobian(2, 1) += static_cast(wCoords[pointIndex][2] * (weight1)); \ - jacobian(0, 2) += static_cast(wCoords[pointIndex][0] * (weight2)); \ - jacobian(1, 2) += static_cast(wCoords[pointIndex][1] * (weight2)); \ - jacobian(2, 2) += static_cast(wCoords[pointIndex][2] * (weight2)) - -template -VTKM_EXEC void JacobianFor3DCell(const WorldCoordType& wCoords, - const vtkm::Vec& pcoords, - vtkm::Matrix& jacobian, - vtkm::CellShapeTagHexahedron) -{ - vtkm::Vec pc(pcoords); - vtkm::Vec rc = vtkm::Vec(JacobianType(1)) - pc; - - jacobian = vtkm::Matrix(JacobianType(0)); - VTKM_DERIVATIVE_WEIGHTS_HEXAHEDRON(pc, rc, VTKM_ACCUM_JACOBIAN_3D); -} - -template -VTKM_EXEC void JacobianFor3DCell(const WorldCoordType& wCoords, - const vtkm::Vec& pcoords, - vtkm::Matrix& jacobian, - vtkm::CellShapeTagWedge) -{ - vtkm::Vec pc(pcoords); - vtkm::Vec rc = vtkm::Vec(JacobianType(1)) - pc; - - jacobian = vtkm::Matrix(JacobianType(0)); - VTKM_DERIVATIVE_WEIGHTS_WEDGE(pc, rc, VTKM_ACCUM_JACOBIAN_3D); -} - -template -VTKM_EXEC void JacobianFor3DCell(const WorldCoordType& wCoords, - const vtkm::Vec& pcoords, - vtkm::Matrix& jacobian, - vtkm::CellShapeTagPyramid) -{ - vtkm::Vec pc(pcoords); - vtkm::Vec rc = vtkm::Vec(JacobianType(1)) - pc; - - jacobian = vtkm::Matrix(JacobianType(0)); - VTKM_DERIVATIVE_WEIGHTS_PYRAMID(pc, rc, VTKM_ACCUM_JACOBIAN_3D); -} - -// Derivatives in quadrilaterals are computed in much the same way as -// hexahedra. Review the documentation for hexahedra derivatives for details -// on the math. The major difference is that the equations are performed in -// a 2D space built with make_SpaceForQuadrilateral. - -#define VTKM_ACCUM_JACOBIAN_2D(pointIndex, weight0, weight1) \ - wcoords2d = space.ConvertCoordToSpace(wCoords[pointIndex]); \ - jacobian(0, 0) += wcoords2d[0] * (weight0); \ - jacobian(1, 0) += wcoords2d[1] * (weight0); \ - jacobian(0, 1) += wcoords2d[0] * (weight1); \ - jacobian(1, 1) += wcoords2d[1] * (weight1) - -template -VTKM_EXEC void JacobianFor2DCell(const WorldCoordType& wCoords, - const vtkm::Vec& pcoords, - const vtkm::exec::internal::Space2D& space, - vtkm::Matrix& jacobian, - vtkm::CellShapeTagQuad) -{ - vtkm::Vec pc(static_cast(pcoords[0]), - static_cast(pcoords[1])); - vtkm::Vec rc = vtkm::Vec(JacobianType(1)) - pc; - - vtkm::Vec wcoords2d; - jacobian = vtkm::Matrix(JacobianType(0)); - VTKM_DERIVATIVE_WEIGHTS_QUAD(pc, rc, VTKM_ACCUM_JACOBIAN_2D); -} - -#if 0 -// This code doesn't work, so I'm bailing on it. Instead, I'm just grabbing a -// triangle and finding the derivative of that. If you can do better, please -// implement it. -template -VTKM_EXEC -void JacobianFor2DCell(const WorldCoordType &wCoords, - const vtkm::Vec &pcoords, - const vtkm::exec::internal::Space2D &space, - vtkm::Matrix &jacobian, - vtkm::CellShapeTagPolygon) -{ - const vtkm::IdComponent numPoints = wCoords.GetNumberOfComponents(); - vtkm::Vec pc(pcoords[0], pcoords[1]); - JacobianType deltaAngle = 2*vtkm::Pi()/numPoints; - - jacobian = vtkm::Matrix(0); - for (vtkm::IdComponent pointIndex = 0; pointIndex < numPoints; pointIndex++) - { - JacobianType angle = pointIndex*deltaAngle; - vtkm::Vec nodePCoords(0.5f*(vtkm::Cos(angle)+1), - 0.5f*(vtkm::Sin(angle)+1)); - - // This is the vector pointing from the user provided parametric coordinate - // to the node at pointIndex in parametric space. - vtkm::Vec pvec = nodePCoords - pc; - - // The weight (the derivative of the interpolation factor) happens to be - // pvec scaled by the cube root of pvec's magnitude. - JacobianType magSqr = vtkm::MagnitudeSquared(pvec); - JacobianType invMag = vtkm::RSqrt(magSqr); - JacobianType scale = invMag*invMag*invMag; - vtkm::Vec weight = scale*pvec; - - vtkm::Vec wcoords2d = - space.ConvertCoordToSpace(wCoords[pointIndex]); - jacobian(0,0) += wcoords2d[0] * weight[0]; - jacobian(1,0) += wcoords2d[1] * weight[0]; - jacobian(0,1) += wcoords2d[0] * weight[1]; - jacobian(1,1) += wcoords2d[1] * weight[1]; - } -} -#endif - -#undef VTKM_ACCUM_JACOBIAN_3D -#undef VTKM_ACCUM_JACOBIAN_2D -} -} // namespace vtkm::exec -#endif //vtk_m_exec_Jacobian_h diff --git a/vtkm/exec/ParametricCoordinates.h b/vtkm/exec/ParametricCoordinates.h index f90285cfb..e92494afb 100644 --- a/vtkm/exec/ParametricCoordinates.h +++ b/vtkm/exec/ParametricCoordinates.h @@ -12,21 +12,35 @@ #include #include -#include -#include #include #include #include -#include #include #include +#include + namespace vtkm { namespace exec { //----------------------------------------------------------------------------- +template +static inline VTKM_EXEC void ParametricCoordinatesCenter(vtkm::IdComponent numPoints, + vtkm::Vec& pcoords, + CellShapeTag, + const vtkm::exec::FunctorBase&) +{ + auto vtkcTag = typename vtkm::internal::CellShapeTagVtkmToVtkc::Type{}; + + (void)numPoints; // Silence compiler warnings. + VTKM_ASSERT(numPoints == vtkcTag.numberOfPoints()); + + pcoords = vtkm::TypeTraits>::ZeroInitialization(); + vtkc::parametricCenter(vtkcTag, pcoords); +} + template static inline VTKM_EXEC void ParametricCoordinatesCenter(vtkm::IdComponent numPoints, vtkm::Vec& pcoords, @@ -35,9 +49,7 @@ static inline VTKM_EXEC void ParametricCoordinatesCenter(vtkm::IdComponent numPo { (void)numPoints; // Silence compiler warnings. VTKM_ASSERT(numPoints == 0); - pcoords[0] = 0; - pcoords[1] = 0; - pcoords[2] = 0; + pcoords = vtkm::TypeTraits>::ZeroInitialization(); } template @@ -48,22 +60,7 @@ static inline VTKM_EXEC void ParametricCoordinatesCenter(vtkm::IdComponent numPo { (void)numPoints; // Silence compiler warnings. VTKM_ASSERT(numPoints == 1); - pcoords[0] = 0; - pcoords[1] = 0; - pcoords[2] = 0; -} - -template -static inline VTKM_EXEC void ParametricCoordinatesCenter(vtkm::IdComponent numPoints, - vtkm::Vec& pcoords, - vtkm::CellShapeTagLine, - const vtkm::exec::FunctorBase&) -{ - (void)numPoints; // Silence compiler warnings. - VTKM_ASSERT(numPoints == 2); - pcoords[0] = 0.5; - pcoords[1] = 0; - pcoords[2] = 0; + pcoords = vtkm::TypeTraits>::ZeroInitialization(); } template @@ -86,19 +83,6 @@ 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::CellShapeTagTriangle, - const vtkm::exec::FunctorBase&) -{ - (void)numPoints; // Silence compiler warnings. - VTKM_ASSERT(numPoints == 3); - pcoords[0] = static_cast(1.0 / 3.0); - pcoords[1] = static_cast(1.0 / 3.0); - pcoords[2] = 0; -} - template static inline VTKM_EXEC void ParametricCoordinatesCenter(vtkm::IdComponent numPoints, vtkm::Vec& pcoords, @@ -114,82 +98,13 @@ static inline VTKM_EXEC void ParametricCoordinatesCenter(vtkm::IdComponent numPo case 2: ParametricCoordinatesCenter(numPoints, pcoords, vtkm::CellShapeTagLine(), worklet); break; - case 3: - ParametricCoordinatesCenter(numPoints, pcoords, vtkm::CellShapeTagTriangle(), worklet); - break; default: - pcoords[0] = 0.5; - pcoords[1] = 0.5; - pcoords[2] = 0; + pcoords = vtkm::TypeTraits>::ZeroInitialization(); + vtkc::parametricCenter(vtkc::Polygon(numPoints), pcoords); break; } } -template -static inline VTKM_EXEC void ParametricCoordinatesCenter(vtkm::IdComponent numPoints, - vtkm::Vec& pcoords, - vtkm::CellShapeTagQuad, - const vtkm::exec::FunctorBase&) -{ - (void)numPoints; // Silence compiler warnings. - VTKM_ASSERT(numPoints == 4); - pcoords[0] = 0.5; - pcoords[1] = 0.5; - pcoords[2] = 0; -} - -template -static inline VTKM_EXEC void ParametricCoordinatesCenter(vtkm::IdComponent numPoints, - vtkm::Vec& pcoords, - vtkm::CellShapeTagTetra, - const vtkm::exec::FunctorBase&) -{ - (void)numPoints; // Silence compiler warnings. - VTKM_ASSERT(numPoints == 4); - pcoords[0] = 0.25; - pcoords[1] = 0.25; - pcoords[2] = 0.25; -} - -template -static inline VTKM_EXEC void ParametricCoordinatesCenter(vtkm::IdComponent numPoints, - vtkm::Vec& pcoords, - vtkm::CellShapeTagHexahedron, - const vtkm::exec::FunctorBase&) -{ - (void)numPoints; // Silence compiler warnings. - VTKM_ASSERT(numPoints == 8); - pcoords[0] = 0.5; - pcoords[1] = 0.5; - pcoords[2] = 0.5; -} - -template -static inline VTKM_EXEC void ParametricCoordinatesCenter(vtkm::IdComponent numPoints, - vtkm::Vec& pcoords, - vtkm::CellShapeTagWedge, - const vtkm::exec::FunctorBase&) -{ - (void)numPoints; // Silence compiler warnings. - VTKM_ASSERT(numPoints == 6); - pcoords[0] = static_cast(1.0 / 3.0); - pcoords[1] = static_cast(1.0 / 3.0); - pcoords[2] = 0.5; -} - -template -static inline VTKM_EXEC void ParametricCoordinatesCenter(vtkm::IdComponent numPoints, - vtkm::Vec& pcoords, - vtkm::CellShapeTagPyramid, - const vtkm::exec::FunctorBase&) -{ - (void)numPoints; // Silence compiler warnings. - VTKM_ASSERT(numPoints == 5); - pcoords[0] = 0.5; - pcoords[1] = 0.5; - pcoords[2] = static_cast(0.2); -} - //----------------------------------------------------------------------------- /// Returns the parametric center of the given cell shape with the given number /// of points. @@ -220,12 +135,29 @@ static inline VTKM_EXEC vtkm::Vec3f ParametricCoordinatesCenter( CellShapeTag shape, const vtkm::exec::FunctorBase& worklet) { - vtkm::Vec3f pcoords; + vtkm::Vec3f pcoords(0.0f); ParametricCoordinatesCenter(numPoints, pcoords, shape, worklet); return pcoords; } //----------------------------------------------------------------------------- +template +static inline VTKM_EXEC void ParametricCoordinatesPoint(vtkm::IdComponent numPoints, + vtkm::IdComponent pointIndex, + vtkm::Vec& pcoords, + CellShapeTag, + const vtkm::exec::FunctorBase&) +{ + auto vtkcTag = typename vtkm::internal::CellShapeTagVtkmToVtkc::Type{}; + + (void)numPoints; // Silence compiler warnings. + VTKM_ASSUME(numPoints == vtkcTag.numberOfPoints()); + VTKM_ASSUME((pointIndex >= 0) && (pointIndex < numPoints)); + + pcoords = vtkm::TypeTraits>::ZeroInitialization(); + vtkc::parametricPoint(vtkcTag, pointIndex, pcoords); +} + template static inline VTKM_EXEC void ParametricCoordinatesPoint(vtkm::IdComponent, vtkm::IdComponent, @@ -248,25 +180,7 @@ static inline VTKM_EXEC void ParametricCoordinatesPoint(vtkm::IdComponent numPoi (void)pointIndex; // Silence compiler warnings. VTKM_ASSUME(numPoints == 1); VTKM_ASSUME(pointIndex == 0); - pcoords[0] = 0; - pcoords[1] = 0; - pcoords[2] = 0; -} - -template -static inline VTKM_EXEC void ParametricCoordinatesPoint(vtkm::IdComponent numPoints, - vtkm::IdComponent pointIndex, - vtkm::Vec& pcoords, - vtkm::CellShapeTagLine, - const vtkm::exec::FunctorBase&) -{ - (void)numPoints; // Silence compiler warnings. - VTKM_ASSUME(numPoints == 2); - VTKM_ASSUME((pointIndex >= 0) && (pointIndex < 2)); - - pcoords[0] = static_cast(pointIndex); - pcoords[1] = 0; - pcoords[2] = 0; + pcoords = vtkm::TypeTraits>::ZeroInitialization(); } template @@ -292,35 +206,6 @@ 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::CellShapeTagTriangle, - const vtkm::exec::FunctorBase&) -{ - (void)numPoints; // Silence compiler warnings. - VTKM_ASSUME(numPoints == 3); - VTKM_ASSUME((pointIndex >= 0) && (pointIndex < 3)); - - switch (pointIndex) - { - case 0: - pcoords[0] = 0; - pcoords[1] = 0; - break; - case 1: - pcoords[0] = 1; - pcoords[1] = 0; - break; - case 2: - pcoords[0] = 0; - pcoords[1] = 1; - break; - } - pcoords[2] = 0; -} - template static inline VTKM_EXEC void ParametricCoordinatesPoint(vtkm::IdComponent numPoints, vtkm::IdComponent pointIndex, @@ -340,234 +225,10 @@ static inline VTKM_EXEC void ParametricCoordinatesPoint(vtkm::IdComponent numPoi case 2: ParametricCoordinatesPoint(numPoints, pointIndex, pcoords, vtkm::CellShapeTagLine(), worklet); return; - case 3: - ParametricCoordinatesPoint( - numPoints, pointIndex, pcoords, vtkm::CellShapeTagTriangle(), worklet); + default: + pcoords = vtkm::TypeTraits>::ZeroInitialization(); + vtkc::parametricPoint(vtkc::Polygon(numPoints), pointIndex, pcoords); return; - case 4: - ParametricCoordinatesPoint(numPoints, pointIndex, pcoords, vtkm::CellShapeTagQuad(), worklet); - return; - } - - // If we are here, then numPoints >= 5. - - const ParametricCoordType angle = - static_cast(pointIndex * 2 * vtkm::Pi() / numPoints); - - pcoords[0] = 0.5f * (vtkm::Cos(angle) + 1); - pcoords[1] = 0.5f * (vtkm::Sin(angle) + 1); - pcoords[2] = 0; -} - -template -static inline VTKM_EXEC void ParametricCoordinatesPoint(vtkm::IdComponent numPoints, - vtkm::IdComponent pointIndex, - vtkm::Vec& pcoords, - vtkm::CellShapeTagQuad, - const vtkm::exec::FunctorBase&) -{ - (void)numPoints; // Silence compiler warnings. - VTKM_ASSUME(numPoints == 4); - VTKM_ASSUME((pointIndex >= 0) && (pointIndex < 4)); - - switch (pointIndex) - { - case 0: - pcoords[0] = 0; - pcoords[1] = 0; - break; - case 1: - pcoords[0] = 1; - pcoords[1] = 0; - break; - case 2: - pcoords[0] = 1; - pcoords[1] = 1; - break; - case 3: - pcoords[0] = 0; - pcoords[1] = 1; - break; - } - pcoords[2] = 0; -} - -template -static inline VTKM_EXEC void ParametricCoordinatesPoint(vtkm::IdComponent numPoints, - vtkm::IdComponent pointIndex, - vtkm::Vec& pcoords, - vtkm::CellShapeTagTetra, - const vtkm::exec::FunctorBase&) -{ - (void)numPoints; // Silence compiler warnings. - VTKM_ASSUME(numPoints == 4); - VTKM_ASSUME((pointIndex >= 0) && (pointIndex < 4)); - - switch (pointIndex) - { - case 0: - pcoords[0] = 0; - pcoords[1] = 0; - pcoords[2] = 0; - break; - case 1: - pcoords[0] = 1; - pcoords[1] = 0; - pcoords[2] = 0; - break; - case 2: - pcoords[0] = 0; - pcoords[1] = 1; - pcoords[2] = 0; - break; - case 3: - pcoords[0] = 0; - pcoords[1] = 0; - pcoords[2] = 1; - break; - } -} - -template -static inline VTKM_EXEC void ParametricCoordinatesPoint(vtkm::IdComponent numPoints, - vtkm::IdComponent pointIndex, - vtkm::Vec& pcoords, - vtkm::CellShapeTagHexahedron, - const vtkm::exec::FunctorBase&) -{ - (void)numPoints; // Silence compiler warnings. - VTKM_ASSUME(numPoints == 8); - VTKM_ASSUME((pointIndex >= 0) && (pointIndex < 8)); - - switch (pointIndex) - { - case 0: - pcoords[0] = 0; - pcoords[1] = 0; - pcoords[2] = 0; - break; - case 1: - pcoords[0] = 1; - pcoords[1] = 0; - pcoords[2] = 0; - break; - case 2: - pcoords[0] = 1; - pcoords[1] = 1; - pcoords[2] = 0; - break; - case 3: - pcoords[0] = 0; - pcoords[1] = 1; - pcoords[2] = 0; - break; - case 4: - pcoords[0] = 0; - pcoords[1] = 0; - pcoords[2] = 1; - break; - case 5: - pcoords[0] = 1; - pcoords[1] = 0; - pcoords[2] = 1; - break; - case 6: - pcoords[0] = 1; - pcoords[1] = 1; - pcoords[2] = 1; - break; - case 7: - pcoords[0] = 0; - pcoords[1] = 1; - pcoords[2] = 1; - break; - } -} - -template -static inline VTKM_EXEC void ParametricCoordinatesPoint(vtkm::IdComponent numPoints, - vtkm::IdComponent pointIndex, - vtkm::Vec& pcoords, - vtkm::CellShapeTagWedge, - const vtkm::exec::FunctorBase&) -{ - (void)numPoints; // Silence compiler warnings. - VTKM_ASSUME(numPoints == 6); - VTKM_ASSUME((pointIndex >= 0) && (pointIndex < 6)); - - switch (pointIndex) - { - case 0: - pcoords[0] = 0; - pcoords[1] = 0; - pcoords[2] = 0; - break; - case 1: - pcoords[0] = 0; - pcoords[1] = 1; - pcoords[2] = 0; - break; - case 2: - pcoords[0] = 1; - pcoords[1] = 0; - pcoords[2] = 0; - break; - case 3: - pcoords[0] = 0; - pcoords[1] = 0; - pcoords[2] = 1; - break; - case 4: - pcoords[0] = 0; - pcoords[1] = 1; - pcoords[2] = 1; - break; - case 5: - pcoords[0] = 1; - pcoords[1] = 0; - pcoords[2] = 1; - break; - } -} - -template -static inline VTKM_EXEC void ParametricCoordinatesPoint(vtkm::IdComponent numPoints, - vtkm::IdComponent pointIndex, - vtkm::Vec& pcoords, - vtkm::CellShapeTagPyramid, - const vtkm::exec::FunctorBase&) -{ - (void)numPoints; // Silence compiler warnings. - VTKM_ASSUME(numPoints == 5); - VTKM_ASSUME((pointIndex >= 0) && (pointIndex < 5)); - - switch (pointIndex) - { - case 0: - pcoords[0] = 0; - pcoords[1] = 0; - pcoords[2] = 0; - break; - case 1: - pcoords[0] = 1; - pcoords[1] = 0; - pcoords[2] = 0; - break; - case 2: - pcoords[0] = 1; - pcoords[1] = 1; - pcoords[2] = 0; - break; - case 3: - pcoords[0] = 0; - pcoords[1] = 1; - pcoords[2] = 0; - break; - case 4: - pcoords[0] = 0.5; - pcoords[1] = 0.5; - pcoords[2] = 1; - break; } } @@ -603,12 +264,34 @@ static inline VTKM_EXEC vtkm::Vec3f ParametricCoordinatesPoint( CellShapeTag shape, const vtkm::exec::FunctorBase& worklet) { - vtkm::Vec3f pcoords; + vtkm::Vec3f pcoords(0.0f); ParametricCoordinatesPoint(numPoints, pointIndex, pcoords, shape, worklet); return pcoords; } //----------------------------------------------------------------------------- +namespace internal +{ + +template +static inline VTKM_EXEC typename WorldCoordVector::ComponentType +ParametricCoordinatesToWorldCoordinatesImpl(VtkcCellShapeTag tag, + const WorldCoordVector& pointWCoords, + const PCoordType& pcoords, + const vtkm::exec::FunctorBase& worklet) +{ + typename WorldCoordVector::ComponentType wcoords(0); + auto status = vtkc::parametricToWorld( + tag, vtkc::makeFieldAccessorNestedSOA(pointWCoords, 3), pcoords, wcoords); + if (status != vtkc::ErrorCode::SUCCESS) + { + worklet.RaiseError(vtkc::errorString(status)); + } + return wcoords; +} + +} // namespace internal + template static inline VTKM_EXEC typename WorldCoordVector::ComponentType ParametricCoordinatesToWorldCoordinates(const WorldCoordVector& pointWCoords, @@ -616,10 +299,143 @@ ParametricCoordinatesToWorldCoordinates(const WorldCoordVector& pointWCoords, CellShapeTag shape, const vtkm::exec::FunctorBase& worklet) { - return vtkm::exec::CellInterpolate(pointWCoords, pcoords, shape, worklet); + auto numPoints = pointWCoords.GetNumberOfComponents(); + return internal::ParametricCoordinatesToWorldCoordinatesImpl( + vtkm::internal::make_VtkcCellShapeTag(shape, numPoints), pointWCoords, pcoords, worklet); +} + +template +static inline VTKM_EXEC typename WorldCoordVector::ComponentType +ParametricCoordinatesToWorldCoordinates(const WorldCoordVector& pointWCoords, + const vtkm::Vec& pcoords, + vtkm::CellShapeTagEmpty empty, + const vtkm::exec::FunctorBase& worklet) +{ + return vtkm::exec::CellInterpolate(pointWCoords, pcoords, empty, worklet); +} + +template +static inline VTKM_EXEC typename WorldCoordVector::ComponentType +ParametricCoordinatesToWorldCoordinates(const WorldCoordVector& pointWCoords, + const vtkm::Vec& pcoords, + vtkm::CellShapeTagPolyLine polyLine, + const vtkm::exec::FunctorBase& worklet) +{ + return vtkm::exec::CellInterpolate(pointWCoords, pcoords, polyLine, worklet); +} + +template +static inline VTKM_EXEC typename WorldCoordVector::ComponentType +ParametricCoordinatesToWorldCoordinates(const WorldCoordVector& pointWCoords, + const vtkm::Vec& pcoords, + vtkm::CellShapeTagPolygon, + const vtkm::exec::FunctorBase& worklet) +{ + auto numPoints = pointWCoords.GetNumberOfComponents(); + switch (numPoints) + { + case 1: + return ParametricCoordinatesToWorldCoordinates( + pointWCoords, pcoords, vtkm::CellShapeTagVertex{}, worklet); + case 2: + return ParametricCoordinatesToWorldCoordinates( + pointWCoords, pcoords, vtkm::CellShapeTagLine{}, worklet); + default: + return internal::ParametricCoordinatesToWorldCoordinatesImpl( + vtkc::Polygon(numPoints), pointWCoords, pcoords, worklet); + } +} + +template +static inline VTKM_EXEC typename WorldCoordVector::ComponentType +ParametricCoordinatesToWorldCoordinates(const vtkm::VecAxisAlignedPointCoordinates<2>& pointWCoords, + const vtkm::Vec& pcoords, + vtkm::CellShapeTagQuad, + const vtkm::exec::FunctorBase& worklet) +{ + return internal::ParametricCoordinatesToWorldCoordinatesImpl( + vtkc::Pixel{}, pointWCoords, pcoords, worklet); +} + +template +static inline VTKM_EXEC typename WorldCoordVector::ComponentType +ParametricCoordinatesToWorldCoordinates(const vtkm::VecAxisAlignedPointCoordinates<3>& pointWCoords, + const vtkm::Vec& pcoords, + vtkm::CellShapeTagHexahedron, + const vtkm::exec::FunctorBase& worklet) +{ + return internal::ParametricCoordinatesToWorldCoordinatesImpl( + vtkc::Voxel{}, pointWCoords, pcoords, worklet); } //----------------------------------------------------------------------------- +/// Returns the world coordinate corresponding to the given parametric coordinate of a cell. +/// +template +static inline VTKM_EXEC typename WorldCoordVector::ComponentType +ParametricCoordinatesToWorldCoordinates(const WorldCoordVector& pointWCoords, + const vtkm::Vec& pcoords, + vtkm::CellShapeTagGeneric shape, + const vtkm::exec::FunctorBase& worklet) +{ + typename WorldCoordVector::ComponentType wcoords(0); + switch (shape.Id) + { + vtkmGenericCellShapeMacro(wcoords = ParametricCoordinatesToWorldCoordinates( + pointWCoords, pcoords, CellShapeTag(), worklet)); + default: + worklet.RaiseError("Bad shape given to ParametricCoordinatesPoint."); + break; + } + return wcoords; +} + +//----------------------------------------------------------------------------- +namespace internal +{ + +template +static inline VTKM_EXEC typename WorldCoordVector::ComponentType +WorldCoordinatesToParametricCoordinatesImpl(VtkcCellShapeTag tag, + const WorldCoordVector& pointWCoords, + const typename WorldCoordVector::ComponentType& wcoords, + bool& success, + const vtkm::exec::FunctorBase& worklet) +{ + VTKM_ASSERT(pointWCoords.GetNumberOfComponents() == tag.numberOfPoints()); + + auto pcoords = vtkm::TypeTraits::ZeroInitialization(); + auto status = vtkc::worldToParametric( + tag, vtkc::makeFieldAccessorNestedSOA(pointWCoords, 3), wcoords, pcoords); + + success = true; + if (status != vtkc::ErrorCode::SUCCESS) + { + worklet.RaiseError(vtkc::errorString(status)); + success = false; + } + return pcoords; +} + +} // namespace internal + +template +static inline VTKM_EXEC typename WorldCoordVector::ComponentType +WorldCoordinatesToParametricCoordinates(const WorldCoordVector& pointWCoords, + const typename WorldCoordVector::ComponentType& wcoords, + CellShapeTag shape, + bool& success, + const vtkm::exec::FunctorBase& worklet) +{ + auto numPoints = pointWCoords.GetNumberOfComponents(); + return internal::WorldCoordinatesToParametricCoordinatesImpl( + vtkm::internal::make_VtkcCellShapeTag(shape, numPoints), + pointWCoords, + wcoords, + success, + worklet); +} + template static inline VTKM_EXEC typename WorldCoordVector::ComponentType WorldCoordinatesToParametricCoordinates(const WorldCoordVector&, @@ -647,36 +463,6 @@ WorldCoordinatesToParametricCoordinates(const WorldCoordVector& pointWCoords, return typename WorldCoordVector::ComponentType(0, 0, 0); } -template -static inline VTKM_EXEC typename WorldCoordVector::ComponentType -WorldCoordinatesToParametricCoordinates(const WorldCoordVector& pointWCoords, - const typename WorldCoordVector::ComponentType& wcoords, - vtkm::CellShapeTagLine, - bool& success, - const vtkm::exec::FunctorBase& vtkmNotUsed(worklet)) -{ - VTKM_ASSERT(pointWCoords.GetNumberOfComponents() == 2); - success = true; - - // Because this is a line, there is only one valid parametric coordinate. Let - // vec be the vector from the first point to the second point - // (pointWCoords[1] - pointWCoords[0]), which is the direction of the line. - // Dot(vec,wcoords-pointWCoords[0])/mag(vec) is the orthoginal projection of - // wcoords on the line and represents the distance between the orthoginal - // projection and pointWCoords[0]. The parametric coordinate is the fraction - // of this over the length of the segment, which is mag(vec). Thus, the - // parametric coordinate is Dot(vec,wcoords-pointWCoords[0])/mag(vec)^2. - - using Vector3 = typename WorldCoordVector::ComponentType; - using T = typename Vector3::ComponentType; - - Vector3 vec = pointWCoords[1] - pointWCoords[0]; - T numerator = vtkm::Dot(vec, wcoords - pointWCoords[0]); - T denominator = vtkm::MagnitudeSquared(vec); - - return Vector3(numerator / denominator, 0, 0); -} - template static inline VTKM_EXEC typename WorldCoordVector::ComponentType WorldCoordinatesToParametricCoordinates(const WorldCoordVector& pointWCoords, @@ -688,14 +474,10 @@ WorldCoordinatesToParametricCoordinates(const WorldCoordVector& pointWCoords, vtkm::IdComponent numPoints = pointWCoords.GetNumberOfComponents(); VTKM_ASSERT(pointWCoords.GetNumberOfComponents() >= 1); - switch (numPoints) + if (numPoints == 1) { - case 1: - return WorldCoordinatesToParametricCoordinates( - pointWCoords, wcoords, vtkm::CellShapeTagVertex(), success, worklet); - case 2: - return WorldCoordinatesToParametricCoordinates( - pointWCoords, wcoords, vtkm::CellShapeTagLine(), success, worklet); + return WorldCoordinatesToParametricCoordinates( + pointWCoords, wcoords, vtkm::CellShapeTagVertex(), success, worklet); } using Vector3 = typename WorldCoordVector::ComponentType; @@ -720,28 +502,18 @@ WorldCoordinatesToParametricCoordinates(const WorldCoordVector& pointWCoords, //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; + idx = 1; } + vtkm::Vec line(pointWCoords[idx - 1], pointWCoords[idx]); + auto lpc = WorldCoordinatesToParametricCoordinates( + line, wcoords, vtkm::CellShapeTagLine{}, success, worklet); + //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; + T polyLineParam = static_cast(idx - 1) * dParam + lpc[0] * dParam; return Vector3(polyLineParam, 0, 0); } @@ -750,448 +522,50 @@ template static inline VTKM_EXEC typename WorldCoordVector::ComponentType WorldCoordinatesToParametricCoordinates(const WorldCoordVector& pointWCoords, const typename WorldCoordVector::ComponentType& wcoords, - vtkm::CellShapeTagTriangle, + vtkm::CellShapeTagPolygon, bool& success, - const vtkm::exec::FunctorBase&) + const vtkm::exec::FunctorBase& worklet) { - success = true; - vtkm::exec::internal::FastVec local(pointWCoords); - return vtkm::exec::internal::ReverseInterpolateTriangle(local.Get(), wcoords); + auto numPoints = pointWCoords.GetNumberOfComponents(); + switch (numPoints) + { + case 1: + return WorldCoordinatesToParametricCoordinates( + pointWCoords, wcoords, vtkm::CellShapeTagVertex{}, success, worklet); + case 2: + return WorldCoordinatesToParametricCoordinates( + pointWCoords, wcoords, vtkm::CellShapeTagLine{}, success, worklet); + default: + return internal::WorldCoordinatesToParametricCoordinatesImpl( + vtkc::Polygon(numPoints), pointWCoords, wcoords, success, worklet); + } } -//----------------------------------------------------------------------------- -namespace detail -{ - -template -class JacobianFunctorQuad -{ - using T = typename WorldCoordVector::ComponentType::ComponentType; - using Vector2 = vtkm::Vec; - using Matrix2x2 = vtkm::Matrix; - using SpaceType = vtkm::exec::internal::Space2D; - - const WorldCoordVector* PointWCoords; - const SpaceType* Space; - -public: - VTKM_EXEC - JacobianFunctorQuad(const WorldCoordVector* pointWCoords, const SpaceType* space) - : PointWCoords(pointWCoords) - , Space(space) - { - } - - VTKM_EXEC - Matrix2x2 operator()(const Vector2& pcoords) const - { - Matrix2x2 jacobian; - vtkm::exec::JacobianFor2DCell(*this->PointWCoords, - vtkm::Vec(pcoords[0], pcoords[1], 0), - *this->Space, - jacobian, - CellShapeTag()); - return jacobian; - } -}; - -template -class CoordinatesFunctorQuad -{ - using T = typename WorldCoordVector::ComponentType::ComponentType; - using Vector2 = vtkm::Vec; - using Vector3 = vtkm::Vec; - using SpaceType = vtkm::exec::internal::Space2D; - - const WorldCoordVector* PointWCoords; - const SpaceType* Space; - const vtkm::exec::FunctorBase* Worklet; - -public: - VTKM_EXEC - CoordinatesFunctorQuad(const WorldCoordVector* pointWCoords, - const SpaceType* space, - const vtkm::exec::FunctorBase* worklet) - : PointWCoords(pointWCoords) - , Space(space) - , Worklet(worklet) - { - } - - VTKM_EXEC - Vector2 operator()(Vector2 pcoords) const - { - Vector3 pcoords3D(pcoords[0], pcoords[1], 0); - Vector3 wcoords = vtkm::exec::ParametricCoordinatesToWorldCoordinates( - *this->PointWCoords, pcoords3D, CellShapeTag(), *this->Worklet); - return this->Space->ConvertCoordToSpace(wcoords); - } -}; - -template -static inline VTKM_EXEC typename WorldCoordVector::ComponentType -WorldCoordinatesToParametricCoordinatesQuad(const WorldCoordVector& pointWCoords, - const typename WorldCoordVector::ComponentType& wcoords, - bool& success, - const vtkm::exec::FunctorBase& worklet) -{ - using T = typename WorldCoordVector::ComponentType::ComponentType; - using Vector2 = vtkm::Vec; - using Vector3 = vtkm::Vec; - - // We have an underdetermined system in 3D, so create a 2D space in the - // plane that the polygon sits. - vtkm::exec::internal::Space2D space(pointWCoords[0], pointWCoords[1], pointWCoords[3]); - - auto result = vtkm::NewtonsMethod( - JacobianFunctorQuad(&pointWCoords, &space), - CoordinatesFunctorQuad( - &pointWCoords, &space, &worklet), - space.ConvertCoordToSpace(wcoords), - Vector2(0.5f, 0.5f)); - - success = result.Valid; - return Vector3(result.Solution[0], result.Solution[1], 0); -} -} // namespace detail - static inline VTKM_EXEC vtkm::Vec3f WorldCoordinatesToParametricCoordinates( const vtkm::VecAxisAlignedPointCoordinates<2>& pointWCoords, const vtkm::Vec3f& wcoords, vtkm::CellShapeTagQuad, bool& success, - const FunctorBase&) + const FunctorBase& worklet) { - success = true; - return (wcoords - pointWCoords.GetOrigin()) / pointWCoords.GetSpacing(); + return internal::WorldCoordinatesToParametricCoordinatesImpl( + vtkc::Pixel{}, pointWCoords, wcoords, success, worklet); } -template -static inline VTKM_EXEC typename WorldCoordVector::ComponentType -WorldCoordinatesToParametricCoordinates(const WorldCoordVector& pointWCoords, - const typename WorldCoordVector::ComponentType& wcoords, - vtkm::CellShapeTagQuad, - bool& success, - const vtkm::exec::FunctorBase& worklet) -{ - vtkm::exec::internal::FastVec local(pointWCoords); - return detail::WorldCoordinatesToParametricCoordinatesQuad( - local.Get(), wcoords, success, worklet); -} - -//----------------------------------------------------------------------------- -template -static inline VTKM_EXEC typename WorldCoordVector::ComponentType -WorldCoordinatesToParametricCoordinates(const WorldCoordVector& pointWCoords, - const typename WorldCoordVector::ComponentType& wcoords, - vtkm::CellShapeTagPolygon, - bool& success, - const vtkm::exec::FunctorBase& worklet) -{ - const vtkm::IdComponent numPoints = pointWCoords.GetNumberOfComponents(); - VTKM_ASSERT(numPoints > 0); - switch (numPoints) - { - case 1: - return WorldCoordinatesToParametricCoordinates( - pointWCoords, wcoords, vtkm::CellShapeTagVertex(), success, worklet); - case 2: - return WorldCoordinatesToParametricCoordinates( - pointWCoords, wcoords, vtkm::CellShapeTagLine(), success, worklet); - case 3: - return WorldCoordinatesToParametricCoordinates( - pointWCoords, wcoords, vtkm::CellShapeTagTriangle(), success, worklet); - case 4: - return WorldCoordinatesToParametricCoordinates( - pointWCoords, wcoords, vtkm::CellShapeTagQuad(), success, worklet); - } - - // If we are here, then there are 5 or more points on this polygon. - - // Arrange the points such that they are on the circle circumscribed in the - // unit square from 0 to 1. That is, the point are on the circle centered at - // coordinate 0.5,0.5 with radius 0.5. The polygon is divided into regions - // defined by they triangle fan formed by the points around the center. This - // is C0 continuous but not necessarily C1 continuous. It is also possible to - // have a non 1 to 1 mapping between parametric coordinates world coordinates - // if the polygon is not planar or convex. - - using WCoordType = typename WorldCoordVector::ComponentType; - - // Find the position of the center point. - WCoordType wcoordCenter = pointWCoords[0]; - for (vtkm::IdComponent pointIndex = 1; pointIndex < numPoints; pointIndex++) - { - wcoordCenter = wcoordCenter + pointWCoords[pointIndex]; - } - wcoordCenter = wcoordCenter * WCoordType(1.0f / static_cast(numPoints)); - - // Find the normal vector to the polygon. If the polygon is planar, convex, - // and in general position, any three points will give a normal in the same - // direction. Although not perfectly robust, we can reduce the effect of - // non-planar, non-convex, or degenerate polygons by picking three points - // topologically far from each other. Note that we do not care about the - // length of the normal in this case. - WCoordType polygonNormal; - { - WCoordType vec1 = pointWCoords[numPoints / 3] - pointWCoords[0]; - WCoordType vec2 = pointWCoords[2 * numPoints / 3] - pointWCoords[1]; - polygonNormal = vtkm::Cross(vec1, vec2); - } - - // Find which triangle wcoords is located in. We do this by defining the - // equations for the planes through the radial edges and perpendicular to the - // polygon. The point is in the triangle if it is on the correct side of both - // planes. - vtkm::IdComponent firstPointIndex; - vtkm::IdComponent secondPointIndex = 0; - bool foundTriangle = false; - for (firstPointIndex = 0; firstPointIndex < numPoints - 1; firstPointIndex++) - { - WCoordType vecInPlane = pointWCoords[firstPointIndex] - wcoordCenter; - WCoordType planeNormal = vtkm::Cross(polygonNormal, vecInPlane); - typename WCoordType::ComponentType planeOffset = vtkm::Dot(planeNormal, wcoordCenter); - if (vtkm::Dot(planeNormal, wcoords) < planeOffset) - { - // wcoords on wrong side of plane, thus outside of triangle - continue; - } - - secondPointIndex = firstPointIndex + 1; - vecInPlane = pointWCoords[secondPointIndex] - wcoordCenter; - planeNormal = vtkm::Cross(polygonNormal, vecInPlane); - planeOffset = vtkm::Dot(planeNormal, wcoordCenter); - if (vtkm::Dot(planeNormal, wcoords) > planeOffset) - { - // wcoords on wrong side of plane, thus outside of triangle - continue; - } - - foundTriangle = true; - break; - } - if (!foundTriangle) - { - // wcoord was outside of all triangles we checked. It must be inside the - // one triangle we did not check (the one between the first and last - // polygon points). - firstPointIndex = numPoints - 1; - secondPointIndex = 0; - } - - // Build a structure containing the points of the triangle wcoords is in and - // use the triangle version of this function to find the parametric - // coordinates. - vtkm::Vec triangleWCoords; - triangleWCoords[0] = wcoordCenter; - triangleWCoords[1] = pointWCoords[firstPointIndex]; - triangleWCoords[2] = pointWCoords[secondPointIndex]; - - WCoordType trianglePCoords = WorldCoordinatesToParametricCoordinates( - triangleWCoords, wcoords, vtkm::CellShapeTagTriangle(), success, worklet); - - // trianglePCoords is in the triangle's parameter space rather than the - // polygon's parameter space. We can find the polygon's parameter space by - // repurposing the ParametricCoordinatesToWorldCoordinates by using the - // polygon parametric coordinates as a proxy for world coordinates. - triangleWCoords[0] = WCoordType(0.5f, 0.5f, 0); - ParametricCoordinatesPoint( - numPoints, firstPointIndex, triangleWCoords[1], vtkm::CellShapeTagPolygon(), worklet); - ParametricCoordinatesPoint( - numPoints, secondPointIndex, triangleWCoords[2], vtkm::CellShapeTagPolygon(), worklet); - return ParametricCoordinatesToWorldCoordinates( - triangleWCoords, trianglePCoords, vtkm::CellShapeTagTriangle(), worklet); -} - -//----------------------------------------------------------------------------- -namespace detail -{ - -template -static inline VTKM_EXEC typename WorldCoordVector::ComponentType -WorldCoordinatesToParametricCoordinatesTetra( - const WorldCoordVector& pointWCoords, - const typename WorldCoordVector::ComponentType& wcoords) -{ - // We solve the world to parametric coordinates problem for tetrahedra - // similarly to that for triangles. Before understanding this code, you - // should understand the triangle code (in ReverseInterpolateTriangle in - // CellInterpolate.h). Go ahead. Read it now. - // - // The tetrahedron code is an obvious extension of the triangle code by - // considering the parallelpiped formed by wcoords and p0 of the triangle - // and the three adjacent faces. This parallelpiped is equivalent to the - // axis-aligned cuboid anchored at the origin of parametric space. - // - // Just like the triangle, we compute the parametric coordinate for each axis - // by intersecting a plane with each edge emanating from p0. The plane is - // defined by the one that goes through wcoords (duh) and is parallel to the - // plane formed by the other two edges emanating from p0 (as dictated by the - // aforementioned parallelpiped). - // - // In review, by parameterizing the line by fraction of distance the distance - // from p0 to the adjacent point (which is itself the parametric coordinate - // we are after), we get the following definition for the intersection. - // - // d = Dot((wcoords - p0), planeNormal)/Dot((p1-p0), planeNormal) - // - - const auto vec0 = pointWCoords[1] - pointWCoords[0]; - const auto vec1 = pointWCoords[2] - pointWCoords[0]; - const auto vec2 = pointWCoords[3] - pointWCoords[0]; - const auto coordVec = wcoords - pointWCoords[0]; - - typename WorldCoordVector::ComponentType pcoords; - auto planeNormal = vtkm::Cross(vec1, vec2); - pcoords[0] = vtkm::Dot(coordVec, planeNormal) / vtkm::Dot(vec0, planeNormal); - - planeNormal = vtkm::Cross(vec0, vec2); - pcoords[1] = vtkm::Dot(coordVec, planeNormal) / vtkm::Dot(vec1, planeNormal); - - planeNormal = vtkm::Cross(vec0, vec1); - pcoords[2] = vtkm::Dot(coordVec, planeNormal) / vtkm::Dot(vec2, planeNormal); - - return pcoords; -} -} // detail - -template -static inline VTKM_EXEC typename WorldCoordVector::ComponentType -WorldCoordinatesToParametricCoordinates(const WorldCoordVector& pointWCoords, - const typename WorldCoordVector::ComponentType& wcoords, - vtkm::CellShapeTagTetra, - bool& success, - const vtkm::exec::FunctorBase&) -{ - success = true; - vtkm::exec::internal::FastVec local(pointWCoords); - return detail::WorldCoordinatesToParametricCoordinatesTetra(local.Get(), wcoords); -} - -//----------------------------------------------------------------------------- -namespace detail -{ - -template -class JacobianFunctor3DCell -{ - using T = typename WorldCoordVector::ComponentType::ComponentType; - using Vector3 = vtkm::Vec; - using Matrix3x3 = vtkm::Matrix; - - const WorldCoordVector* PointWCoords; - -public: - VTKM_EXEC - JacobianFunctor3DCell(const WorldCoordVector* pointWCoords) - : PointWCoords(pointWCoords) - { - } - - VTKM_EXEC - Matrix3x3 operator()(const Vector3& pcoords) const - { - Matrix3x3 jacobian; - vtkm::exec::JacobianFor3DCell(*this->PointWCoords, pcoords, jacobian, CellShapeTag()); - return jacobian; - } -}; - -template -class CoordinatesFunctor3DCell -{ - using T = typename WorldCoordVector::ComponentType::ComponentType; - using Vector3 = vtkm::Vec; - - const WorldCoordVector* PointWCoords; - const vtkm::exec::FunctorBase* Worklet; - -public: - VTKM_EXEC - CoordinatesFunctor3DCell(const WorldCoordVector* pointWCoords, - const vtkm::exec::FunctorBase* worklet) - : PointWCoords(pointWCoords) - , Worklet(worklet) - { - } - - VTKM_EXEC - Vector3 operator()(Vector3 pcoords) const - { - return vtkm::exec::ParametricCoordinatesToWorldCoordinates( - *this->PointWCoords, pcoords, CellShapeTag(), *this->Worklet); - } -}; - -template -static inline VTKM_EXEC typename WorldCoordVector::ComponentType -WorldCoordinatesToParametricCoordinates3D(const WorldCoordVector& pointWCoords, - const typename WorldCoordVector::ComponentType& wcoords, - CellShapeTag, - bool& success, - const vtkm::exec::FunctorBase& worklet) -{ - auto result = vtkm::NewtonsMethod( - JacobianFunctor3DCell(&pointWCoords), - CoordinatesFunctor3DCell(&pointWCoords, &worklet), - wcoords, - typename WorldCoordVector::ComponentType(0.5f, 0.5f, 0.5f)); - success = result.Valid; - return result.Solution; -} -} // detail - static inline VTKM_EXEC vtkm::Vec3f WorldCoordinatesToParametricCoordinates( const vtkm::VecAxisAlignedPointCoordinates<3>& pointWCoords, const vtkm::Vec3f& wcoords, vtkm::CellShapeTagHexahedron, bool& success, - const FunctorBase&) + const FunctorBase& worklet) { - success = true; - return (wcoords - pointWCoords.GetOrigin()) / pointWCoords.GetSpacing(); -} - -template -static inline VTKM_EXEC typename WorldCoordVector::ComponentType -WorldCoordinatesToParametricCoordinates(const WorldCoordVector& pointWCoords, - const typename WorldCoordVector::ComponentType& wcoords, - vtkm::CellShapeTagHexahedron, - bool& success, - const vtkm::exec::FunctorBase& worklet) -{ - vtkm::exec::internal::FastVec local(pointWCoords); - return detail::WorldCoordinatesToParametricCoordinates3D( - local.Get(), wcoords, vtkm::CellShapeTagHexahedron(), success, worklet); -} - -template -static inline VTKM_EXEC typename WorldCoordVector::ComponentType -WorldCoordinatesToParametricCoordinates(const WorldCoordVector& pointWCoords, - const typename WorldCoordVector::ComponentType& wcoords, - vtkm::CellShapeTagWedge, - bool& success, - const vtkm::exec::FunctorBase& worklet) -{ - vtkm::exec::internal::FastVec local(pointWCoords); - return detail::WorldCoordinatesToParametricCoordinates3D( - local.Get(), wcoords, vtkm::CellShapeTagWedge(), success, worklet); -} - -template -static inline VTKM_EXEC typename WorldCoordVector::ComponentType -WorldCoordinatesToParametricCoordinates(const WorldCoordVector& pointWCoords, - const typename WorldCoordVector::ComponentType& wcoords, - vtkm::CellShapeTagPyramid, - bool& success, - const vtkm::exec::FunctorBase& worklet) -{ - vtkm::exec::internal::FastVec local(pointWCoords); - return detail::WorldCoordinatesToParametricCoordinates3D( - local.Get(), wcoords, vtkm::CellShapeTagPyramid(), success, worklet); + return internal::WorldCoordinatesToParametricCoordinatesImpl( + vtkc::Voxel{}, pointWCoords, wcoords, success, worklet); } //----------------------------------------------------------------------------- +/// Returns the world paramteric corresponding to the given world coordinate for a cell. +/// template static inline VTKM_EXEC typename WorldCoordVector::ComponentType WorldCoordinatesToParametricCoordinates(const WorldCoordVector& pointWCoords, diff --git a/vtkm/thirdparty/vtkc/CMakeLists.txt b/vtkm/thirdparty/vtkc/CMakeLists.txt new file mode 100644 index 000000000..c142cc8a6 --- /dev/null +++ b/vtkm/thirdparty/vtkc/CMakeLists.txt @@ -0,0 +1,28 @@ +##============================================================================ +## 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. +##============================================================================ +add_library(vtkm_vtkc INTERFACE) + +vtkm_get_kit_name(kit_name kit_dir) + +# vtkc needs C++11 +target_compile_features(vtkm_vtkc INTERFACE cxx_std_11) + +target_include_directories(vtkm_vtkc INTERFACE + $ + $) + +install(TARGETS vtkm_vtkc + EXPORT ${VTKm_EXPORT_NAME}) + +## Install headers +if(NOT VTKm_INSTALL_ONLY_LIBRARIES) + install(DIRECTORY vtkmvtkc + DESTINATION ${VTKm_INSTALL_INCLUDE_DIR}/${kit_dir}/) +endif() diff --git a/vtkm/worklet/testing/UnitTestParticleAdvection.cxx b/vtkm/worklet/testing/UnitTestParticleAdvection.cxx index 4f684a5df..4794c918e 100644 --- a/vtkm/worklet/testing/UnitTestParticleAdvection.cxx +++ b/vtkm/worklet/testing/UnitTestParticleAdvection.cxx @@ -317,7 +317,7 @@ void ValidateEvaluator(const EvalType& eval, Status status = statusPortal.Get(index); vtkm::Vec3f result = resultsPortal.Get(index); VTKM_TEST_ASSERT(status == Status::SUCCESS, "Error in evaluator for " + msg); - VTKM_TEST_ASSERT(result == vec, "Error in evaluator result for " + msg); + VTKM_TEST_ASSERT(test_equal(result, vec), "Error in evaluator result for " + msg); } pointsHandle.ReleaseResources(); evalStatus.ReleaseResources(); @@ -370,10 +370,10 @@ void ValidateIntegrator(const IntegratorType& integrator, vtkm::Vec3f result = resultsPortal.Get(index); VTKM_TEST_ASSERT(status != Status::FAIL, "Error in evaluator for " + msg); if (status == Status::OUTSIDE_SPATIAL_BOUNDS) - VTKM_TEST_ASSERT(result == pointsPortal.Get(index), + VTKM_TEST_ASSERT(test_equal(result, pointsPortal.Get(index)), "Error in evaluator result for [OUTSIDE SPATIAL]" + msg); else - VTKM_TEST_ASSERT(result == expStepResults[(size_t)index], + VTKM_TEST_ASSERT(test_equal(result, expStepResults[(size_t)index]), "Error in evaluator result for " + msg); } pointsHandle.ReleaseResources(); diff --git a/vtkm/worklet/testing/UnitTestTemporalAdvection.cxx b/vtkm/worklet/testing/UnitTestTemporalAdvection.cxx index 6966f4dba..d62a9bab3 100644 --- a/vtkm/worklet/testing/UnitTestTemporalAdvection.cxx +++ b/vtkm/worklet/testing/UnitTestTemporalAdvection.cxx @@ -84,7 +84,7 @@ void ValidateEvaluator(const EvalType& eval, vtkm::Vec result = resultsPortal.Get(index); vtkm::Vec expected = validityPortal.Get(index); VTKM_TEST_ASSERT(status == Status::SUCCESS, "Error in evaluator for " + msg); - VTKM_TEST_ASSERT(result == expected, "Error in evaluator result for " + msg); + VTKM_TEST_ASSERT(test_equal(result, expected), "Error in evaluator result for " + msg); } evalStatus.ReleaseResources(); evalResults.ReleaseResources();