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

Use Light-Weight Cell Library

Browse Source
This commit is contained in:
Sujin Philip 2019-09-24 21:22:10 -04:00
parent e56b34fc98
commit e08d862f94
12 changed files with 452 additions and 2589 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
vtkm/CMakeLists.txt
View File

@ -76,6 +76,8 @@ if(VTKm_ENABLE_LOGGING)
endif()
add_subdirectory(thirdparty/optionparser)
add_subdirectory(thirdparty/taotuple)
add_subdirectory(thirdparty/vtkc)
add_subdirectory(testing)
add_subdirectory(internal)

69
vtkm/CellShape.h
View File

@ -13,6 +13,17 @@
#include <vtkm/StaticAssert.h>
#include <vtkm/Types.h>
#include <vtkc/Polygon.h>
#include <vtkc/Shapes.h>
// Vtk-c does not have tags for Empty and PolyLine. Define dummy tags here to
// avoid compilation errors. These tags are not used anywhere.
namespace vtkc
{
struct Empty;
struct PolyLine;
}
namespace vtkm
{
@ -22,21 +33,21 @@ namespace vtkm
enum CellShapeIdEnum
{
// Linear cells
CELL_SHAPE_EMPTY = 0,
CELL_SHAPE_VERTEX = 1,
CELL_SHAPE_EMPTY = vtkc::ShapeId::EMPTY,
CELL_SHAPE_VERTEX = vtkc::ShapeId::VERTEX,
//CELL_SHAPE_POLY_VERTEX = 2,
CELL_SHAPE_LINE = 3,
CELL_SHAPE_LINE = vtkc::ShapeId::LINE,
CELL_SHAPE_POLY_LINE = 4,
CELL_SHAPE_TRIANGLE = 5,
CELL_SHAPE_TRIANGLE = vtkc::ShapeId::TRIANGLE,
//CELL_SHAPE_TRIANGLE_STRIP = 6,
CELL_SHAPE_POLYGON = 7,
CELL_SHAPE_POLYGON = vtkc::ShapeId::POLYGON,
//CELL_SHAPE_PIXEL = 8,
CELL_SHAPE_QUAD = 9,
CELL_SHAPE_TETRA = 10,
CELL_SHAPE_QUAD = vtkc::ShapeId::QUAD,
CELL_SHAPE_TETRA = vtkc::ShapeId::TETRA,
//CELL_SHAPE_VOXEL = 11,
CELL_SHAPE_HEXAHEDRON = 12,
CELL_SHAPE_WEDGE = 13,
CELL_SHAPE_PYRAMID = 14,
CELL_SHAPE_HEXAHEDRON = vtkc::ShapeId::HEXAHEDRON,
CELL_SHAPE_WEDGE = vtkc::ShapeId::WEDGE,
CELL_SHAPE_PYRAMID = vtkc::ShapeId::PYRAMID,
NUMBER_OF_CELL_SHAPES
};
@ -58,6 +69,10 @@ struct CellShapeTagCheck : std::false_type
{
};
/// Convert VTK-m tag to VTK-c tag
template <typename VtkmCellShapeTag>
struct CellShapeTagVtkmToVtkc;
} // namespace internal
/// Checks that the argument is a proper cell shape tag. This is a handy
@ -94,6 +109,11 @@ struct CellShapeIdToTag
struct CellShapeTagCheck<vtkm::CellShapeTag##name> : std::true_type \
{ \
}; \
template <> \
struct CellShapeTagVtkmToVtkc<vtkm::CellShapeTag##name> \
{ \
using Type = vtkc::name; \
}; \
} \
static inline VTKM_EXEC_CONT const char* GetCellShapeName(vtkm::CellShapeTag##name) \
{ \
@ -140,6 +160,35 @@ struct CellShapeTagGeneric
vtkm::UInt8 Id;
};
namespace internal
{
template <typename VtkmCellShapeTag>
VTKM_EXEC_CONT inline typename CellShapeTagVtkmToVtkc<VtkmCellShapeTag>::Type make_VtkcCellShapeTag(
const VtkmCellShapeTag&,
vtkm::IdComponent numPoints = 0)
{
using VtkcCellShapeTag = typename CellShapeTagVtkmToVtkc<VtkmCellShapeTag>::Type;
static_cast<void>(numPoints); // unused
return VtkcCellShapeTag{};
}
VTKM_EXEC_CONT
inline vtkc::Polygon make_VtkcCellShapeTag(const vtkm::CellShapeTagPolygon&,
vtkm::IdComponent numPoints = 0)
{
return vtkc::Polygon(numPoints);
}
VTKM_EXEC_CONT
inline vtkc::Cell make_VtkcCellShapeTag(const vtkm::CellShapeTagGeneric& tag,
vtkm::IdComponent numPoints = 0)
{
return vtkc::Cell(static_cast<std::int8_t>(tag.Id), numPoints);
}
} // namespace internal
#define vtkmGenericCellShapeMacroCase(cellShapeId, call) \
case vtkm::cellShapeId: \
{ \

2
vtkm/cont/CMakeLists.txt
View File

@ -211,7 +211,7 @@ if(TARGET vtkm::openmp)
endif()
target_link_libraries(vtkm_cont PUBLIC vtkm_compiler_flags ${backends})
target_link_libraries(vtkm_cont PUBLIC vtkm_taotuple vtkm_optionparser vtkm_diy)
target_link_libraries(vtkm_cont PUBLIC vtkm_taotuple vtkm_optionparser vtkm_diy vtkm_vtkc)
if(TARGET vtkm_loguru)
target_link_libraries(vtkm_cont PRIVATE vtkm_loguru)
endif()

1
vtkm/exec/CMakeLists.txt
View File

@ -29,7 +29,6 @@ set(headers
ExecutionWholeArray.h
FieldNeighborhood.h
FunctorBase.h
Jacobian.h
ParametricCoordinates.h
PointLocator.h
PointLocatorUniformGrid.h

1130
vtkm/exec/CellDerivative.h
View File

File diff suppressed because it is too large Load Diff

70
vtkm/exec/CellInside.h
View File

@ -13,84 +13,32 @@
#include <vtkm/CellShape.h>
#include <vtkm/Types.h>
#include <vtkc/vtkc.h>
namespace vtkm
{
namespace exec
{
template <typename T, typename CellShapeTag>
static inline VTKM_EXEC bool CellInside(const vtkm::Vec<T, 3>& pcoords, CellShapeTag)
{
using VtkcTagType = typename vtkm::internal::CellShapeTagVtkmToVtkc<CellShapeTag>::Type;
return vtkc::cellInside(VtkcTagType{}, pcoords);
}
template <typename T>
static inline VTKM_EXEC bool CellInside(const vtkm::Vec<T, 3>&, vtkm::CellShapeTagEmpty)
{
return false;
}
template <typename T>
static inline VTKM_EXEC bool CellInside(const vtkm::Vec<T, 3>& pcoords, vtkm::CellShapeTagVertex)
{
return pcoords[0] == T(0) && pcoords[1] == T(0) && pcoords[2] == T(0);
}
template <typename T>
static inline VTKM_EXEC bool CellInside(const vtkm::Vec<T, 3>& pcoords, vtkm::CellShapeTagLine)
{
return pcoords[0] >= T(0) && pcoords[0] <= T(1);
}
template <typename T>
static inline VTKM_EXEC bool CellInside(const vtkm::Vec<T, 3>& pcoords, vtkm::CellShapeTagPolyLine)
{
return pcoords[0] >= T(0) && pcoords[0] <= T(1);
}
template <typename T>
static inline VTKM_EXEC bool CellInside(const vtkm::Vec<T, 3>& pcoords, vtkm::CellShapeTagTriangle)
{
return pcoords[0] >= T(0) && pcoords[1] >= T(0) && (pcoords[0] + pcoords[1] <= T(1));
}
template <typename T>
static inline VTKM_EXEC bool CellInside(const vtkm::Vec<T, 3>& pcoords, vtkm::CellShapeTagPolygon)
{
return ((pcoords[0] - T(0.5)) * (pcoords[0] - T(0.5))) +
((pcoords[1] - T(0.5)) * (pcoords[1] - T(0.5))) <=
T(0.25);
}
template <typename T>
static inline VTKM_EXEC bool CellInside(const vtkm::Vec<T, 3>& pcoords, vtkm::CellShapeTagQuad)
{
return pcoords[0] >= T(0) && pcoords[0] <= T(1) && pcoords[1] >= T(0) && pcoords[1] <= T(1);
}
template <typename T>
static inline VTKM_EXEC bool CellInside(const vtkm::Vec<T, 3>& pcoords, vtkm::CellShapeTagTetra)
{
return pcoords[0] >= T(0) && pcoords[1] >= T(0) && pcoords[2] >= T(0) &&
(pcoords[0] + pcoords[1] + pcoords[2] <= T(1));
}
template <typename T>
static inline VTKM_EXEC bool CellInside(const vtkm::Vec<T, 3>& pcoords,
vtkm::CellShapeTagHexahedron)
{
return pcoords[0] >= T(0) && pcoords[0] <= T(1) && pcoords[1] >= T(0) && pcoords[1] <= T(1) &&
pcoords[2] >= T(0) && pcoords[2] <= T(1);
}
template <typename T>
static inline VTKM_EXEC bool CellInside(const vtkm::Vec<T, 3>& pcoords, vtkm::CellShapeTagWedge)
{
return pcoords[0] >= T(0) && pcoords[1] >= T(0) && pcoords[2] >= T(0) && pcoords[2] <= T(1) &&
(pcoords[0] + pcoords[1] <= T(1));
}
template <typename T>
static inline VTKM_EXEC bool CellInside(const vtkm::Vec<T, 3>& pcoords, vtkm::CellShapeTagPyramid)
{
return pcoords[0] >= T(0) && pcoords[0] <= T(1) && pcoords[1] >= T(0) && pcoords[1] <= T(1) &&
pcoords[2] >= T(0) && pcoords[2] <= T(1);
}
/// Checks if the parametric coordinates `pcoords` are on the inside for the
/// specified cell type.
///

402
vtkm/exec/CellInterpolate.h
View File

@ -12,11 +12,11 @@
#include <vtkm/Assert.h>
#include <vtkm/CellShape.h>
#include <vtkm/Math.h>
#include <vtkm/VecAxisAlignedPointCoordinates.h>
#include <vtkm/VectorAnalysis.h>
#include <vtkm/exec/FunctorBase.h>
#include <vtkc/vtkc.h>
#if (defined(VTKM_GCC) || defined(VTKM_CLANG))
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wconversion"
@ -30,97 +30,30 @@ namespace exec
namespace internal
{
// This is really the WorldCoorindatesToParametericCoordinates function, but
// moved to this header file because it is required to interpolate in a
// polygon, which is divided into triangles.
template <typename WorldCoordVector>
VTKM_EXEC typename WorldCoordVector::ComponentType ReverseInterpolateTriangle(
const WorldCoordVector& pointWCoords,
const typename WorldCoordVector::ComponentType& wcoords)
template <typename VtkcCellShapeTag, typename FieldVecType, typename ParametricCoordType>
VTKM_EXEC typename FieldVecType::ComponentType CellInterpolateImpl(
VtkcCellShapeTag tag,
const FieldVecType& field,
const ParametricCoordType& pcoords,
const vtkm::exec::FunctorBase& worklet)
{
VTKM_ASSERT(pointWCoords.GetNumberOfComponents() == 3);
VTKM_ASSERT(tag.numberOfPoints() == field.GetNumberOfComponents());
// We will solve the world to parametric coordinates problem geometrically.
// Consider the parallelogram formed by wcoords and p0 of the triangle and
// the two adjacent edges. This parallelogram is equivalent to the
// axis-aligned rectangle anchored at the origin of parametric space.
//
// p2 |\ (1,0) |\ //
// | \ | \ //
// | \ | \ //
// | \ | \ //
// | \ | \ //
// | \ | (u,v) \ //
// | --- \ |-------* \ //
// | ---*wcoords | | \ //
// | | \ | | \ //
// p0 *--- | \ (0,0) *------------------\ (1,0) //
// ---| \ //
// x-- \ //
// --- \ //
// ---\ p1 //
//
// In this diagram, the distance between p0 and the point marked x divided by
// the length of the edge it is on is equal, by proportionality, to the u
// parametric coordinate. (The v coordinate follows the other edge
// accordingly.) Thus, if we can find the intersection at x (or more
// specifically the distance between p0 and x), then we can find that
// parametric coordinate.
//
// Because the triangle is in 3-space, we are actually going to intersect the
// edge with a plane that is parallel to the opposite edge of p0 and
// perpendicular to the triangle. This is partially because it is easy to
// find the intersection between a plane and a line and partially because the
// computation will work for points not on the plane. (The result is
// equivalent to a point projected on the plane.)
//
// First, we define an implicit plane as:
//
// Dot((p - wcoords), planeNormal) = 0
//
// where planeNormal is the normal to the plane (easily computed from the
// triangle), and p is any point in the plane. Next, we define the parametric
// form of the line:
//
// p(d) = (p1 - p0)d + p0
//
// Where d is the fraction of distance from p0 toward p1. Note that d is
// actually equal to the parametric coordinate we are trying to find. Once we
// compute it, we are done. We can skip the part about finding the actual
// coordinates of the intersection.
//
// Solving for the intersection is as simple as substituting the line's
// definition of p(d) into p for the plane equation. With some basic algebra
// you get:
//
// d = Dot((wcoords - p0), planeNormal)/Dot((p1-p0), planeNormal)
//
// From here, the u coordinate is simply d. The v coordinate follows
// similarly.
//
using Vector3 = typename WorldCoordVector::ComponentType;
using T = typename Vector3::ComponentType;
Vector3 pcoords(T(0));
Vector3 triangleNormal = vtkm::TriangleNormal(pointWCoords[0], pointWCoords[1], pointWCoords[2]);
for (vtkm::IdComponent dimension = 0; dimension < 2; dimension++)
using FieldValueType = typename FieldVecType::ComponentType;
IdComponent numComponents = vtkm::VecTraits<FieldValueType>::GetNumberOfComponents(field[0]);
FieldValueType result(0);
auto status =
vtkc::interpolate(tag, vtkc::makeFieldAccessorNestedSOA(field, numComponents), pcoords, result);
if (status != vtkc::ErrorCode::SUCCESS)
{
Vector3 p0 = pointWCoords[0];
Vector3 p1 = pointWCoords[dimension + 1];
Vector3 p2 = pointWCoords[2 - dimension];
Vector3 planeNormal = vtkm::Cross(triangleNormal, p2 - p0);
T d = vtkm::Dot(wcoords - p0, planeNormal) / vtkm::Dot(p1 - p0, planeNormal);
pcoords[dimension] = d;
worklet.RaiseError(vtkc::errorString(status));
}
return pcoords;
}
return result;
}
} // namespace internal
//-----------------------------------------------------------------------------
/// \brief Interpolate a point field in a cell.
///
/// Given the point field values for each node and the parametric coordinates
@ -145,6 +78,19 @@ VTKM_EXEC typename FieldVecType::ComponentType CellInterpolate(
return result;
}
//-----------------------------------------------------------------------------
template <typename FieldVecType, typename ParametricCoordType, typename CellShapeTag>
VTKM_EXEC typename FieldVecType::ComponentType CellInterpolate(
const FieldVecType& pointFieldValues,
const vtkm::Vec<ParametricCoordType, 3>& pcoords,
CellShapeTag tag,
const vtkm::exec::FunctorBase& worklet)
{
auto vtkcTag =
vtkm::internal::make_VtkcCellShapeTag(tag, pointFieldValues.GetNumberOfComponents());
return internal::CellInterpolateImpl(vtkcTag, pointFieldValues, pcoords, worklet);
}
//-----------------------------------------------------------------------------
template <typename FieldVecType, typename ParametricCoordType>
VTKM_EXEC typename FieldVecType::ComponentType CellInterpolate(
@ -157,45 +103,6 @@ VTKM_EXEC typename FieldVecType::ComponentType CellInterpolate(
return typename FieldVecType::ComponentType();
}
//-----------------------------------------------------------------------------
template <typename FieldVecType, typename ParametricCoordType>
VTKM_EXEC typename FieldVecType::ComponentType CellInterpolate(
const FieldVecType& pointFieldValues,
const vtkm::Vec<ParametricCoordType, 3>,
vtkm::CellShapeTagVertex,
const vtkm::exec::FunctorBase& vtkmNotUsed(worklet))
{
VTKM_ASSERT(pointFieldValues.GetNumberOfComponents() == 1);
return pointFieldValues[0];
}
//-----------------------------------------------------------------------------
template <typename FieldVecType, typename ParametricCoordType>
VTKM_EXEC typename FieldVecType::ComponentType CellInterpolate(
const FieldVecType& pointFieldValues,
const vtkm::Vec<ParametricCoordType, 3>& parametricCoords,
vtkm::CellShapeTagLine,
const vtkm::exec::FunctorBase& vtkmNotUsed(worklet))
{
VTKM_ASSERT(pointFieldValues.GetNumberOfComponents() == 2);
return vtkm::Lerp(pointFieldValues[0], pointFieldValues[1], parametricCoords[0]);
}
template <typename ParametricCoordType>
VTKM_EXEC vtkm::Vec3f CellInterpolate(const vtkm::VecAxisAlignedPointCoordinates<1>& field,
const vtkm::Vec<ParametricCoordType, 3>& pcoords,
vtkm::CellShapeTagLine,
const vtkm::exec::FunctorBase&)
{
using T = vtkm::Vec3f;
const T& origin = field.GetOrigin();
const T& spacing = field.GetSpacing();
return T(
origin[0] + static_cast<vtkm::FloatDefault>(pcoords[0]) * spacing[0], origin[1], origin[2]);
}
//-----------------------------------------------------------------------------
template <typename FieldVecType, typename ParametricCoordType>
VTKM_EXEC typename FieldVecType::ComponentType CellInterpolate(
@ -207,12 +114,9 @@ VTKM_EXEC typename FieldVecType::ComponentType CellInterpolate(
const vtkm::IdComponent numPoints = field.GetNumberOfComponents();
VTKM_ASSERT(numPoints >= 1);
switch (numPoints)
if (numPoints == 1)
{
case 1:
return CellInterpolate(field, pcoords, vtkm::CellShapeTagVertex(), worklet);
case 2:
return CellInterpolate(field, pcoords, vtkm::CellShapeTagLine(), worklet);
return CellInterpolate(field, pcoords, vtkm::CellShapeTagVertex(), worklet);
}
using T = ParametricCoordType;
@ -220,50 +124,13 @@ VTKM_EXEC typename FieldVecType::ComponentType CellInterpolate(
T dt = 1 / static_cast<T>(numPoints - 1);
vtkm::IdComponent idx = static_cast<vtkm::IdComponent>(pcoords[0] / dt);
if (idx == numPoints - 1)
return field[numPoints - 1];
T t = (pcoords[0] - static_cast<T>(idx) * dt) / dt;
return vtkm::Lerp(field[idx], field[idx + 1], t);
}
template <typename ParametricCoordType>
VTKM_EXEC vtkm::Vec3f CellInterpolate(const vtkm::VecAxisAlignedPointCoordinates<1>& field,
const vtkm::Vec<ParametricCoordType, 3>& pcoords,
vtkm::CellShapeTagPolyLine,
const vtkm::exec::FunctorBase& worklet)
{
const vtkm::IdComponent numPoints = field.GetNumberOfComponents();
VTKM_ASSERT(numPoints >= 1);
switch (numPoints)
{
case 1:
return CellInterpolate(field, pcoords, vtkm::CellShapeTagVertex(), worklet);
case 2:
return CellInterpolate(field, pcoords, vtkm::CellShapeTagLine(), worklet);
return field[numPoints - 1];
}
using T = vtkm::Vec3f;
const T& origin = field.GetOrigin();
const T& spacing = field.GetSpacing();
return T(
origin[0] + static_cast<vtkm::FloatDefault>(pcoords[0]) * spacing[0], origin[1], origin[2]);
}
//-----------------------------------------------------------------------------
template <typename FieldVecType, typename ParametricCoordType>
VTKM_EXEC typename FieldVecType::ComponentType CellInterpolate(
const FieldVecType& field,
const vtkm::Vec<ParametricCoordType, 3>& pcoords,
vtkm::CellShapeTagTriangle,
const vtkm::exec::FunctorBase& vtkmNotUsed(worklet))
{
VTKM_ASSERT(field.GetNumberOfComponents() == 3);
using T = typename FieldVecType::ComponentType;
return static_cast<T>((field[0] * (1 - pcoords[0] - pcoords[1])) + (field[1] * pcoords[0]) +
(field[2] * pcoords[1]));
T pc = (pcoords[0] - static_cast<T>(idx) * dt) / dt;
return internal::CellInterpolateImpl(
vtkc::Line{}, vtkm::make_Vec(field[idx], field[idx + 1]), &pc, worklet);
}
//-----------------------------------------------------------------------------
@ -282,208 +149,29 @@ VTKM_EXEC typename FieldVecType::ComponentType CellInterpolate(
return CellInterpolate(field, pcoords, vtkm::CellShapeTagVertex(), worklet);
case 2:
return CellInterpolate(field, pcoords, vtkm::CellShapeTagLine(), worklet);
case 3:
return CellInterpolate(field, pcoords, vtkm::CellShapeTagTriangle(), worklet);
case 4:
return CellInterpolate(field, pcoords, vtkm::CellShapeTagQuad(), worklet);
default:
return internal::CellInterpolateImpl(vtkc::Polygon(numPoints), field, pcoords, worklet);
}
// If we are here, then there are 5 or more points on this polygon.
// Arrange the points such that they are on the circle circumscribed in the
// unit square from 0 to 1. That is, the point are on the circle centered at
// coordinate 0.5,0.5 with radius 0.5. The polygon is divided into regions
// defined by they triangle fan formed by the points around the center. This
// is C0 continuous but not necessarily C1 continuous. It is also possible to
// have a non 1 to 1 mapping between parametric coordinates world coordinates
// if the polygon is not planar or convex.
using FieldType = typename FieldVecType::ComponentType;
// Find the interpolation for the center point.
FieldType fieldCenter = field[0];
for (vtkm::IdComponent pointIndex = 1; pointIndex < numPoints; pointIndex++)
{
fieldCenter = fieldCenter + field[pointIndex];
}
fieldCenter = fieldCenter * FieldType(1.0f / static_cast<float>(numPoints));
if ((vtkm::Abs(pcoords[0] - 0.5f) < 4 * vtkm::Epsilon<ParametricCoordType>()) &&
(vtkm::Abs(pcoords[1] - 0.5f) < 4 * vtkm::Epsilon<ParametricCoordType>()))
{
return fieldCenter;
}
ParametricCoordType angle = vtkm::ATan2(pcoords[1] - 0.5f, pcoords[0] - 0.5f);
if (angle < 0)
{
angle += 2 * vtkm::Pi<ParametricCoordType>();
}
const ParametricCoordType deltaAngle = 2 * vtkm::Pi<ParametricCoordType>() / numPoints;
vtkm::IdComponent firstPointIndex =
static_cast<vtkm::IdComponent>(vtkm::Floor(angle / deltaAngle));
vtkm::IdComponent secondPointIndex = firstPointIndex + 1;
if (secondPointIndex == numPoints)
{
secondPointIndex = 0;
}
// Transform pcoords for polygon into pcoords for triangle.
vtkm::Vec<vtkm::Vec<ParametricCoordType, 3>, 3> polygonCoords;
polygonCoords[0][0] = 0.5f;
polygonCoords[0][1] = 0.5f;
polygonCoords[0][2] = 0;
polygonCoords[1][0] =
0.5f * (vtkm::Cos(deltaAngle * static_cast<ParametricCoordType>(firstPointIndex)) + 1);
polygonCoords[1][1] =
0.5f * (vtkm::Sin(deltaAngle * static_cast<ParametricCoordType>(firstPointIndex)) + 1);
polygonCoords[1][2] = 0.0f;
polygonCoords[2][0] =
0.5f * (vtkm::Cos(deltaAngle * static_cast<ParametricCoordType>(secondPointIndex)) + 1);
polygonCoords[2][1] =
0.5f * (vtkm::Sin(deltaAngle * static_cast<ParametricCoordType>(secondPointIndex)) + 1);
polygonCoords[2][2] = 0.0f;
vtkm::Vec<ParametricCoordType, 3> trianglePCoords =
vtkm::exec::internal::ReverseInterpolateTriangle(polygonCoords, pcoords);
// Set up parameters for triangle that pcoords is in
vtkm::Vec<FieldType, 3> triangleField;
triangleField[0] = fieldCenter;
triangleField[1] = field[firstPointIndex];
triangleField[2] = field[secondPointIndex];
// Now use the triangle interpolate
return vtkm::exec::CellInterpolate(
triangleField, trianglePCoords, vtkm::CellShapeTagTriangle(), worklet);
}
//-----------------------------------------------------------------------------
template <typename FieldVecType, typename ParametricCoordType>
VTKM_EXEC typename FieldVecType::ComponentType CellInterpolate(
const FieldVecType& field,
const vtkm::Vec<ParametricCoordType, 3>& pcoords,
vtkm::CellShapeTagQuad,
const vtkm::exec::FunctorBase& vtkmNotUsed(worklet))
{
VTKM_ASSERT(field.GetNumberOfComponents() == 4);
using T = typename FieldVecType::ComponentType;
T bottomInterp = vtkm::Lerp(field[0], field[1], pcoords[0]);
T topInterp = vtkm::Lerp(field[3], field[2], pcoords[0]);
return vtkm::Lerp(bottomInterp, topInterp, pcoords[1]);
}
template <typename ParametricCoordType>
VTKM_EXEC vtkm::Vec3f CellInterpolate(const vtkm::VecAxisAlignedPointCoordinates<2>& field,
const vtkm::Vec<ParametricCoordType, 3>& pcoords,
vtkm::CellShapeTagQuad,
const vtkm::exec::FunctorBase&)
const vtkm::exec::FunctorBase& worklet)
{
using T = vtkm::Vec3f;
const T& origin = field.GetOrigin();
const T& spacing = field.GetSpacing();
return T(origin[0] + static_cast<vtkm::FloatDefault>(pcoords[0]) * spacing[0],
origin[1] + static_cast<vtkm::FloatDefault>(pcoords[1]) * spacing[1],
origin[2]);
return internal::CellInterpolateImpl(vtkc::Pixel{}, field, pcoords, worklet);
}
//-----------------------------------------------------------------------------
template <typename FieldVecType, typename ParametricCoordType>
VTKM_EXEC typename FieldVecType::ComponentType CellInterpolate(
const FieldVecType& field,
const vtkm::Vec<ParametricCoordType, 3>& pcoords,
vtkm::CellShapeTagTetra,
const vtkm::exec::FunctorBase& vtkmNotUsed(worklet))
{
VTKM_ASSERT(field.GetNumberOfComponents() == 4);
using T = typename FieldVecType::ComponentType;
return static_cast<T>((field[0] * (1 - pcoords[0] - pcoords[1] - pcoords[2])) +
(field[1] * pcoords[0]) + (field[2] * pcoords[1]) +
(field[3] * pcoords[2]));
}
//-----------------------------------------------------------------------------
template <typename FieldVecType, typename ParametricCoordType>
VTKM_EXEC typename FieldVecType::ComponentType CellInterpolate(
const FieldVecType& field,
const vtkm::Vec<ParametricCoordType, 3>& pcoords,
vtkm::CellShapeTagHexahedron,
const vtkm::exec::FunctorBase& vtkmNotUsed(worklet))
{
VTKM_ASSERT(field.GetNumberOfComponents() == 8);
using T = typename FieldVecType::ComponentType;
T bottomFrontInterp = vtkm::Lerp(field[0], field[1], pcoords[0]);
T bottomBackInterp = vtkm::Lerp(field[3], field[2], pcoords[0]);
T topFrontInterp = vtkm::Lerp(field[4], field[5], pcoords[0]);
T topBackInterp = vtkm::Lerp(field[7], field[6], pcoords[0]);
T bottomInterp = vtkm::Lerp(bottomFrontInterp, bottomBackInterp, pcoords[1]);
T topInterp = vtkm::Lerp(topFrontInterp, topBackInterp, pcoords[1]);
return vtkm::Lerp(bottomInterp, topInterp, pcoords[2]);
}
template <typename ParametricCoordType>
VTKM_EXEC vtkm::Vec3f CellInterpolate(const vtkm::VecAxisAlignedPointCoordinates<3>& field,
const vtkm::Vec<ParametricCoordType, 3>& pcoords,
vtkm::CellShapeTagHexahedron,
const vtkm::exec::FunctorBase&)
const vtkm::exec::FunctorBase& worklet)
{
vtkm::Vec3f pcoordsCast(static_cast<vtkm::FloatDefault>(pcoords[0]),
static_cast<vtkm::FloatDefault>(pcoords[1]),
static_cast<vtkm::FloatDefault>(pcoords[2]));
return field.GetOrigin() + pcoordsCast * field.GetSpacing();
}
//-----------------------------------------------------------------------------
template <typename FieldVecType, typename ParametricCoordType>
VTKM_EXEC typename FieldVecType::ComponentType CellInterpolate(
const FieldVecType& field,
const vtkm::Vec<ParametricCoordType, 3>& pcoords,
vtkm::CellShapeTagWedge,
const vtkm::exec::FunctorBase& vtkmNotUsed(worklet))
{
VTKM_ASSERT(field.GetNumberOfComponents() == 6);
using T = typename FieldVecType::ComponentType;
T bottomInterp = static_cast<T>((field[0] * (1 - pcoords[0] - pcoords[1])) +
(field[1] * pcoords[1]) + (field[2] * pcoords[0]));
T topInterp = static_cast<T>((field[3] * (1 - pcoords[0] - pcoords[1])) +
(field[4] * pcoords[1]) + (field[5] * pcoords[0]));
return vtkm::Lerp(bottomInterp, topInterp, pcoords[2]);
}
//-----------------------------------------------------------------------------
template <typename FieldVecType, typename ParametricCoordType>
VTKM_EXEC typename FieldVecType::ComponentType CellInterpolate(
const FieldVecType& field,
const vtkm::Vec<ParametricCoordType, 3>& pcoords,
vtkm::CellShapeTagPyramid,
const vtkm::exec::FunctorBase& vtkmNotUsed(worklet))
{
VTKM_ASSERT(field.GetNumberOfComponents() == 5);
using T = typename FieldVecType::ComponentType;
T frontInterp = vtkm::Lerp(field[0], field[1], pcoords[0]);
T backInterp = vtkm::Lerp(field[3], field[2], pcoords[0]);
T baseInterp = vtkm::Lerp(frontInterp, backInterp, pcoords[1]);
return vtkm::Lerp(baseInterp, field[4], pcoords[2]);
return internal::CellInterpolateImpl(vtkc::Voxel{}, field, pcoords, worklet);
}
}
} // namespace vtkm::exec

239
vtkm/exec/Jacobian.h
View File

@ -1,239 +0,0 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
//
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#ifndef vtk_m_exec_Jacobian_h
#define vtk_m_exec_Jacobian_h
#include <vtkm/Assert.h>
#include <vtkm/CellShape.h>
#include <vtkm/Math.h>
#include <vtkm/Matrix.h>
#include <vtkm/VectorAnalysis.h>
namespace vtkm
{
namespace exec
{
namespace internal
{
template <typename T>
struct Space2D
{
using Vec3 = vtkm::Vec<T, 3>;
using Vec2 = vtkm::Vec<T, 2>;
Vec3 Origin;
Vec3 Basis0;
Vec3 Basis1;
VTKM_EXEC
Space2D(const Vec3& origin, const Vec3& pointFirst, const Vec3& pointLast)
{
this->Origin = origin;
this->Basis0 = vtkm::Normal(pointFirst - this->Origin);
Vec3 n = vtkm::Cross(this->Basis0, pointLast - this->Origin);
this->Basis1 = vtkm::Normal(vtkm::Cross(this->Basis0, n));
}
VTKM_EXEC
Vec2 ConvertCoordToSpace(const Vec3 coord) const
{
Vec3 vec = coord - this->Origin;
return Vec2(vtkm::Dot(vec, this->Basis0), vtkm::Dot(vec, this->Basis1));
}
template <typename U>
VTKM_EXEC vtkm::Vec<U, 3> ConvertVecFromSpace(const vtkm::Vec<U, 2> vec) const
{
return vec[0] * this->Basis0 + vec[1] * this->Basis1;
}
};
} //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_WEDGE(pc, rc, call) \
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(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]))
//-----------------------------------------------------------------------------
// This returns the Jacobian of a hexahedron's (or other 3D cell's) coordinates
// with respect to parametric coordinates. Explicitly, this is (d is partial
// derivative):
//
// | |
// | dx/du dx/dv dx/dw |
// | |
// | dy/du dy/dv dy/dw |
// | |
// | dz/du dz/dv dz/dw |
// | |
//
#define VTKM_ACCUM_JACOBIAN_3D(pointIndex, weight0, weight1, weight2) \
jacobian(0, 0) += static_cast<JacobianType>(wCoords[pointIndex][0] * (weight0)); \
jacobian(1, 0) += static_cast<JacobianType>(wCoords[pointIndex][1] * (weight0)); \
jacobian(2, 0) += static_cast<JacobianType>(wCoords[pointIndex][2] * (weight0)); \
jacobian(0, 1) += static_cast<JacobianType>(wCoords[pointIndex][0] * (weight1)); \
jacobian(1, 1) += static_cast<JacobianType>(wCoords[pointIndex][1] * (weight1)); \
jacobian(2, 1) += static_cast<JacobianType>(wCoords[pointIndex][2] * (weight1)); \
jacobian(0, 2) += static_cast<JacobianType>(wCoords[pointIndex][0] * (weight2)); \
jacobian(1, 2) += static_cast<JacobianType>(wCoords[pointIndex][1] * (weight2)); \
jacobian(2, 2) += static_cast<JacobianType>(wCoords[pointIndex][2] * (weight2))
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);
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);
}
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::CellShapeTagWedge)
{
vtkm::Vec<JacobianType, 3> pc(pcoords);
vtkm::Vec<JacobianType, 3> rc = vtkm::Vec<JacobianType, 3>(JacobianType(1)) - pc;
jacobian = vtkm::Matrix<JacobianType, 3, 3>(JacobianType(0));
VTKM_DERIVATIVE_WEIGHTS_WEDGE(pc, rc, VTKM_ACCUM_JACOBIAN_3D);
}
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::CellShapeTagPyramid)
{
vtkm::Vec<JacobianType, 3> pc(pcoords);
vtkm::Vec<JacobianType, 3> rc = vtkm::Vec<JacobianType, 3>(JacobianType(1)) - pc;
jacobian = vtkm::Matrix<JacobianType, 3, 3>(JacobianType(0));
VTKM_DERIVATIVE_WEIGHTS_PYRAMID(pc, rc, VTKM_ACCUM_JACOBIAN_3D);
}
// Derivatives in quadrilaterals are computed in much the same way as
// hexahedra. Review the documentation for hexahedra derivatives for details
// on the math. The major difference is that the equations are performed in
// a 2D space built with make_SpaceForQuadrilateral.
#define VTKM_ACCUM_JACOBIAN_2D(pointIndex, weight0, weight1) \
wcoords2d = space.ConvertCoordToSpace(wCoords[pointIndex]); \
jacobian(0, 0) += wcoords2d[0] * (weight0); \
jacobian(1, 0) += wcoords2d[1] * (weight0); \
jacobian(0, 1) += wcoords2d[0] * (weight1); \
jacobian(1, 1) += wcoords2d[1] * (weight1)
template <typename WorldCoordType,
typename ParametricCoordType,
typename SpaceType,
typename JacobianType>
VTKM_EXEC void JacobianFor2DCell(const WorldCoordType& wCoords,
const vtkm::Vec<ParametricCoordType, 3>& pcoords,
const vtkm::exec::internal::Space2D<SpaceType>& space,
vtkm::Matrix<JacobianType, 2, 2>& jacobian,
vtkm::CellShapeTagQuad)
{
vtkm::Vec<JacobianType, 2> pc(static_cast<JacobianType>(pcoords[0]),
static_cast<JacobianType>(pcoords[1]));
vtkm::Vec<JacobianType, 2> rc = vtkm::Vec<JacobianType, 2>(JacobianType(1)) - pc;
vtkm::Vec<SpaceType, 2> wcoords2d;
jacobian = vtkm::Matrix<JacobianType, 2, 2>(JacobianType(0));
VTKM_DERIVATIVE_WEIGHTS_QUAD(pc, rc, VTKM_ACCUM_JACOBIAN_2D);
}
#if 0
// This code doesn't work, so I'm bailing on it. Instead, I'm just grabbing a
// triangle and finding the derivative of that. If you can do better, please
// implement it.
template<typename WorldCoordType,
typename ParametricCoordType,
typename JacobianType>
VTKM_EXEC
void JacobianFor2DCell(const WorldCoordType &wCoords,
const vtkm::Vec<ParametricCoordType,3> &pcoords,
const vtkm::exec::internal::Space2D<JacobianType> &space,
vtkm::Matrix<JacobianType,2,2> &jacobian,
vtkm::CellShapeTagPolygon)
{
const vtkm::IdComponent numPoints = wCoords.GetNumberOfComponents();
vtkm::Vec<JacobianType,2> pc(pcoords[0], pcoords[1]);
JacobianType deltaAngle = 2*vtkm::Pi<JacobianType>()/numPoints;
jacobian = vtkm::Matrix<JacobianType,2,2>(0);
for (vtkm::IdComponent pointIndex = 0; pointIndex < numPoints; pointIndex++)
{
JacobianType angle = pointIndex*deltaAngle;
vtkm::Vec<JacobianType,2> nodePCoords(0.5f*(vtkm::Cos(angle)+1),
0.5f*(vtkm::Sin(angle)+1));
// This is the vector pointing from the user provided parametric coordinate
// to the node at pointIndex in parametric space.
vtkm::Vec<JacobianType,2> pvec = nodePCoords - pc;
// The weight (the derivative of the interpolation factor) happens to be
// pvec scaled by the cube root of pvec's magnitude.
JacobianType magSqr = vtkm::MagnitudeSquared(pvec);
JacobianType invMag = vtkm::RSqrt(magSqr);
JacobianType scale = invMag*invMag*invMag;
vtkm::Vec<JacobianType,2> weight = scale*pvec;
vtkm::Vec<JacobianType,2> wcoords2d =
space.ConvertCoordToSpace(wCoords[pointIndex]);
jacobian(0,0) += wcoords2d[0] * weight[0];
jacobian(1,0) += wcoords2d[1] * weight[0];
jacobian(0,1) += wcoords2d[0] * weight[1];
jacobian(1,1) += wcoords2d[1] * weight[1];
}
}
#endif
#undef VTKM_ACCUM_JACOBIAN_3D
#undef VTKM_ACCUM_JACOBIAN_2D
}
} // namespace vtkm::exec
#endif //vtk_m_exec_Jacobian_h

