Shuenhoy/vtk-m
1
0
Fork 0
mirror of https://gitlab.kitware.com/vtk/vtk-m synced 2024-09-08 13:23:51 +00:00
Code Issues 2 Actions Packages Projects Releases Wiki Activity

Fix Jacobian cases for wedge and pyramid

Browse Source 
When the implementation for the cell derivative functions was made, the
special cases for creating the Jacobian for wedge and pyramid cell
shapes was not working. Instead, it just used the hexahedra case for a
degenerate cell. This fixes the issues with the special cases.

There were multiple issues to be fixed. There were some complaints
by the compiler about types. There was a mistake in the pyramid table.
But the biggest issue was a problem with macro expansions. It was
the classic tale of forgetting that you have to encase parameters
in parenthesis to make sure that operator precedence comes out as
expected.
This commit is contained in:
Kenneth Moreland 2018-01-31 14:16:27 -07:00
parent bfc66729ff
commit 43da67f747
2 changed files with 42 additions and 143 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

48
vtkm/exec/CellDerivative.h
View File

@ -45,9 +45,9 @@ namespace exec
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
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].
@ -64,14 +64,7 @@ VTKM_EXEC vtkm::Vec<typename FieldVecType::ComponentType, 3> ParametricDerivativ
GradientType rc = GradientType(FieldType(1)) - pc;
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]);
VTKM_ACCUM_PARAMETRIC_DERIVATIVE_3D(3, -pc[1] * rc[2], rc[0] * rc[2], -rc[0] * pc[1]);
VTKM_ACCUM_PARAMETRIC_DERIVATIVE_3D(4, -rc[1] * pc[2], -rc[0] * pc[2], rc[0] * rc[1]);
VTKM_ACCUM_PARAMETRIC_DERIVATIVE_3D(5, rc[1] * pc[2], -pc[0] * pc[2], pc[0] * rc[1]);
VTKM_ACCUM_PARAMETRIC_DERIVATIVE_3D(6, pc[1] * pc[2], pc[0] * pc[2], pc[0] * pc[1]);
VTKM_ACCUM_PARAMETRIC_DERIVATIVE_3D(7, -pc[1] * pc[2], rc[0] * pc[2], rc[0] * pc[1]);
VTKM_DERIVATIVE_WEIGHTS_HEXAHEDRON(pc, rc, VTKM_ACCUM_PARAMETRIC_DERIVATIVE_3D);
return parametricDerivative;
}
@ -81,22 +74,16 @@ VTKM_EXEC vtkm::Vec<typename FieldVecType::ComponentType, 3> ParametricDerivativ
const vtkm::Vec<ParametricCoordType, 3>& pcoords,
vtkm::CellShapeTagWedge)
{
#if 0
// This is not working. Just leverage the hexahedron code that is working.
using FieldType = typename FieldVecType::ComponentType;
using GradientType = vtkm::Vec<FieldType,3>;
using GradientType = vtkm::Vec<FieldType, 3>;
GradientType pc(pcoords);
GradientType rc = GradientType(1) - pc;
GradientType rc = GradientType(FieldType(1)) - pc;
GradientType parametricDerivative(0);
GradientType parametricDerivative(FieldType(0));
VTKM_DERIVATIVE_WEIGHTS_WEDGE(pc, rc, VTKM_ACCUM_PARAMETRIC_DERIVATIVE_3D);
return parametricDerivative;
#else
return ParametricDerivative(
vtkm::exec::internal::PermuteWedgeToHex(field), pcoords, vtkm::CellShapeTagHexahedron());
#endif
}
template <typename FieldVecType, typename ParametricCoordType>
@ -105,29 +92,23 @@ VTKM_EXEC vtkm::Vec<typename FieldVecType::ComponentType, 3> ParametricDerivativ
const vtkm::Vec<ParametricCoordType, 3>& pcoords,
vtkm::CellShapeTagPyramid)
{
#if 0
// This is not working. Just leverage the hexahedron code that is working.
using FieldType = typename FieldVecType::ComponentType;
using GradientType = vtkm::Vec<FieldType,3>;
using GradientType = vtkm::Vec<FieldType, 3>;
GradientType pc(pcoords);
GradientType rc = GradientType(1) - pc;
GradientType rc = GradientType(FieldType(1)) - pc;
GradientType parametricDerivative(0);
GradientType parametricDerivative(FieldType(0));
VTKM_DERIVATIVE_WEIGHTS_PYRAMID(pc, rc, VTKM_ACCUM_PARAMETRIC_DERIVATIVE_3D);
return parametricDerivative;
#else
return ParametricDerivative(
vtkm::exec::internal::PermutePyramidToHex(field), pcoords, vtkm::CellShapeTagHexahedron());
#endif
}
#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
parametricDerivative[0] += field[pointIndex] * (weight0); \
parametricDerivative[1] += field[pointIndex] * (weight1)
template <typename FieldVecType, typename ParametricCoordType>
VTKM_EXEC vtkm::Vec<typename FieldVecType::ComponentType, 2> ParametricDerivative(
@ -142,10 +123,7 @@ VTKM_EXEC vtkm::Vec<typename FieldVecType::ComponentType, 2> ParametricDerivativ
GradientType rc = GradientType(FieldType(1)) - pc;
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]);
VTKM_ACCUM_PARAMETRIC_DERIVATIVE_2D(3, -pc[1], rc[0]);
VTKM_DERIVATIVE_WEIGHTS_QUAD(pc, rc, VTKM_ACCUM_PARAMETRIC_DERIVATIVE_2D);
return parametricDerivative;
}

