From 2e05c1dd03070cb3663e1e1e14c6ae3bcf1d5e4d Mon Sep 17 00:00:00 2001 From: Kenneth Moreland Date: Mon, 5 Dec 2016 17:18:46 -0700 Subject: [PATCH] Support derivatives of vectors The cell derivative/gradient functions were all designed with scalars in mind. Although the field type is templated and you could pass in a vector type for the field, many of these classes would perform the computation incorrectly. These changes specifically support derivatives of vector types. --- vtkm/exec/CellDerivative.h | 484 +++++++++++++++++-- vtkm/exec/testing/UnitTestCellDerivative.cxx | 183 +++++-- 2 files changed, 564 insertions(+), 103 deletions(-) diff --git a/vtkm/exec/CellDerivative.h b/vtkm/exec/CellDerivative.h index b410d97c6..8b60e903f 100644 --- a/vtkm/exec/CellDerivative.h +++ b/vtkm/exec/CellDerivative.h @@ -34,6 +34,51 @@ namespace vtkm { namespace exec { +namespace detail { + +template +struct BaseComponentOfImpl; + +template +struct BaseComponentOfImpl +{ +private: + using ComponentType = typename vtkm::VecTraits::ComponentType; +public: + using type = + typename BaseComponentOfImpl< + ComponentType, + typename vtkm::TypeTraits::DimensionalityTag + >::type; +}; + +template +struct BaseComponentOfImpl + : BaseComponentOfImpl +{ }; + +template +struct BaseComponentOfImpl +{ + using type = ScalarType; +}; + +} // namespace detail + +/// Finds the base component type of a Vec. If you have a Vec of Vecs, it will +/// descend all Vecs until you get to the scalar type. +/// +// If this becomes useful outside of CellDerivative, then this should probably +// be moved to a different header file. +template +struct BaseComponentOf +{ + using type = + typename detail::BaseComponentOfImpl< + VecType, + typename vtkm::TypeTraits::DimensionalityTag>::type; +}; + // 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 @@ -254,6 +299,105 @@ CellDerivativeFor3DCell(const FieldVecType &field, valid); } +template +VTKM_EXEC +vtkm::Vec +CellDerivativeFor2DCellFinish( + const vtkm::Vec &field, + const vtkm::Matrix &LUFactorization, + const vtkm::Vec &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::Vec &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; + + 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::Vec &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 &pcoords, CellShapeTag) { - typedef typename FieldVecType::ComponentType FieldType; + using FieldType = typename FieldVecType::ComponentType; + using BaseFieldType = typename BaseComponentOf::type; // We have an underdetermined system in 3D, so create a 2D space in the // plane that the polygon sits. - vtkm::exec::internal::Space2D space( + 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::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]. - vtkm::Vec parametricDerivative = - ParametricDerivative(field,pcoords,CellShapeTag()); + // 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 @@ -298,10 +445,27 @@ CellDerivativeFor2DCell(const FieldVecType &field, // world space. bool valid; // Ignored. - vtkm::Vec gradient2D = - vtkm::SolveLinearSystem(jacobianTranspose, parametricDerivative, valid); - - return space.ConvertVecFromSpace(gradient2D); + // 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::Vec 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 @@ -379,6 +543,65 @@ CellDerivative(const FieldVecType &field, } //----------------------------------------------------------------------------- +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; + 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 @@ -393,18 +616,17 @@ CellDerivative(const FieldVecType &field, VTKM_ASSERT(field.GetNumberOfComponents() == 2); VTKM_ASSERT(wCoords.GetNumberOfComponents() == 2); - typedef typename FieldVecType::ComponentType FieldType; - typedef vtkm::Vec GradientType; + using FieldType = typename FieldVecType::ComponentType; + using BaseComponentType = typename BaseComponentOf::type; - // 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. + FieldType deltaField(field[1] - field[0]); + vtkm::Vec vec(wCoords[1] - wCoords[0]); - FieldType deltaField = field[1] - field[0]; - GradientType vec = wCoords[1] - wCoords[0]; - - return (deltaField/vtkm::MagnitudeSquared(vec))*vec; + return detail::CellDerivativeLineImpl( + deltaField, + vec, + vtkm::MagnitudeSquared(vec), + typename vtkm::TypeTraits::DimensionalityTag()); } template((field[1]-field[0])/wCoords.GetSpacing()[0], 0, 0); + return vtkm::Vec( + (field[1]-field[0])/wCoords.GetSpacing()[0], T(0), T(0)); } //----------------------------------------------------------------------------- +namespace detail { + +template +VTKM_EXEC +vtkm::Vec +TriangleDerivativeFinish(const vtkm::Vec &field, + const vtkm::Matrix &LUFactorization, + const vtkm::Vec &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::Vec &LUPermutation, + vtkm::TypeTraitsVectorTag) +{ + using FieldTraits = vtkm::VecTraits; + using FieldComponentType = typename FieldTraits::ComponentType; + + vtkm::Vec gradient; + + 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::Vec &LUPermutation, + vtkm::TypeTraitsMatrixTag) +{ + return TriangleDerivativeFinish( + field, LUFactorization, LUPermutation, vtkm::TypeTraitsVectorTag()); +} + template VTKM_EXEC @@ -432,7 +724,7 @@ vtkm::Vec TriangleDerivative(const vtkm::Vec &field, const vtkm::Vec &wCoords) { - typedef vtkm::Vec GradientType; + using BaseComponentType = typename BaseComponentOf::type; // The scalar values of the three points in a triangle completely specify a // linear field (with constant gradient) assuming the field is constant in @@ -458,33 +750,41 @@ TriangleDerivative(const vtkm::Vec &field, // are the differences in points and normal, b has the scalar differences, // and x is really the gradient g. - GradientType v0 = wCoords[1] - wCoords[0]; - GradientType v1 = wCoords[2] - wCoords[0]; - GradientType n = vtkm::Cross(v0, v1); + 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::Matrix A; vtkm::MatrixSetRow(A, 0, v0); vtkm::MatrixSetRow(A, 1, v1); vtkm::MatrixSetRow(A, 2, n); - GradientType b(field[1] - field[0], field[2] - field[0], ValueType(0)); - - // If we want to later change this method to take the gradient of multiple - // values (for example, to find the Jacobian of a vector field), then there - // are more efficient ways to solve them all than independently solving this - // equation for each component of the field. You could find the inverse of - // matrix A. Or you could alter the functions in vtkm/Matrix.h to - // simultaneously solve multiple equations. - // 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; - return vtkm::SolveLinearSystem(A, b, 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::Vec 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 field; inputField.CopyInto(field); vtkm::Vec wpoints; wCoords.CopyInto(wpoints); - return TriangleDerivative(field, wpoints); + return detail::TriangleDerivative(field, wpoints); } //----------------------------------------------------------------------------- @@ -587,7 +887,7 @@ PolygonDerivative(const FieldVecType &field, // Now use the triangle derivative. pcoords is actually invalid for the // triangle, but that does not matter as the derivative for a triangle does // not depend on it. - return TriangleDerivative(triangleField, triangleWCoords); + return detail::TriangleDerivative(triangleField, triangleWCoords); } //----------------------------------------------------------------------------- @@ -690,6 +990,76 @@ CellDerivative(const FieldVecType &field, } //----------------------------------------------------------------------------- +namespace detail { + +template +VTKM_EXEC +vtkm::Vec +TetraDerivativeFinish(const vtkm::Vec &field, + const vtkm::Matrix &LUFactorization, + const vtkm::Vec &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::Vec &LUPermutation, + vtkm::TypeTraitsVectorTag) +{ + using FieldTraits = vtkm::VecTraits; + using FieldComponentType = typename FieldTraits::ComponentType; + + vtkm::Vec gradient; + + 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::Vec &LUPermutation, + vtkm::TypeTraitsMatrixTag) +{ + return TetraDerivativeFinish( + field, LUFactorization, LUPermutation, vtkm::TypeTraitsVectorTag()); +} + template VTKM_EXEC @@ -697,7 +1067,7 @@ vtkm::Vec TetraDerivative(const vtkm::Vec &field, const vtkm::Vec &wCoords) { - typedef vtkm::Vec GradientType; + using BaseComponentType = typename BaseComponentOf::type; // 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 @@ -722,33 +1092,41 @@ TetraDerivative(const vtkm::Vec &field, // are the differences in points and normal, b has the scalar differences, // and x is really the gradient g. - GradientType v0 = wCoords[1] - wCoords[0]; - GradientType v1 = wCoords[2] - wCoords[0]; - GradientType v2 = wCoords[3] - wCoords[0]; + 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::Matrix A; vtkm::MatrixSetRow(A, 0, v0); vtkm::MatrixSetRow(A, 1, v1); vtkm::MatrixSetRow(A, 2, v2); - GradientType b(field[1]-field[0], field[2]-field[0], field[3]-field[0]); - - // If we want to later change this method to take the gradient of multiple - // values (for example, to find the Jacobian of a vector field), then there - // are more efficient ways to solve them all than independently solving this - // equation for each component of the field. You could find the inverse of - // matrix A. Or you could alter the functions in vtkm/Matrix.h to - // simultaneously solve multiple equations. - // 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; - return vtkm::SolveLinearSystem(A, b, 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::Vec 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 field; inputField.CopyInto(field); vtkm::Vec wpoints; wCoords.CopyInto(wpoints); - return TetraDerivative(field, wpoints); + return detail::TetraDerivative(field, wpoints); } //----------------------------------------------------------------------------- diff --git a/vtkm/exec/testing/UnitTestCellDerivative.cxx b/vtkm/exec/testing/UnitTestCellDerivative.cxx index c46614985..bd21a4c56 100644 --- a/vtkm/exec/testing/UnitTestCellDerivative.cxx +++ b/vtkm/exec/testing/UnitTestCellDerivative.cxx @@ -51,14 +51,18 @@ vtkm::Vec WorldToParametric(const vtkm::Vec &wcoord) /// Simple structure describing a linear field. Has a convienience class /// for getting values. +template struct LinearField { - vtkm::Vec Gradient; - vtkm::Float64 OriginValue; + vtkm::Vec Gradient; + FieldType OriginValue; template - T GetValue(vtkm::Vec coordinates) const { - return vtkm::dot(coordinates, - vtkm::Vec(this->Gradient)) + T(this->OriginValue); + FieldType GetValue(vtkm::Vec coordinates) const + { + return ( ( coordinates[0]*this->Gradient[0] + + coordinates[1]*this->Gradient[1] + + coordinates[2]*this->Gradient[2] ) + + this->OriginValue ); } }; @@ -92,7 +96,7 @@ struct TestDerivativeFunctor template void DoTestWithWCoords(CellShapeTag shape, const WCoordsVecType worldCoordinates, - LinearField field, + LinearField field, vtkm::Vec expectedGradient) const { // Stuff to fake running in the execution environment. @@ -159,7 +163,7 @@ struct TestDerivativeFunctor template void DoTest(CellShapeTag shape, vtkm::IdComponent numPoints, - LinearField field, + LinearField field, vtkm::Vec expectedGradient) const { // Stuff to fake running in the execution environment. @@ -192,9 +196,12 @@ struct TestDerivativeFunctor vtkm::IdComponent numPoints, vtkm::IdComponent topDim) const { - LinearField field; + LinearField field; vtkm::Vec expectedGradient; + using FieldTraits = vtkm::VecTraits; + using FieldComponentType = typename FieldTraits::ComponentType; + // Correct topDim for polygons with fewer than 3 points. if (topDim > numPoints-1) { @@ -202,38 +209,70 @@ struct TestDerivativeFunctor } std::cout << "Simple field, " << numPoints << " points" << std::endl; - field.OriginValue = 0.0; - field.Gradient = vtkm::make_Vec(1.0, 1.0, 1.0); - expectedGradient[0] = static_cast((topDim > 0) ? field.Gradient[0] : 0); - expectedGradient[1] = static_cast((topDim > 1) ? field.Gradient[1] : 0); - expectedGradient[2] = static_cast((topDim > 2) ? field.Gradient[2] : 0); + for (vtkm::IdComponent fieldComponent = 0; + fieldComponent < FieldTraits::GetNumberOfComponents(FieldType()); + fieldComponent++) + { + FieldTraits::SetComponent(field.OriginValue, fieldComponent, 0.0); + } + field.Gradient = vtkm::make_Vec(FieldType(1.0f), FieldType(1.0f), FieldType(1.0f)); + expectedGradient[0] = ((topDim > 0) ? field.Gradient[0] : FieldType(0)); + expectedGradient[1] = ((topDim > 1) ? field.Gradient[1] : FieldType(0)); + expectedGradient[2] = ((topDim > 2) ? field.Gradient[2] : FieldType(0)); this->DoTest(shape, numPoints, field, expectedGradient); std::cout << "Uneven gradient, " << numPoints << " points" << std::endl; - field.OriginValue = -7.0; - field.Gradient = vtkm::make_Vec(0.25, 14.0, 11.125); - expectedGradient[0] = static_cast((topDim > 0) ? field.Gradient[0] : 0); - expectedGradient[1] = static_cast((topDim > 1) ? field.Gradient[1] : 0); - expectedGradient[2] = static_cast((topDim > 2) ? field.Gradient[2] : 0); + for (vtkm::IdComponent fieldComponent = 0; + fieldComponent < FieldTraits::GetNumberOfComponents(FieldType()); + fieldComponent++) + { + FieldTraits::SetComponent(field.OriginValue, + fieldComponent, + FieldComponentType(-7.0f)); + } + field.Gradient = vtkm::make_Vec(FieldType(0.25f), FieldType(14.0f), FieldType(11.125f)); + expectedGradient[0] = ((topDim > 0) ? field.Gradient[0] : FieldType(0)); + expectedGradient[1] = ((topDim > 1) ? field.Gradient[1] : FieldType(0)); + expectedGradient[2] = ((topDim > 2) ? field.Gradient[2] : FieldType(0)); this->DoTest(shape, numPoints, field, expectedGradient); std::cout << "Negative gradient directions, " << numPoints << " points" << std::endl; - field.OriginValue = 5.0; - field.Gradient = vtkm::make_Vec(-11.125, -0.25, 14.0); - expectedGradient[0] = static_cast((topDim > 0) ? field.Gradient[0] : 0); - expectedGradient[1] = static_cast((topDim > 1) ? field.Gradient[1] : 0); - expectedGradient[2] = static_cast((topDim > 2) ? field.Gradient[2] : 0); + for (vtkm::IdComponent fieldComponent = 0; + fieldComponent < FieldTraits::GetNumberOfComponents(FieldType()); + fieldComponent++) + { + FieldTraits::SetComponent(field.OriginValue, + fieldComponent, + FieldComponentType(5.0f)); + } + field.Gradient = vtkm::make_Vec(FieldType(-11.125f), FieldType(-0.25f), FieldType(14.0f)); + expectedGradient[0] = ((topDim > 0) ? field.Gradient[0] : FieldType(0)); + expectedGradient[1] = ((topDim > 1) ? field.Gradient[1] : FieldType(0)); + expectedGradient[2] = ((topDim > 2) ? field.Gradient[2] : FieldType(0)); this->DoTest(shape, numPoints, field, expectedGradient); std::cout << "Random linear field, " << numPoints << " points" << std::endl; - std::uniform_real_distribution randomDist(-20.0, 20.0); - field.OriginValue = randomDist(g_RandomGenerator); - field.Gradient = vtkm::make_Vec(randomDist(g_RandomGenerator), - randomDist(g_RandomGenerator), - randomDist(g_RandomGenerator)); - expectedGradient[0] = static_cast((topDim > 0) ? field.Gradient[0] : 0); - expectedGradient[1] = static_cast((topDim > 1) ? field.Gradient[1] : 0); - expectedGradient[2] = static_cast((topDim > 2) ? field.Gradient[2] : 0); + std::uniform_real_distribution randomDist(-20.0f, 20.0f); + for (vtkm::IdComponent fieldComponent = 0; + fieldComponent < FieldTraits::GetNumberOfComponents(FieldType()); + fieldComponent++) + { + FieldTraits::SetComponent(field.OriginValue, + fieldComponent, + randomDist(g_RandomGenerator)); + FieldTraits::SetComponent(field.Gradient[0], + fieldComponent, + randomDist(g_RandomGenerator)); + FieldTraits::SetComponent(field.Gradient[1], + fieldComponent, + randomDist(g_RandomGenerator)); + FieldTraits::SetComponent(field.Gradient[2], + fieldComponent, + randomDist(g_RandomGenerator)); + } + expectedGradient[0] = ((topDim > 0) ? field.Gradient[0] : FieldType(0)); + expectedGradient[1] = ((topDim > 1) ? field.Gradient[1] : FieldType(0)); + expectedGradient[2] = ((topDim > 2) ? field.Gradient[2] : FieldType(0)); this->DoTest(shape, numPoints, field, expectedGradient); } @@ -285,38 +324,82 @@ void TestDerivative() vtkm::testing::Testing::TryAllCellShapes(TestDerivativeFunctor()); std::cout << "======== Float64 ==========================" << std::endl; vtkm::testing::Testing::TryAllCellShapes(TestDerivativeFunctor()); + std::cout << "======== Vec ===================" << std::endl; + vtkm::testing::Testing::TryAllCellShapes(TestDerivativeFunctor >()); + std::cout << "======== Vec ===================" << std::endl; + vtkm::testing::Testing::TryAllCellShapes(TestDerivativeFunctor >()); std::uniform_real_distribution randomDist(-20.0, 20.0); - LinearField field; - field.OriginValue = randomDist(g_RandomGenerator); - field.Gradient = vtkm::make_Vec(randomDist(g_RandomGenerator), - randomDist(g_RandomGenerator), - randomDist(g_RandomGenerator)); - vtkm::Vec expectedGradient = field.Gradient; - - TestDerivativeFunctor testFunctor; vtkm::Vec origin = vtkm::Vec(0.25f, 0.25f, 0.25f); vtkm::Vec spacing = vtkm::Vec(2.0f, 2.0f, 2.0f); + + LinearField scalarField; + scalarField.OriginValue = randomDist(g_RandomGenerator); + scalarField.Gradient = vtkm::make_Vec(randomDist(g_RandomGenerator), + randomDist(g_RandomGenerator), + randomDist(g_RandomGenerator)); + vtkm::Vec expectedScalarGradient = scalarField.Gradient; + + TestDerivativeFunctor testFunctorScalar; std::cout << "======== Uniform Point Coordinates 3D =====" << std::endl; - testFunctor.DoTestWithWCoords( + testFunctorScalar.DoTestWithWCoords( vtkm::CellShapeTagHexahedron(), vtkm::VecRectilinearPointCoordinates<3>(origin, spacing), - field, - expectedGradient); + scalarField, + expectedScalarGradient); std::cout << "======== Uniform Point Coordinates 2D =====" << std::endl; - expectedGradient[2] = 0.0; - testFunctor.DoTestWithWCoords( + expectedScalarGradient[2] = 0.0; + testFunctorScalar.DoTestWithWCoords( vtkm::CellShapeTagQuad(), vtkm::VecRectilinearPointCoordinates<2>(origin, spacing), - field, - expectedGradient); + scalarField, + expectedScalarGradient); std::cout << "======== Uniform Point Coordinates 1D =====" << std::endl; - expectedGradient[1] = 0.0; - testFunctor.DoTestWithWCoords( + expectedScalarGradient[1] = 0.0; + testFunctorScalar.DoTestWithWCoords( vtkm::CellShapeTagLine(), vtkm::VecRectilinearPointCoordinates<1>(origin, spacing), - field, - expectedGradient); + scalarField, + expectedScalarGradient); + + LinearField > vectorField; + vectorField.OriginValue = vtkm::make_Vec(randomDist(g_RandomGenerator), + randomDist(g_RandomGenerator), + randomDist(g_RandomGenerator)); + vectorField.Gradient = vtkm::make_Vec( + vtkm::make_Vec(randomDist(g_RandomGenerator), + randomDist(g_RandomGenerator), + randomDist(g_RandomGenerator)), + vtkm::make_Vec(randomDist(g_RandomGenerator), + randomDist(g_RandomGenerator), + randomDist(g_RandomGenerator)), + vtkm::make_Vec(randomDist(g_RandomGenerator), + randomDist(g_RandomGenerator), + randomDist(g_RandomGenerator))); + vtkm::Vec,3> expectedVectorGradient = + vectorField.Gradient; + + TestDerivativeFunctor > testFunctorVector; + std::cout << "======== Uniform Point Coordinates 3D =====" << std::endl; + testFunctorVector.DoTestWithWCoords( + vtkm::CellShapeTagHexahedron(), + vtkm::VecRectilinearPointCoordinates<3>(origin, spacing), + vectorField, + expectedVectorGradient); + std::cout << "======== Uniform Point Coordinates 2D =====" << std::endl; + expectedVectorGradient[2] = vtkm::Vec(0.0); + testFunctorVector.DoTestWithWCoords( + vtkm::CellShapeTagQuad(), + vtkm::VecRectilinearPointCoordinates<2>(origin, spacing), + vectorField, + expectedVectorGradient); + std::cout << "======== Uniform Point Coordinates 1D =====" << std::endl; + expectedVectorGradient[1] = vtkm::Vec(0.0); + testFunctorVector.DoTestWithWCoords( + vtkm::CellShapeTagLine(), + vtkm::VecRectilinearPointCoordinates<1>(origin, spacing), + vectorField, + expectedVectorGradient); } } // anonymous namespace