1090
vtkm/exec/ParametricCoordinates.h
View File

File diff suppressed because it is too large Load Diff

28
vtkm/thirdparty/vtkc/CMakeLists.txt vendored Normal file
View File

@ -0,0 +1,28 @@
##============================================================================
## Copyright (c) Kitware, Inc.
## All rights reserved.
## See LICENSE.txt for details.
##
## This software is distributed WITHOUT ANY WARRANTY; without even
## the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
## PURPOSE. See the above copyright notice for more information.
##============================================================================
add_library(vtkm_vtkc INTERFACE)
vtkm_get_kit_name(kit_name kit_dir)
# vtkc needs C++11
target_compile_features(vtkm_vtkc INTERFACE cxx_std_11)
target_include_directories(vtkm_vtkc INTERFACE
$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/vtkmvtkc>
$<INSTALL_INTERFACE:${VTKm_INSTALL_INCLUDE_DIR}/vtkm/thirdparty/vtkc/vtkmvtkc>)
install(TARGETS vtkm_vtkc
EXPORT ${VTKm_EXPORT_NAME})
## Install headers
if(NOT VTKm_INSTALL_ONLY_LIBRARIES)
install(DIRECTORY vtkmvtkc
DESTINATION ${VTKm_INSTALL_INCLUDE_DIR}/${kit_dir}/)
endif()