137
vtkm/exec/Jacobian.h
View File

@ -68,98 +68,38 @@ struct Space2D
}
};
// Given a series of point values for a wedge, return a new series of point
// for a hexahedron that has the same interpolation within the wedge.
template <typename FieldVecType>
VTKM_EXEC vtkm::Vec<typename FieldVecType::ComponentType, 8> PermuteWedgeToHex(
const FieldVecType& field)
{
vtkm::Vec<typename FieldVecType::ComponentType, 8> hexField;
hexField[0] = field[0];
hexField[1] = field[2];
hexField[2] = field[2] + field[1] - field[0];
hexField[3] = field[1];
hexField[4] = field[3];
hexField[5] = field[5];
hexField[6] = field[5] + field[4] - field[3];
hexField[7] = field[4];
return hexField;
}
// Given a series of point values for a pyramid, return a new series of point
// for a hexahedron that has the same interpolation within the pyramid.
template <typename FieldVecType>
VTKM_EXEC vtkm::Vec<typename FieldVecType::ComponentType, 8> PermutePyramidToHex(
const FieldVecType& field)
{
using T = typename FieldVecType::ComponentType;
vtkm::Vec<T, 8> hexField;
T baseCenter = T(0.25f) * (field[0] + field[1] + field[2] + field[3]);
hexField[0] = field[0];
hexField[1] = field[1];
hexField[2] = field[2];
hexField[3] = field[3];
hexField[4] = field[4] + (field[0] - baseCenter);
hexField[5] = field[4] + (field[1] - baseCenter);
hexField[6] = field[4] + (field[2] - baseCenter);
hexField[7] = field[4] + (field[3] - baseCenter);
return hexField;
}
} //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_VOXEL(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], rc[0] * rc[2], -rc[0] * pc[1]); \
call(3, pc[1] * rc[2], pc[0] * rc[2], -pc[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], rc[0] * pc[2], rc[0] * pc[1]); \
call(7, pc[1] * pc[2], pc[0] * pc[2], pc[0] * pc[1])
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], -1.0f + pc[0] + pc[1]); \
call(1, 0.0f, rc[2], -pc[1]); \
call(2, rc[2], 0.0f, -pc[0]); \
call(3, -pc[2], -pc[2], 1.0f - pc[0] - pc[1]); \
call(4, 0.0f, pc[2], pc[1]); \
call(5, pc[2], 0.0f, pc[0])
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(3, 0.0f, 0.0f, 1.0f)
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])
#define VTKM_DERIVATIVE_WEIGHTS_PIXEL(pc, rc, call) \
call(0, -rc[1], -rc[0]); \
call(1, rc[1], -pc[0]); \
call(2, -pc[1], rc[0]); \
call(3, pc[1], pc[0])
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
@ -205,17 +145,11 @@ VTKM_EXEC void JacobianFor3DCell(const WorldCoordType& wCoords,
vtkm::Matrix<JacobianType, 3, 3>& jacobian,
vtkm::CellShapeTagWedge)
{
#if 0
// This is not working. Just leverage the hexahedron code that is working.
vtkm::Vec<JacobianType,3> pc(pcoords);
vtkm::Vec<JacobianType,3> rc = vtkm::Vec<JacobianType,3>(1) - pc;
vtkm::Vec<JacobianType, 3> pc(pcoords);
vtkm::Vec<JacobianType, 3> rc = vtkm::Vec<JacobianType, 3>(JacobianType(1)) - pc;
jacobian = vtkm::Matrix<JacobianType,3,3>(0);
jacobian = vtkm::Matrix<JacobianType, 3, 3>(JacobianType(0));
VTKM_DERIVATIVE_WEIGHTS_WEDGE(pc, rc, VTKM_ACCUM_JACOBIAN_3D);
#else
JacobianFor3DCell(
internal::PermuteWedgeToHex(wCoords), pcoords, jacobian, vtkm::CellShapeTagHexahedron());
#endif
}
template <typename WorldCoordType, typename ParametricCoordType, typename JacobianType>
@ -224,17 +158,11 @@ VTKM_EXEC void JacobianFor3DCell(const WorldCoordType& wCoords,
vtkm::Matrix<JacobianType, 3, 3>& jacobian,
vtkm::CellShapeTagPyramid)
{
#if 0
// This is not working. Just leverage the hexahedron code that is working.
vtkm::Vec<JacobianType,3> pc(pcoords);
vtkm::Vec<JacobianType,3> rc = vtkm::Vec<JacobianType,3>(1) - pc;
vtkm::Vec<JacobianType, 3> pc(pcoords);
vtkm::Vec<JacobianType, 3> rc = vtkm::Vec<JacobianType, 3>(JacobianType(1)) - pc;
jacobian = vtkm::Matrix<JacobianType,3,3>(0);
jacobian = vtkm::Matrix<JacobianType, 3, 3>(JacobianType(0));
VTKM_DERIVATIVE_WEIGHTS_PYRAMID(pc, rc, VTKM_ACCUM_JACOBIAN_3D);
#else
JacobianFor3DCell(
internal::PermutePyramidToHex(wCoords), pcoords, jacobian, vtkm::CellShapeTagHexahedron());
#endif
}
// Derivatives in quadrilaterals are computed in much the same way as
@ -316,13 +244,6 @@ void JacobianFor2DCell(const WorldCoordType &wCoords,
#undef VTKM_ACCUM_JACOBIAN_3D
#undef VTKM_ACCUM_JACOBIAN_2D
#undef VTKM_DERIVATIVE_WEIGHTS_HEXAHEDRON
#undef VTKM_DERIVATIVE_WEIGHTS_VOXEL
#undef VTKM_DERIVATIVE_WEIGHTS_WEDGE
#undef VTKM_DERIVATIVE_WEIGHTS_PYRAMID
#undef VTKM_DERIVATIVE_WEIGHTS_QUAD
#undef VTKM_DERIVATIVE_WEIGHTS_PIXEL
}
} // namespace vtkm::exec
#endif //vtk_m_exec_Jacobian_h