diff --git a/vtkm/Math.h b/vtkm/Math.h index 9e14ba2fd..472939f3b 100644 --- a/vtkm/Math.h +++ b/vtkm/Math.h @@ -1550,6 +1550,35 @@ struct FloatLimits } }; +template +struct FloatLimits< vtkm::Vec > +{ + typedef vtkm::detail::IEEE754Bits32 BitsType; + + VTKM_EXEC_CONT_EXPORT + static vtkm::Vec Nan() { + BitsType nan = {VTKM_NAN_BITS_32}; + return vtkm::Vec(nan.scalar); + } + + VTKM_EXEC_CONT_EXPORT + static vtkm::Vec Infinity() { + BitsType inf = {VTKM_INF_BITS_32}; + return vtkm::Vec(inf.scalar); + } + + VTKM_EXEC_CONT_EXPORT + static vtkm::Vec NegativeInfinity() { + BitsType neginf = {VTKM_NEG_INF_BITS_32}; + return vtkm::Vec(neginf.scalar); + } + + VTKM_EXEC_CONT_EXPORT + static vtkm::Vec Epsilon() { + return vtkm::Vec(VTKM_EPSILON_32); + } +}; + template<> struct FloatLimits { @@ -1579,6 +1608,35 @@ struct FloatLimits } }; +template +struct FloatLimits< vtkm::Vec > +{ + typedef vtkm::detail::IEEE754Bits64 BitsType; + + VTKM_EXEC_CONT_EXPORT + static vtkm::Vec Nan() { + BitsType nan = {VTKM_NAN_BITS_64}; + return vtkm::Vec(nan.scalar); + } + + VTKM_EXEC_CONT_EXPORT + static vtkm::Vec Infinity() { + BitsType inf = {VTKM_INF_BITS_64}; + return vtkm::Vec(inf.scalar); + } + + VTKM_EXEC_CONT_EXPORT + static vtkm::Vec NegativeInfinity() { + BitsType neginf = {VTKM_NEG_INF_BITS_64}; + return vtkm::Vec(neginf.scalar); + } + + VTKM_EXEC_CONT_EXPORT + static vtkm::Vec Epsilon() { + return vtkm::Vec(VTKM_EPSILON_64); + } +}; + #undef VTKM_NAN_BITS_32 #undef VTKM_INF_BITS_32 #undef VTKM_NEG_INF_BITS_32 diff --git a/vtkm/Math.h.in b/vtkm/Math.h.in index e3f19e89d..63d0ad323 100644 --- a/vtkm/Math.h.in +++ b/vtkm/Math.h.in @@ -566,6 +566,35 @@ struct FloatLimits } }; +template +struct FloatLimits< vtkm::Vec > +{ + typedef vtkm::detail::IEEE754Bits32 BitsType; + + VTKM_EXEC_CONT_EXPORT + static vtkm::Vec Nan() { + BitsType nan = {VTKM_NAN_BITS_32}; + return vtkm::Vec(nan.scalar); + } + + VTKM_EXEC_CONT_EXPORT + static vtkm::Vec Infinity() { + BitsType inf = {VTKM_INF_BITS_32}; + return vtkm::Vec(inf.scalar); + } + + VTKM_EXEC_CONT_EXPORT + static vtkm::Vec NegativeInfinity() { + BitsType neginf = {VTKM_NEG_INF_BITS_32}; + return vtkm::Vec(neginf.scalar); + } + + VTKM_EXEC_CONT_EXPORT + static vtkm::Vec Epsilon() { + return vtkm::Vec(VTKM_EPSILON_32); + } +}; + template<> struct FloatLimits { @@ -595,6 +624,35 @@ struct FloatLimits } }; +template +struct FloatLimits< vtkm::Vec > +{ + typedef vtkm::detail::IEEE754Bits64 BitsType; + + VTKM_EXEC_CONT_EXPORT + static vtkm::Vec Nan() { + BitsType nan = {VTKM_NAN_BITS_64}; + return vtkm::Vec(nan.scalar); + } + + VTKM_EXEC_CONT_EXPORT + static vtkm::Vec Infinity() { + BitsType inf = {VTKM_INF_BITS_64}; + return vtkm::Vec(inf.scalar); + } + + VTKM_EXEC_CONT_EXPORT + static vtkm::Vec NegativeInfinity() { + BitsType neginf = {VTKM_NEG_INF_BITS_64}; + return vtkm::Vec(neginf.scalar); + } + + VTKM_EXEC_CONT_EXPORT + static vtkm::Vec Epsilon() { + return vtkm::Vec(VTKM_EPSILON_64); + } +}; + #undef VTKM_NAN_BITS_32 #undef VTKM_INF_BITS_32 #undef VTKM_NEG_INF_BITS_32 diff --git a/vtkm/Matrix.h b/vtkm/Matrix.h index 2a842cb04..80deaed6a 100644 --- a/vtkm/Matrix.h +++ b/vtkm/Matrix.h @@ -387,7 +387,7 @@ void MatrixLUPFactor(vtkm::Matrix &A, { permutation[index] = index; } - inversionParity = 1; + inversionParity = T(1); valid = true; for (vtkm::IdComponent rowIndex = 0; rowIndex < Size; rowIndex++) diff --git a/vtkm/Types.h b/vtkm/Types.h index aa7d23ce4..54de1fd7c 100644 --- a/vtkm/Types.h +++ b/vtkm/Types.h @@ -491,6 +491,16 @@ public: return result; } + VTKM_EXEC_CONT_EXPORT + DerivedClass& operator+=(const DerivedClass& other) + { + for (vtkm::IdComponent i = 0; i < Size; ++i) + { + this->Components[i] += other[i]; + } + return *static_cast(this); + } + VTKM_EXEC_CONT_EXPORT DerivedClass operator-(const DerivedClass& other) const { @@ -502,6 +512,16 @@ public: return result; } + VTKM_EXEC_CONT_EXPORT + DerivedClass& operator-=(const DerivedClass& other) + { + for (vtkm::IdComponent i = 0; i < Size; ++i) + { + this->Components[i] -= other[i]; + } + return *static_cast(this); + } + VTKM_EXEC_CONT_EXPORT DerivedClass operator*(const DerivedClass& other) const { @@ -513,6 +533,16 @@ public: return result; } + VTKM_EXEC_CONT_EXPORT + DerivedClass& operator*=(const DerivedClass& other) + { + for (vtkm::IdComponent i = 0; i < Size; ++i) + { + this->Components[i] *= other[i]; + } + return *static_cast(this); + } + VTKM_EXEC_CONT_EXPORT DerivedClass operator/(const DerivedClass& other) const { @@ -524,6 +554,16 @@ public: return result; } + VTKM_EXEC_CONT_EXPORT + DerivedClass& operator/=(const DerivedClass& other) + { + for (vtkm::IdComponent i = 0; i < Size; ++i) + { + this->Components[i] /= other[i]; + } + return *static_cast(this); + } + #if (defined(VTKM_GCC) || defined(VTKM_CLANG)) #pragma GCC diagnostic pop #endif // gcc || clang diff --git a/vtkm/exec/CellDerivative.h b/vtkm/exec/CellDerivative.h index 0dd511303..ebe2f8ff2 100644 --- a/vtkm/exec/CellDerivative.h +++ b/vtkm/exec/CellDerivative.h @@ -59,9 +59,9 @@ ParametricDerivative(const FieldVecType &field, typedef vtkm::Vec GradientType; GradientType pc(pcoords); - GradientType rc = GradientType(1) - pc; + GradientType rc = GradientType(FieldType(1)) - pc; - GradientType parametricDerivative(0); + GradientType parametricDerivative(FieldType(0)); VTKM_ACCUM_PARAMETRIC_DERIVATIVE_3D(0, -rc[1]*rc[2], -rc[0]*rc[2], -rc[0]*rc[1]); VTKM_ACCUM_PARAMETRIC_DERIVATIVE_3D(1, rc[1]*rc[2], -pc[0]*rc[2], -pc[0]*rc[1]); VTKM_ACCUM_PARAMETRIC_DERIVATIVE_3D(2, pc[1]*rc[2], pc[0]*rc[2], -pc[0]*pc[1]); @@ -146,9 +146,9 @@ ParametricDerivative(const FieldVecType &field, GradientType pc(static_cast(pcoords[0]), static_cast(pcoords[1])); - GradientType rc = GradientType(1) - pc; + GradientType rc = GradientType(FieldType(1)) - pc; - GradientType parametricDerivative(0); + GradientType parametricDerivative(FieldType(0)); VTKM_ACCUM_PARAMETRIC_DERIVATIVE_2D(0, -rc[1], -rc[0]); VTKM_ACCUM_PARAMETRIC_DERIVATIVE_2D(1, rc[1], -pc[0]); VTKM_ACCUM_PARAMETRIC_DERIVATIVE_2D(2, pc[1], pc[0]); @@ -375,7 +375,7 @@ CellDerivative(const FieldVecType &field, VTKM_ASSERT(wCoords.GetNumberOfComponents() == 1); typedef vtkm::Vec GradientType; - return GradientType(0,0,0); + return vtkm::TypeTraits::ZeroInitialization(); } //----------------------------------------------------------------------------- @@ -467,7 +467,7 @@ TriangleDerivative(const vtkm::Vec &field, vtkm::MatrixSetRow(A, 1, v1); vtkm::MatrixSetRow(A, 2, n); - GradientType b(field[1] - field[0], field[2] - field[0], 0); + 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 @@ -677,7 +677,7 @@ CellDerivative(const FieldVecType &field, VecT pc(static_cast(pcoords[0]), static_cast(pcoords[1])); - VecT rc = VecT(1) - pc; + VecT rc = VecT(T(1)) - pc; VecT sum = field[0]*VecT(-rc[1], -rc[0]); sum = sum + field[1]*VecT( rc[1], -pc[0]); @@ -686,7 +686,7 @@ CellDerivative(const FieldVecType &field, return vtkm::Vec(sum[0]/wCoords.GetSpacing()[0], sum[1]/wCoords.GetSpacing()[1], - 0); + T(0)); } //----------------------------------------------------------------------------- @@ -813,7 +813,7 @@ CellDerivative(const FieldVecType &field, VecT pc(static_cast(pcoords[0]), static_cast(pcoords[1]), static_cast(pcoords[2])); - VecT rc = VecT(1) - pc; + 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]); diff --git a/vtkm/exec/Jacobian.h b/vtkm/exec/Jacobian.h index 576ace5b1..d319a35be 100644 --- a/vtkm/exec/Jacobian.h +++ b/vtkm/exec/Jacobian.h @@ -194,9 +194,9 @@ void JacobianFor3DCell(const WorldCoordType &wCoords, vtkm::CellShapeTagHexahedron) { vtkm::Vec pc(pcoords); - vtkm::Vec rc = vtkm::Vec(1) - pc; + vtkm::Vec rc = vtkm::Vec(JacobianType(1)) - pc; - jacobian = vtkm::Matrix(0); + jacobian = vtkm::Matrix(JacobianType(0)); VTKM_DERIVATIVE_WEIGHTS_HEXAHEDRON(pc, rc, VTKM_ACCUM_JACOBIAN_3D); } @@ -272,10 +272,10 @@ void JacobianFor2DCell(const WorldCoordType &wCoords, { vtkm::Vec pc(static_cast(pcoords[0]), static_cast(pcoords[1])); - vtkm::Vec rc = vtkm::Vec(1) - pc; + vtkm::Vec rc = vtkm::Vec(JacobianType(1)) - pc; vtkm::Vec wcoords2d; - jacobian = vtkm::Matrix(0); + jacobian = vtkm::Matrix(JacobianType(0)); VTKM_DERIVATIVE_WEIGHTS_QUAD(pc, rc, VTKM_ACCUM_JACOBIAN_2D); } diff --git a/vtkm/worklet/Gradient.h b/vtkm/worklet/Gradient.h index ac931a029..b0fa0681b 100644 --- a/vtkm/worklet/Gradient.h +++ b/vtkm/worklet/Gradient.h @@ -35,14 +35,20 @@ namespace worklet { struct GradientOutTypes : vtkm::ListTagBase< vtkm::Vec, - vtkm::Vec + vtkm::Vec, + vtkm::Vec< vtkm::Vec, 3>, + vtkm::Vec< vtkm::Vec, 3>, + vtkm::Vec< vtkm::Vec, 3>, + vtkm::Vec< vtkm::Vec, 3>, + vtkm::Vec< vtkm::Vec, 3>, + vtkm::Vec< vtkm::Vec, 3> > { }; struct CellGradient : vtkm::worklet::WorkletMapPointToCell { typedef void ControlSignature(CellSetIn, FieldInPoint pointCoordinates, - FieldInPoint inputField, + FieldInPoint inputField, FieldOutCell outputField); typedef void ExecutionSignature(CellShape, PointCount, _2, _3, _4); @@ -64,8 +70,8 @@ struct CellGradient : vtkm::worklet::WorkletMapPointToCell using InValueType = typename FieldInVecType::ComponentType; using InDimensionTag = typename TypeTraits::DimensionalityTag; - //grad the dimension tag for the output component type - using OutValueType = typename VecTraits::ComponentType; + //grab the dimension tag for the output component type + using OutValueType = typename FieldOutType::ComponentType; using OutDimensionTag = typename TypeTraits::DimensionalityTag; //Verify that input and output dimension tags match @@ -99,6 +105,20 @@ struct CellGradient : vtkm::worklet::WorkletMapPointToCell std::false_type) const { //this is invalid + // std::cout << "calling invalid compute" << std::endl; + // using InValueType = typename FieldInVecType::ComponentType; + // using InDimensionTag = typename TypeTraits::DimensionalityTag; + + // using OutValueType = typename VecTraits::ComponentType; + // using OutDimensionTag = typename TypeTraits::DimensionalityTag; + + // std::cout << typeid(InDimensionTag).name() << '\n'; + // std::cout << typeid(FieldInVecType).name() << '\n'; + + // std::cout << typeid(OutDimensionTag).name() << '\n'; + // std::cout << typeid(FieldOutType).name() << '\n'; + + // std::cout << std::endl; } }; diff --git a/vtkm/worklet/testing/UnitTestCellGradient.cxx b/vtkm/worklet/testing/UnitTestCellGradient.cxx index 5d27cfe86..01de8de0e 100644 --- a/vtkm/worklet/testing/UnitTestCellGradient.cxx +++ b/vtkm/worklet/testing/UnitTestCellGradient.cxx @@ -78,6 +78,62 @@ void TestCellGradientUniform3D() test_equal(result.GetPortalConstControl().Get(i), expected[i]), "Wrong result for CellGradient worklet on 3D uniform data"); } + +} + +void TestCellGradientUniform3DWithVectorField() +{ + std::cout << "Testing CellGradient Worklet with a vector field on 3D strucutred data" << std::endl; + vtkm::cont::testing::MakeTestDataSet testDataSet; + vtkm::cont::DataSet dataSet = testDataSet.Make3DUniformDataSet0(); + + //Verify that we can compute the gradient of a 3 component vector + const int nVerts = 18; + vtkm::Float64 vars[nVerts] = {10.1, 20.1, 30.1, 40.1, 50.2, + 60.2, 70.2, 80.2, 90.3, 100.3, + 110.3, 120.3, 130.4, 140.4, + 150.4, 160.4, 170.5, 180.5}; + std::vector< vtkm::Vec > vec(18); + for(std::size_t i=0; i < vec.size(); ++i) + { + vec[i] = vtkm::make_Vec(vars[i],vars[i],vars[i]); + } + vtkm::cont::ArrayHandle< vtkm::Vec > input = + vtkm::cont::make_ArrayHandle(vec); + + //we need to add Vec3 array to the dataset + vtkm::cont::ArrayHandle< vtkm::Vec< vtkm::Vec, 3> > result; + vtkm::worklet::DispatcherMapTopology dispatcher; + dispatcher.Invoke(dataSet.GetCellSet(), + dataSet.GetCoordinateSystem(), + input, + result); + + vtkm::Vec< vtkm::Vec, 3> expected[4] = { + { {10.025,10.025,10.025}, {30.075,30.075,30.075}, {60.125,60.125,60.125} }, + { {10.025,10.025,10.025}, {30.075,30.075,30.075}, {60.125,60.125,60.125} }, + { {10.025,10.025,10.025}, {30.075,30.075,30.075}, {60.175,60.175,60.175} }, + { {10.025,10.025,10.025}, {30.075,30.075,30.075}, {60.175,60.175,60.175} } + }; + for (int i = 0; i < 4; ++i) + { + vtkm::Vec< vtkm::Vec, 3> e = expected[i]; + vtkm::Vec< vtkm::Vec, 3> r = result.GetPortalConstControl().Get(i); + + std::cout << "0: " << e[0] << " " << r[0] << std::endl; + std::cout << "1: " << e[1] << " " << r[1] << std::endl; + std::cout << "2: " << e[2] << " " << r[2] << std::endl; + VTKM_TEST_ASSERT( + test_equal(e[0],r[0]), + "Wrong result for vec field CellGradient worklet on 3D uniform data"); + VTKM_TEST_ASSERT( + test_equal(e[1],r[1]), + "Wrong result for vec field CellGradient worklet on 3D uniform data"); + VTKM_TEST_ASSERT( + test_equal(e[2],r[2]), + "Wrong result for vec field CellGradient worklet on 3D uniform data"); + } + } void TestCellGradientExplicit() @@ -109,6 +165,7 @@ void TestCellGradient() { TestCellGradientUniform2D(); TestCellGradientUniform3D(); + TestCellGradientUniform3DWithVectorField(); TestCellGradientExplicit(); }