6
vtkm/worklet/testing/UnitTestParticleAdvection.cxx
View File

@ -317,7 +317,7 @@ void ValidateEvaluator(const EvalType& eval,
Status status = statusPortal.Get(index);
vtkm::Vec3f result = resultsPortal.Get(index);
VTKM_TEST_ASSERT(status == Status::SUCCESS, "Error in evaluator for " + msg);
VTKM_TEST_ASSERT(result == vec, "Error in evaluator result for " + msg);
VTKM_TEST_ASSERT(test_equal(result, vec), "Error in evaluator result for " + msg);
}
pointsHandle.ReleaseResources();
evalStatus.ReleaseResources();
@ -370,10 +370,10 @@ void ValidateIntegrator(const IntegratorType& integrator,
vtkm::Vec3f result = resultsPortal.Get(index);
VTKM_TEST_ASSERT(status != Status::FAIL, "Error in evaluator for " + msg);
if (status == Status::OUTSIDE_SPATIAL_BOUNDS)
VTKM_TEST_ASSERT(result == pointsPortal.Get(index),
VTKM_TEST_ASSERT(test_equal(result, pointsPortal.Get(index)),
"Error in evaluator result for [OUTSIDE SPATIAL]" + msg);
else
VTKM_TEST_ASSERT(result == expStepResults[(size_t)index],
VTKM_TEST_ASSERT(test_equal(result, expStepResults[(size_t)index]),
"Error in evaluator result for " + msg);
}
pointsHandle.ReleaseResources();

2
vtkm/worklet/testing/UnitTestTemporalAdvection.cxx
View File

@ -84,7 +84,7 @@ void ValidateEvaluator(const EvalType& eval,
vtkm::Vec<ScalarType, 3> result = resultsPortal.Get(index);
vtkm::Vec<ScalarType, 3> expected = validityPortal.Get(index);
VTKM_TEST_ASSERT(status == Status::SUCCESS, "Error in evaluator for " + msg);
VTKM_TEST_ASSERT(result == expected, "Error in evaluator result for " + msg);
VTKM_TEST_ASSERT(test_equal(result, expected), "Error in evaluator result for " + msg);
}
evalStatus.ReleaseResources();
evalResults.ReleaseResources();