Skip to content

vtkm::exec::JacobianFor3DCell possible bugs

Currently our dashboards that test with double as the default floating point type have raised some warnings with JacobianFor3DCell.

The warnings listed below are with clang but reproduce with better messaging the issue seen at https://open.cdash.org/viewBuildError.php?type=1&buildid=5351890:

vtkm/Types.h:631:44: warning: implicit conversion loses floating-point precision: 'const ComponentType' (aka 'const double') to 'const float' [-Wconversion]
       this->Components[i] = static_cast<T>(src[i]);
                             ~~~~~~~~~~~    ^~~~~~
vtkm/Types.h:894:7: note: in instantiation of function template specialization 'vtkm::detail::VecBase<vtkm::Vec<float, 3>, 3, vtkm::Vec<vtkm::Vec<float, 3>, 3> >::VecBase<double, vtkm::Vec<double, 3> >' requested here
     : Superclass(src)
       ^
vtkm/exec/Jacobian.h:135:30: note: in instantiation of function template specialization 'vtkm::Vec<vtkm::Vec<float, 3>, 3>::Vec<double>' requested here
   vtkm::Vec<JacobianType, 3> pc(pcoords);
                              ^
vtkm/exec/CellDerivative.h:200:15: note: in instantiation of function template specialization 'vtkm::exec::JacobianFor3DCell<vtkm::Vec<vtkm::Vec<double, 3>, 8>, double, vtkm::Vec<float, 3> >' requested here
   vtkm::exec::JacobianFor3DCell(wCoords, pcoords, jacobianTranspose, CellShapeTag());

In a first pass I didn't understand why this was causing a warning. An explicit cast of double to float should not cause a conversion warning. The secret is realizing that while src is a vtkm::Vec<double,3> the Components are a vtkm::Vec<vtkm::Vec<float, 3>, 3>. So in effect the warning is that we are trying to cast a double to a vtkm::Vec<float, 3> and that is a problem!

Now into the actual bug. If we look at the calling site we see the code:

template <typename WorldCoordType, typename ParametricCoordType, typename JacobianType>
VTKM_EXEC void JacobianFor3DCell(const WorldCoordType& wCoords,
                                 const vtkm::Vec<ParametricCoordType, 3>& pcoords,
                                 vtkm::Matrix<JacobianType, 3, 3>& jacobian,
                                 vtkm::CellShapeTagHexahedron)
{
  vtkm::Vec<JacobianType, 3> pc(pcoords); //<----- this line!
  vtkm::Vec<JacobianType, 3> rc = vtkm::Vec<JacobianType, 3>(JacobianType(1)) - pc;

  jacobian = vtkm::Matrix<JacobianType, 3, 3>(JacobianType(0));
  VTKM_DERIVATIVE_WEIGHTS_HEXAHEDRON(pc, rc, VTKM_ACCUM_JACOBIAN_3D);
}

To my limited understanding the desire is to have pc be equal to pcoords when JacobianType are scalar values. When JacobianType is actually a vtkm::Vec I don't know what the desired behavior should be. Should each component of pc be equal to the contents of pcoords?

For example if pcoords is {0, 1, 2} and JacobainType is Vec<T,3> should pc look like:

{ {0, 1, 2}, {0, 1, 2}, {0, 1, 2} }

or should it look like:

{ {0, 0, 0}, {0, 0, 0}, {0, 0, 0} }

Currently we are getting the second behavior. Personally I believe people would expect the first behavior.

Edit:

Thinking about this issue more, I don't think the correct solution is at the Jacobian level, I think the issue is with vtkm::Vec and defining behavior when constructing a vtkm::Vec of vtkm::Vec.

Edited by Robert Maynard