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

Merge topic 'fix-special-jacobian-cases'

Browse Source 
43da67f7 Fix Jacobian cases for wedge and pyramid

Acked-by: Kitware Robot <kwrobot@kitware.com>
Acked-by: Sujin Philip <sujin.philip@kitware.com>
Merge-request: !1077
This commit is contained in:
Kenneth Moreland 2018-02-08 16:56:08 +00:00 committed by Kitware Robot
commit 6a908cc9da
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