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 'use-vtkc'

Browse Source 
e08d862f9 Use Light-Weight Cell Library
e56b34fc9 Merge branch 'upstream-vtkc' into use-vtkc
e909edd61 vtkc 2019-09-24 (cdc72582)
2a6aadc69 Add update file for VTK-c

Acked-by: Kenneth Moreland <kmorel@sandia.gov>
Merge-request: !1852
This commit is contained in:
Brad King 2019-09-25 09:19:26 -04:00
commit 9a906e4cbf
34 changed files with 5086 additions and 2589 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
.gitattributes vendored
View File

@ -17,3 +17,4 @@ data/** filter=lfs diff=lfs merge=lfs -text
vtkm/thirdparty/diy/vtkmdiy/** -format.clang-format -whitespace
vtkm/thirdparty/optionparser/vtkmoptionparser/** -format.clang-format -whitespace
vtkm/thirdparty/lodepng/vtkmlodepng/** -format.clang-format -whitespace
vtkm/thirdparty/vtkc/vtkmvtkc/** -format.clang-format -whitespace

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()

25
vtkm/thirdparty/vtkc/update.sh vendored Executable file
View File

@ -0,0 +1,25 @@
#!/usr/bin/env bash
set -e
set -x
shopt -s dotglob
readonly name="vtkc"
readonly ownership="VTK-c Upstream <kwrobot@kitware.com>"
readonly subtree="vtkm/thirdparty/$name/vtkm$name"
readonly repo="https://gitlab.kitware.com/sujin.philip/vtk-c.git"
readonly tag="master"
readonly paths="
vtkc
LICENSE.md
README.md
"
extract_source () {
git_archive
pushd "${extractdir}/${name}-reduced"
rm -rf vtkc/testing
popd
}
. "${BASH_SOURCE%/*}/../update-common.sh"

43
vtkm/thirdparty/vtkc/vtkmvtkc/LICENSE.md vendored Normal file
View File

@ -0,0 +1,43 @@
# VTK-c License Version 1.0 #
Copyright (c) 2019
Kitware Inc.,
National Technology & Engineering Solutions of Sandia, LLC (NTESS),
UT-Battelle, LLC.,
Los Alamos National Security, LLC.,
All rights reserved.
Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government
retains certain rights in this software.
Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
Laboratory (LANL), the U.S. Government retains certain rights in
this software.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the
distribution.
* Neither the name of Kitware nor the names of any contributors may
be used to endorse or promote products derived from this software
without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHORS OR
CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

65
vtkm/thirdparty/vtkc/vtkmvtkc/README.md vendored Normal file
View File

@ -0,0 +1,65 @@
# VTK-c #
VTK-c is a lightweight collection of cell types and cell functionality
that is designed be used scientific visualization libraries of
VTK-m, and VTK.
You can find out more about the design of VTK-c in [DESIGN.md].
## Contributing ##
There are many ways to contribute to [VTK-c]:
+ Submit new or add to discussions of a feature requests or bugs on the
[VTK-c Issue Tracker].
+ Submit a Pull Request
+ See [CONTRIBUTING.md] for detailed instructions on how to create a
Pull Request.
+ See the [VTK-c Coding Conventions] that must be followed for
contributed code.
## Compiler Requirements ##
+ C++11 Compiler. VTK-c has been confirmed to work with the following
+ GCC 4.8+
+ Clang 3.3+
+ XCode 5.0+
+ MSVC 2015+
## Example##
Below is a simple example of using VTK-c to get derivatives and
parametric coordinates for different cell types:
```cpp
#incude <vtkc/vtkc.h>
std::array<float, 3> pcoords;
auto status = vtkc::parametricCenter(vtkc::Hexahedron{}, pcoords);
std::vector<std::array<float, 3>> points = { ... };
std::vector<std::array<double, 4>> field = { ... };
std::array<double, 4> derivs[3];
status = vtkc::derivative(vtkc::Hexahedron{},
vtkc::makeFieldAccessorNestedSOAConst(points, 3),
vtkc::makeFieldAccessorNestedSOAConst(field, 4),
pcoords,
derivs[0],
derivs[1],
derivs[2]);
```
## License ##
VTK-c is distributed under the OSI-approved BSD 3-clause License.
See [LICENSE.md] for details.
[VTK-c]: https://gitlab.kitware.com/sujin.philip/vtk-c/
[VTK-c Issue Tracker]: https://gitlab.kitware.com/sujin.philip/vtk-c/issues
[CONTRIBUTING.md]: CONTRIBUTING.md
[DESIGN.md]: docs/Design.md
[LICENSE.md]: LICENSE.md
[VTK-c Coding Conventions]: docs/CodingConventions.md

69
vtkm/thirdparty/vtkc/vtkmvtkc/vtkc/ErrorCode.h vendored Normal file
View File

@ -0,0 +1,69 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.md 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_c_ErrorCode_h
#define vtk_c_ErrorCode_h
#include <vtkc/internal/Config.h>
#include <cstdint>
namespace vtkc
{
enum class ErrorCode : std::int32_t
{
SUCCESS = 0,
INVALID_SHAPE_ID,
INVALID_NUMBER_OF_POINTS,
WRONG_SHAPE_ID_FOR_TAG_TYPE,
INVALID_POINT_ID,
SOLUTION_DID_NOT_CONVERGE,
MATRIX_LUP_FACTORIZATION_FAILED,
DEGENERATE_CELL_DETECTED
};
VTKC_EXEC
inline const char* errorString(ErrorCode code) noexcept
{
switch (code)
{
case ErrorCode::SUCCESS:
return "Success";
case ErrorCode::INVALID_SHAPE_ID:
return "Invalid shape id";
case ErrorCode::INVALID_NUMBER_OF_POINTS:
return "Invalid number of points";
case ErrorCode::WRONG_SHAPE_ID_FOR_TAG_TYPE:
return "Wrong shape id for tag type";
case ErrorCode::INVALID_POINT_ID:
return "Invalid point id";
case ErrorCode::SOLUTION_DID_NOT_CONVERGE:
return "Solution did not converge";
case ErrorCode::MATRIX_LUP_FACTORIZATION_FAILED:
return "LUP factorization failed";
case ErrorCode::DEGENERATE_CELL_DETECTED:
return "Degenerate cell detected";
}
return "Invalid error";
}
} // vtkc
#define VTKC_RETURN_ON_ERROR(call) \
{ \
auto status = call; \
if (status != vtkc::ErrorCode::SUCCESS) \
{ \
return status; \
} \
}
#endif // vtk_c_ErrorCode_h

367
vtkm/thirdparty/vtkc/vtkmvtkc/vtkc/FieldAccessor.h vendored Normal file
View File

@ -0,0 +1,367 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.md 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_c_FieldAccessor_h
#define vtk_c_FieldAccessor_h
#include <vtkc/internal/Config.h>
#include <utility>
namespace vtkc
{
namespace internal
{
//----------------------------------------------------------------------------
/// Does T support the indexing operator?
/// Examples:
/// std::array
/// std::vector
/// T*
/// custom classes with operator[](integer_type)
///
struct NotVecTag {};
template <typename T> inline decltype(std::declval<T>()[0]) VecCheck(T);
inline NotVecTag VecCheck(...);
template <typename T>
using IsVecType =
std::integral_constant<bool, !std::is_same<decltype(VecCheck(std::declval<T>())),
internal::NotVecTag>::value>;
//----------------------------------------------------------------------------
/// Does T support operator()(integer_type)?
struct NotFuncTag {};
template <typename T> inline decltype(std::declval<T>()(0)) FuncCheck(T);
inline NotFuncTag FuncCheck(...);
template <typename T>
using IsFuncType =
std::integral_constant<bool, !std::is_same<decltype(FuncCheck(std::declval<T>())),
internal::NotFuncTag>::value>;
//----------------------------------------------------------------------------
template <typename T> typename T::size_type IndexTypeImpl(T);
IdComponent IndexTypeImpl(...);
template <typename VecT>
using IndexType = decltype(IndexTypeImpl(std::declval<VecT>()));
//----------------------------------------------------------------------------
struct VecTypeIndexer {};
struct VecTypeFunctor {};
struct VecTypeScalar {};
template <typename T>
using VecTypeTag =
typename std::conditional<IsVecType<T>::value,
VecTypeIndexer,
typename std::conditional<IsFuncType<T>::value,
VecTypeFunctor,
VecTypeScalar>::type>::type;
//----------------------------------------------------------------------------
template <typename T>
VTKC_EXEC auto componentImpl(VecTypeIndexer, T&& vec, int idx) -> decltype(vec[0])
{
return vec[static_cast<IndexType<T>>(idx)];
}
template <typename T>
VTKC_EXEC auto componentImpl(VecTypeFunctor, T&& vec, int idx) -> decltype(vec(0))
{
return vec(idx);
}
template <typename T>
VTKC_EXEC T& componentImpl(VecTypeScalar, T&& vec, int)
{
return vec;
}
} // namespace internal
//----------------------------------------------------------------------------
template <typename T>
VTKC_EXEC auto component(T&& vec, int idx)
-> decltype(internal::componentImpl<T>(internal::VecTypeTag<T>{}, std::forward<T>(vec), 0))
{
return internal::componentImpl<T>(internal::VecTypeTag<T>{}, std::forward<T>(vec), idx);
}
//----------------------------------------------------------------------------
template <typename T>
using ComponentType = typename std::decay<decltype(component(std::declval<T>(), 0))>::type;
///============================================================================
/// Since there are different ways fields maybe represented in the clients of
/// this libarary, VTK-c relies on helper classes that implement the
/// \c FieldAccessor "concept" to access the elements of a field.
///
/// These classes should wrap the field and provide the follolwing interface:
///
/// template <typename FieldType>
/// class FieldAccessor
/// {
/// public:
/// /// An alias for the component type of the field
/// using ValueType = ...;
///
/// /// Return the number of components
/// int getNumberOfComponents() const;
///
/// /// Set the value to `tuple` and `comp` to `value`.
/// void setValue(int tuple, int comp, const ValueType& value) const;
///
/// /// Get the value at `tuple` and `comp`.
/// ValueType getValue(int tuple, int comp) const;
///
/// /// Set the tuple at index `tuple`. It is recomended to make this a
/// /// template function and use `vtkc::component` to access
/// /// the components of `value`.
/// template <typename VecType>
/// void setTuple(int tuple, const VecType& value) const;
///
/// /// Get the tuple at index `tuple`. It is recomended to make this a
/// /// template function and use `vtkc::component` to access
/// /// the components of `value`.
/// template <typename VecType>
/// void getTuple(int tuple, VecType& value) const;
/// }
///
/// The set functions are optional and such a class would act as a const variant
/// of \c FieldAccessor.
///
///----------------------------------------------------------------------------
template <typename FieldType>
class FieldAccessorNestedSOA
{
public:
using ValueType = typename std::decay<decltype(component(std::declval<FieldType>()[0], 0))>::type;
VTKC_EXEC FieldAccessorNestedSOA(FieldType& field, int numberOfComponents = 1)
: Field(&field), NumberOfComponents(numberOfComponents)
{
}
VTKC_EXEC int getNumberOfComponents() const
{
return this->NumberOfComponents;
}
VTKC_EXEC void setValue(int tuple, int comp, const ValueType& value) const
{
component((*this->Field)[static_cast<IdxType>(tuple)], comp) = value;
}
VTKC_EXEC ValueType getValue(int tuple, int comp) const
{
return component((*this->Field)[static_cast<IdxType>(tuple)], comp);
}
template <typename VecType>
VTKC_EXEC void setTuple(int tuple, const VecType& value) const
{
for (int i = 0; i < this->NumberOfComponents; ++i)
{
component((*this->Field)[static_cast<IdxType>(tuple)], i) =
static_cast<ValueType>(component(value, i));
}
}
template <typename VecType>
VTKC_EXEC void getTuple(int tuple, VecType& value) const
{
for (int i = 0; i < this->NumberOfComponents; ++i)
{
component(value, i) = static_cast<ComponentType<VecType>>(
component((*this->Field)[static_cast<IdxType>(tuple)], i));
}
}
private:
using IdxType = internal::IndexType<FieldType>;
FieldType* Field;
int NumberOfComponents;
};
template <typename FieldType>
class FieldAccessorNestedSOAConst
{
public:
using ValueType = typename std::decay<decltype(component(std::declval<FieldType>()[0], 0))>::type;
VTKC_EXEC FieldAccessorNestedSOAConst(const FieldType& field, int numberOfComponents = 1)
: Field(&field), NumberOfComponents(numberOfComponents)
{
}
VTKC_EXEC int getNumberOfComponents() const
{
return this->NumberOfComponents;
}
VTKC_EXEC ValueType getValue(int tuple, int comp) const
{
return component((*this->Field)[static_cast<IdxType>(tuple)], comp);
}
template <typename VecType>
VTKC_EXEC void getTuple(int tuple, VecType& value) const
{
for (int i = 0; i < this->NumberOfComponents; ++i)
{
component(value, i) = static_cast<ComponentType<VecType>>(
component((*this->Field)[static_cast<IdxType>(tuple)], i));
}
}
private:
using IdxType = internal::IndexType<FieldType>;
const FieldType* Field;
int NumberOfComponents;
};
template <typename FieldType>
VTKC_EXEC
FieldAccessorNestedSOA<FieldType> makeFieldAccessorNestedSOA(FieldType& field,
int numberOfComponents = 1)
{
return FieldAccessorNestedSOA<FieldType>(field, numberOfComponents);
}
template <typename FieldType>
VTKC_EXEC
FieldAccessorNestedSOAConst<FieldType> makeFieldAccessorNestedSOAConst(const FieldType& field,
int numberOfComponents = 1)
{
return FieldAccessorNestedSOAConst<FieldType>(field, numberOfComponents);
}
///----------------------------------------------------------------------------
template <typename FieldType>
class FieldAccessorFlatSOA
{
public:
using ValueType = typename std::decay<decltype(std::declval<FieldType>()[0])>::type;
VTKC_EXEC FieldAccessorFlatSOA(FieldType& field, int numberOfComponents = 1)
: Field(&field), NumberOfComponents(numberOfComponents)
{
}
VTKC_EXEC int getNumberOfComponents() const
{
return this->NumberOfComponents;
}
VTKC_EXEC void setValue(int tuple, int comp, const ValueType& value) const
{
auto FlatIdx = static_cast<IdxType>(tuple * this->NumberOfComponents + comp);
(*this->Field)[FlatIdx] = value;
}
VTKC_EXEC ValueType getValue(int tuple, int comp) const
{
auto FlatIdx = static_cast<IdxType>(tuple * this->NumberOfComponents + comp);
return (*this->Field)[FlatIdx];
}
template <typename VecType>
VTKC_EXEC void setTuple(int tuple, const VecType& value) const
{
auto start = static_cast<IdxType>(tuple * this->NumberOfComponents);
for (int i = 0; i < this->NumberOfComponents; ++i)
{
(*this->Field)[start++] = static_cast<ValueType>(component(value, i));
}
}
template <typename VecType>
VTKC_EXEC void getTuple(int tuple, VecType& value) const
{
auto start = static_cast<IdxType>(tuple * this->NumberOfComponents);
for (int i = 0; i < this->NumberOfComponents; ++i)
{
component(value, i) = static_cast<ComponentType<VecType>>((*this->Field)[start++]);
}
}
private:
using IdxType = internal::IndexType<FieldType>;
FieldType* Field;
int NumberOfComponents;
};
template <typename FieldType>
class FieldAccessorFlatSOAConst
{
public:
using ValueType = typename std::decay<decltype(std::declval<FieldType>()[0])>::type;
VTKC_EXEC FieldAccessorFlatSOAConst(const FieldType& field, int numberOfComponents = 1)
: Field(&field), NumberOfComponents(numberOfComponents)
{
}
VTKC_EXEC int getNumberOfComponents() const
{
return this->NumberOfComponents;
}
VTKC_EXEC ValueType getValue(int tuple, int comp) const
{
auto FlatIdx = static_cast<IdxType>(tuple * this->NumberOfComponents + comp);
return (*this->Field)[FlatIdx];
}
template <typename VecType>
VTKC_EXEC void getTuple(int tuple, VecType& value) const
{
auto start = static_cast<IdxType>(tuple * this->NumberOfComponents);
for (int i = 0; i < this->NumberOfComponents; ++i)
{
component(value, i) = static_cast<ComponentType<VecType>>((*this->Field)[start++]);
}
}
private:
using IdxType = internal::IndexType<FieldType>;
const FieldType* Field;
int NumberOfComponents;
};
template <typename FieldType>
VTKC_EXEC
FieldAccessorFlatSOA<FieldType> makeFieldAccessorFlatSOA(FieldType& field,
int numberOfComponents)
{
return FieldAccessorFlatSOA<FieldType>(field, numberOfComponents);
}
template <typename FieldType>
VTKC_EXEC
FieldAccessorFlatSOAConst<FieldType> makeFieldAccessorFlatSOAConst(const FieldType& field,
int numberOfComponents)
{
return FieldAccessorFlatSOAConst<FieldType>(field, numberOfComponents);
}
} // namespace vtkc
#endif // vtk_c_FieldAccessor_h

254
vtkm/thirdparty/vtkc/vtkmvtkc/vtkc/Hexahedron.h vendored Normal file
View File

@ -0,0 +1,254 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.md 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_c_Hexahedron_h
#define vtk_c_Hexahedron_h
#include <vtkc/ErrorCode.h>
#include <vtkc/Shapes.h>
#include <vtkc/internal/Common.h>
namespace vtkc
{
class Hexahedron : public Cell
{
public:
constexpr VTKC_EXEC Hexahedron() : Cell(ShapeId::HEXAHEDRON, 8) {}
constexpr VTKC_EXEC explicit Hexahedron(const Cell& cell) : Cell(cell) {}
};
VTKC_EXEC inline vtkc::ErrorCode validate(Hexahedron tag) noexcept
{
if (tag.shape() != ShapeId::HEXAHEDRON && tag.shape() != ShapeId::VOXEL)
{
return ErrorCode::WRONG_SHAPE_ID_FOR_TAG_TYPE;
}
if (tag.numberOfPoints() != 8)
{
return ErrorCode::INVALID_NUMBER_OF_POINTS;
}
return ErrorCode::SUCCESS;
}
template<typename CoordType>
VTKC_EXEC inline vtkc::ErrorCode parametricCenter(Hexahedron, CoordType&& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
component(pcoords, 0) = 0.5f;
component(pcoords, 1) = 0.5f;
component(pcoords, 2) = 0.5f;
return ErrorCode::SUCCESS;
}
template<typename CoordType>
VTKC_EXEC inline vtkc::ErrorCode parametricPoint(
Hexahedron, IdComponent pointId, CoordType&& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
switch (pointId)
{
case 0:
component(pcoords, 0) = 0.0f;
component(pcoords, 1) = 0.0f;
component(pcoords, 2) = 0.0f;
break;
case 1:
component(pcoords, 0) = 1.0f;
component(pcoords, 1) = 0.0f;
component(pcoords, 2) = 0.0f;
break;
case 2:
component(pcoords, 0) = 1.0f;
component(pcoords, 1) = 1.0f;
component(pcoords, 2) = 0.0f;
break;
case 3:
component(pcoords, 0) = 0.0f;
component(pcoords, 1) = 1.0f;
component(pcoords, 2) = 0.0f;
break;
case 4:
component(pcoords, 0) = 0.0f;
component(pcoords, 1) = 0.0f;
component(pcoords, 2) = 1.0f;
break;
case 5:
component(pcoords, 0) = 1.0f;
component(pcoords, 1) = 0.0f;
component(pcoords, 2) = 1.0f;
break;
case 6:
component(pcoords, 0) = 1.0f;
component(pcoords, 1) = 1.0f;
component(pcoords, 2) = 1.0f;
break;
case 7:
component(pcoords, 0) = 0.0f;
component(pcoords, 1) = 1.0f;
component(pcoords, 2) = 1.0f;
break;
default:
return ErrorCode::INVALID_POINT_ID;
}
return ErrorCode::SUCCESS;
}
template<typename CoordType>
VTKC_EXEC inline ComponentType<CoordType> parametricDistance(Hexahedron, const CoordType& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
return internal::findParametricDistance(pcoords, 3);
}
template<typename CoordType>
VTKC_EXEC inline bool cellInside(Hexahedron, const CoordType& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
using T = ComponentType<CoordType>;
return component(pcoords, 0) >= T{0} && component(pcoords, 0) <= T{1} &&
component(pcoords, 1) >= T{0} && component(pcoords, 1) <= T{1} &&
component(pcoords, 2) >= T{0} && component(pcoords, 2) <= T{1};
}
template <typename Values, typename CoordType, typename Result>
VTKC_EXEC inline vtkc::ErrorCode interpolate(
Hexahedron,
const Values& values,
const CoordType& pcoords,
Result&& result) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
using T = internal::ClosestFloatType<typename Values::ValueType>;
for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
{
auto vbf = internal::lerp(static_cast<T>(values.getValue(0, c)),
static_cast<T>(values.getValue(1, c)),
static_cast<T>(component(pcoords, 0)));
auto vbb = internal::lerp(static_cast<T>(values.getValue(3, c)),
static_cast<T>(values.getValue(2, c)),
static_cast<T>(component(pcoords, 0)));
auto vtf = internal::lerp(static_cast<T>(values.getValue(4, c)),
static_cast<T>(values.getValue(5, c)),
static_cast<T>(component(pcoords, 0)));
auto vtb = internal::lerp(static_cast<T>(values.getValue(7, c)),
static_cast<T>(values.getValue(6, c)),
static_cast<T>(component(pcoords, 0)));
auto vb = internal::lerp(vbf, vbb, static_cast<T>(component(pcoords, 1)));
auto vt = internal::lerp(vtf, vtb, static_cast<T>(component(pcoords, 1)));
auto v = internal::lerp(vb, vt, static_cast<T>(component(pcoords, 2)));
component(result, c) = static_cast<ComponentType<Result>>(v);
}
return ErrorCode::SUCCESS;
}
namespace internal
{
template <typename Values, typename CoordType, typename Result>
VTKC_EXEC inline void parametricDerivative(Hexahedron,
const Values& values,
IdComponent comp,
const CoordType& pcoords,
Result&& result) noexcept
{
using T = internal::ClosestFloatType<typename Values::ValueType>;
T p0 = static_cast<T>(component(pcoords, 0));
T p1 = static_cast<T>(component(pcoords, 1));
T p2 = static_cast<T>(component(pcoords, 2));
T rm = T{1} - p0;
T sm = T{1} - p1;
T tm = T{1} - p2;
T dr = (static_cast<T>(values.getValue(0, comp)) * -sm * tm) +
(static_cast<T>(values.getValue(1, comp)) * sm * tm) +
(static_cast<T>(values.getValue(2, comp)) * p1 * tm) +
(static_cast<T>(values.getValue(3, comp)) * -p1 * tm) +
(static_cast<T>(values.getValue(4, comp)) * -sm * p2) +
(static_cast<T>(values.getValue(5, comp)) * sm * p2) +
(static_cast<T>(values.getValue(6, comp)) * p1 * p2) +
(static_cast<T>(values.getValue(7, comp)) * -p1 * p2);
T ds = (static_cast<T>(values.getValue(0, comp)) * -rm * tm) +
(static_cast<T>(values.getValue(1, comp)) * -p0 * tm) +
(static_cast<T>(values.getValue(2, comp)) * p0 * tm) +
(static_cast<T>(values.getValue(3, comp)) * rm * tm) +
(static_cast<T>(values.getValue(4, comp)) * -rm * p2) +
(static_cast<T>(values.getValue(5, comp)) * -p0 * p2) +
(static_cast<T>(values.getValue(6, comp)) * p0 * p2) +
(static_cast<T>(values.getValue(7, comp)) * rm * p2);
T dt = (static_cast<T>(values.getValue(0, comp)) * -rm * sm) +
(static_cast<T>(values.getValue(1, comp)) * -p0 * sm) +
(static_cast<T>(values.getValue(2, comp)) * -p0 * p1) +
(static_cast<T>(values.getValue(3, comp)) * -rm * p1) +
(static_cast<T>(values.getValue(4, comp)) * rm * sm) +
(static_cast<T>(values.getValue(5, comp)) * p0 * sm) +
(static_cast<T>(values.getValue(6, comp)) * p0 * p1) +
(static_cast<T>(values.getValue(7, comp)) * rm * p1);
component(result, 0) = static_cast<ComponentType<Result>>(dr);
component(result, 1) = static_cast<ComponentType<Result>>(ds);
component(result, 2) = static_cast<ComponentType<Result>>(dt);
}
} // internal
template <typename Points, typename Values, typename CoordType, typename Result>
VTKC_EXEC inline vtkc::ErrorCode derivative(
Hexahedron,
const Points& points,
const Values& values,
const CoordType& pcoords,
Result&& dx,
Result&& dy,
Result&& dz) noexcept
{
return internal::derivative3D(Hexahedron{},
points,
values,
pcoords,
std::forward<Result>(dx),
std::forward<Result>(dy),
std::forward<Result>(dz));
}
template <typename Points, typename PCoordType, typename WCoordType>
VTKC_EXEC inline vtkc::ErrorCode parametricToWorld(
Hexahedron,
const Points& points,
const PCoordType& pcoords,
WCoordType&& wcoords) noexcept
{
return interpolate(Hexahedron{}, points, pcoords, std::forward<WCoordType>(wcoords));
}
template <typename Points, typename WCoordType, typename PCoordType>
VTKC_EXEC inline vtkc::ErrorCode worldToParametric(
Hexahedron,
const Points& points,
const WCoordType& wcoords,
PCoordType&& pcoords) noexcept
{
return internal::worldToParametric3D(
Hexahedron{}, points, wcoords, std::forward<PCoordType>(pcoords));
}
} // vtkc
#endif // vtk_c_Hexahedron_h

176
vtkm/thirdparty/vtkc/vtkmvtkc/vtkc/Line.h vendored Normal file
View File

@ -0,0 +1,176 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.md 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_c_Line_h
#define vtk_c_Line_h
#include <vtkc/ErrorCode.h>
#include <vtkc/Shapes.h>
#include <vtkc/internal/Common.h>
namespace vtkc
{
class Line : public Cell
{
public:
constexpr VTKC_EXEC Line() : Cell(ShapeId::LINE, 2) {}
constexpr VTKC_EXEC explicit Line(const Cell& cell) : Cell(cell) {}
};
VTKC_EXEC inline vtkc::ErrorCode validate(Line tag) noexcept
{
if (tag.shape() != ShapeId::LINE)
{
return ErrorCode::WRONG_SHAPE_ID_FOR_TAG_TYPE;
}
if (tag.numberOfPoints() != 2)
{
return ErrorCode::INVALID_NUMBER_OF_POINTS;
}
return ErrorCode::SUCCESS;
}
template<typename CoordType>
VTKC_EXEC inline vtkc::ErrorCode parametricCenter(Line, CoordType&& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
component(pcoords, 0) = 0.5f;
return ErrorCode::SUCCESS;
}
template<typename CoordType>
VTKC_EXEC inline vtkc::ErrorCode parametricPoint(
Line, IdComponent pointId, CoordType&& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
switch (pointId)
{
case 0:
component(pcoords, 0) = 0.0f;
return ErrorCode::SUCCESS;
case 1:
component(pcoords, 0) = 1.0f;
return ErrorCode::SUCCESS;
default:
return ErrorCode::INVALID_POINT_ID;
}
}
template<typename CoordType>
VTKC_EXEC inline ComponentType<CoordType> parametricDistance(Line, const CoordType& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
return internal::findParametricDistance(pcoords, 1);
}
template<typename CoordType>
VTKC_EXEC inline bool cellInside(Line, const CoordType& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
using T = ComponentType<CoordType>;
return component(pcoords, 0) >= T{0} && component(pcoords, 0) <= T{1};
}
template <typename Values, typename CoordType, typename Result>
VTKC_EXEC inline vtkc::ErrorCode interpolate(
Line,
const Values& values,
const CoordType& pcoords,
Result&& result) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
using ProcessingType = internal::ClosestFloatType<typename Values::ValueType>;
using ResultCompType = ComponentType<Result>;
for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
{
auto ival = internal::lerp(static_cast<ProcessingType>(values.getValue(0, c)),
static_cast<ProcessingType>(values.getValue(1, c)),
static_cast<ProcessingType>(component(pcoords, 0)));
component(result, c) = static_cast<ResultCompType>(ival);
}
return ErrorCode::SUCCESS;
}
template <typename Points, typename Values, typename CoordType, typename Result>
VTKC_EXEC inline vtkc::ErrorCode derivative(
Line,
const Points& points,
const Values& values,
const CoordType&,
Result&& dx,
Result&& dy,
Result&& dz) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
using ProcessingType = internal::ClosestFloatType<typename Values::ValueType>;
using ResultCompType = ComponentType<Result>;
ProcessingType dPt[3] = {
static_cast<ProcessingType>(points.getValue(1, 0) - points.getValue(0, 0)),
static_cast<ProcessingType>(points.getValue(1, 1) - points.getValue(0, 1)),
static_cast<ProcessingType>(points.getValue(1, 2) - points.getValue(0, 2)) };
for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
{
auto dv = static_cast<ProcessingType>(values.getValue(1, c) - values.getValue(0, c));
component(dx, c) = (dPt[0] != 0.0f) ? static_cast<ResultCompType>(dv/dPt[0]) : ResultCompType{0};
component(dy, c) = (dPt[1] != 0.0f) ? static_cast<ResultCompType>(dv/dPt[1]) : ResultCompType{0};
component(dz, c) = (dPt[2] != 0.0f) ? static_cast<ResultCompType>(dv/dPt[2]) : ResultCompType{0};
}
return ErrorCode::SUCCESS;
}
template <typename Points, typename PCoordType, typename WCoordType>
VTKC_EXEC inline vtkc::ErrorCode parametricToWorld(
Line,
const Points& points,
const PCoordType& pcoords,
WCoordType&& wcoords) noexcept
{
return interpolate(Line{}, points, pcoords, std::forward<WCoordType>(wcoords));
}
template <typename Points, typename WCoordType, typename PCoordType>
VTKC_EXEC inline vtkc::ErrorCode worldToParametric(
Line,
const Points& points,
const WCoordType& wcoords,
PCoordType&& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(PCoordType);
using T = ComponentType<PCoordType>;
internal::Vector<T, 3> p0(static_cast<T>(points.getValue(0, 0)),
static_cast<T>(points.getValue(0, 1)),
static_cast<T>(points.getValue(0, 2)));
internal::Vector<T, 3> p1(static_cast<T>(points.getValue(1, 0)),
static_cast<T>(points.getValue(1, 1)),
static_cast<T>(points.getValue(1, 2)));
internal::Vector<T, 3> wc(static_cast<T>(component(wcoords, 0)),
static_cast<T>(component(wcoords, 1)),
static_cast<T>(component(wcoords, 2)));
auto v1 = p1 - p0;
auto v2 = wc - p0;
component(pcoords, 0) = internal::dot(v1, v2) / internal::dot(v1, v1);
return ErrorCode::SUCCESS;
}
} //namespace vtkc
#endif //vtk_c_Line_h

178
vtkm/thirdparty/vtkc/vtkmvtkc/vtkc/Pixel.h vendored Normal file
View File

@ -0,0 +1,178 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.md 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_c_Pixel_h
#define vtk_c_Pixel_h
#include <vtkc/ErrorCode.h>
#include <vtkc/Quad.h>
#include <vtkc/internal/Common.h>
namespace vtkc
{
class Pixel : public Quad
{
public:
constexpr VTKC_EXEC Pixel() : Quad(Cell(ShapeId::PIXEL, 4)) {}
constexpr VTKC_EXEC explicit Pixel(const Cell& cell) : Quad(cell) {}
};
namespace internal
{
template <typename Points, typename T>
VTKC_EXEC inline int getPixelSpacing(const Points& points, T spacing[3])
{
int zeros = 0;
for (int i = 0; i < 3; ++i)
{
spacing[i] = static_cast<T>(points.getValue(2, i) - points.getValue(0, i));
if (spacing[i] == T{0})
{
zeros |= 1 << i;
}
}
return zeros;
}
} // internal
template <typename Points, typename Values, typename CoordType, typename Result>
VTKC_EXEC inline vtkc::ErrorCode derivative(
Pixel,
const Points& points,
const Values& values,
const CoordType& pcoords,
Result&& dx,
Result&& dy,
Result&& dz) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
using ProcessingType = internal::ClosestFloatType<typename Values::ValueType>;
using ResultCompType = ComponentType<Result>;
ProcessingType spacing[3];
int zeros = internal::getPixelSpacing(points, spacing);
for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
{
ProcessingType dvdp[2];
internal::parametricDerivative(Quad{}, values, c, pcoords, dvdp);
switch (zeros)
{
case 1: // yz plane
component(dx, c) = ResultCompType{0};
component(dy, c) = static_cast<ResultCompType>(dvdp[0] / spacing[1]);
component(dz, c) = static_cast<ResultCompType>(dvdp[1] / spacing[2]);
break;
case 2: // xz plane
component(dx, c) = static_cast<ResultCompType>(dvdp[0] / spacing[0]);
component(dy, c) = ResultCompType{0};
component(dz, c) = static_cast<ResultCompType>(dvdp[1] / spacing[2]);
break;
case 4: // xy plane
component(dx, c) = static_cast<ResultCompType>(dvdp[0] / spacing[0]);
component(dy, c) = static_cast<ResultCompType>(dvdp[1] / spacing[1]);
component(dz, c) = ResultCompType{0};
break;
default:
return ErrorCode::DEGENERATE_CELL_DETECTED;
}
}
return ErrorCode::SUCCESS;
}
template <typename Points, typename PCoordType, typename WCoordType>
VTKC_EXEC vtkc::ErrorCode parametricToWorld(
Pixel,
const Points& points,
const PCoordType& pcoords,
WCoordType&& wcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(PCoordType);
using T = typename Points::ValueType;
T spacing[3];
auto zeros = internal::getPixelSpacing(points, spacing);
switch (zeros)
{
case 1: // yz plane
component(wcoords, 0) = points.getValue(0, 0);
component(wcoords, 1) = points.getValue(0, 1) +
(spacing[1] * static_cast<T>(component(pcoords, 0)));
component(wcoords, 2) = points.getValue(0, 2) +
(spacing[2] * static_cast<T>(component(pcoords, 1)));
return ErrorCode::SUCCESS;
case 2: // xz plane
component(wcoords, 0) = points.getValue(0, 0) +
(spacing[0] * static_cast<T>(component(pcoords, 0)));
component(wcoords, 1) = points.getValue(0, 1);
component(wcoords, 2) = points.getValue(0, 2) +
(spacing[2] * static_cast<T>(component(pcoords, 1)));
return ErrorCode::SUCCESS;
case 4: // xy plane
component(wcoords, 0) = points.getValue(0, 0) +
(spacing[0] * static_cast<T>(component(pcoords, 0)));
component(wcoords, 1) = points.getValue(0, 1) +
(spacing[1] * static_cast<T>(component(pcoords, 1)));
component(wcoords, 2) = points.getValue(0, 2);
return ErrorCode::SUCCESS;
default:
return ErrorCode::DEGENERATE_CELL_DETECTED;
}
}
template <typename Points, typename WCoordType, typename PCoordType>
VTKC_EXEC vtkc::ErrorCode worldToParametric(
Pixel,
const Points& points,
const WCoordType& wcoords,
PCoordType&& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(PCoordType);
using T = ComponentType<PCoordType>;
T spacing[3];
int zeros = internal::getPixelSpacing(points, spacing);
switch (zeros)
{
case 1: // yz plane
component(pcoords, 0) = static_cast<T>(component(wcoords, 1) - points.getValue(0, 1)) /
spacing[1];
component(pcoords, 1) = static_cast<T>(component(wcoords, 2) - points.getValue(0, 2)) /
spacing[2];
return ErrorCode::SUCCESS;
case 2: // xz plane
component(pcoords, 0) = static_cast<T>(component(wcoords, 0) - points.getValue(0, 0)) /
spacing[0];
component(pcoords, 1) = static_cast<T>(component(wcoords, 2) - points.getValue(0, 2)) /
spacing[2];
return ErrorCode::SUCCESS;
case 4: // xy plane
component(pcoords, 0) = static_cast<T>(component(wcoords, 0) - points.getValue(0, 0)) /
spacing[0];
component(pcoords, 1) = static_cast<T>(component(wcoords, 1) - points.getValue(0, 1)) /
spacing[1];
return ErrorCode::SUCCESS;
default:
return ErrorCode::DEGENERATE_CELL_DETECTED;
}
}
} //namespace vtkc
#endif //vtk_c_Pixel_h

541
vtkm/thirdparty/vtkc/vtkmvtkc/vtkc/Polygon.h vendored Normal file
View File

@ -0,0 +1,541 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.md 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_c_Polygon_h
#define vtk_c_Polygon_h
#include <vtkc/ErrorCode.h>
#include <vtkc/Quad.h>
#include <vtkc/Shapes.h>
#include <vtkc/Triangle.h>
#include <vtkc/internal/Common.h>
namespace vtkc
{
/// \c Polygon with 3 and 4 points behave exactly as \c Triangle and \c Quad
/// respectively. For 5 or more points, the points are arranged 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 the 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.
class Polygon : public Cell
{
public:
constexpr VTKC_EXEC Polygon() : Cell(ShapeId::POLYGON, 3) {}
constexpr VTKC_EXEC explicit Polygon(vtkc::IdComponent numPoints)
: Cell(ShapeId::POLYGON, numPoints)
{
}
constexpr VTKC_EXEC explicit Polygon(const Cell& cell) : Cell(cell) {}
};
VTKC_EXEC inline vtkc::ErrorCode validate(Polygon tag) noexcept
{
if (tag.shape() != ShapeId::POLYGON)
{
return ErrorCode::WRONG_SHAPE_ID_FOR_TAG_TYPE;
}
if (tag.numberOfPoints() < 3)
{
return ErrorCode::INVALID_NUMBER_OF_POINTS;
}
return ErrorCode::SUCCESS;
}
template<typename CoordType>
VTKC_EXEC inline vtkc::ErrorCode parametricCenter(Polygon tag, CoordType&& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
switch (tag.numberOfPoints())
{
case 3:
return parametricCenter(Triangle{}, pcoords);
case 4:
return parametricCenter(Quad{}, pcoords);
default:
component(pcoords, 0) = 0.5f;
component(pcoords, 1) = 0.5f;
return ErrorCode::SUCCESS;
}
}
template<typename CoordType>
VTKC_EXEC inline vtkc::ErrorCode parametricPoint(
Polygon tag, IdComponent pointId, CoordType&& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
if (pointId < 0 || pointId >= tag.numberOfPoints())
{
return ErrorCode::INVALID_POINT_ID;
}
switch (tag.numberOfPoints())
{
case 3:
return parametricPoint(Triangle{}, pointId, pcoords);
case 4:
return parametricPoint(Quad{}, pointId, pcoords);
default:
{
using T = ComponentType<CoordType>;
constexpr double two_pi = 2.0 * 3.14159265359;
auto angle = (static_cast<T>(pointId) * static_cast<T>(two_pi)) / static_cast<T>(tag.numberOfPoints());
component(pcoords, 0) = 0.5f * (VTKC_MATH_CALL(cos, (angle)) + 1.0f);
component(pcoords, 1) = 0.5f * (VTKC_MATH_CALL(sin, (angle)) + 1.0f);
return ErrorCode::SUCCESS;
}
}
}
template<typename CoordType>
VTKC_EXEC inline ComponentType<CoordType> parametricDistance(Polygon tag, const CoordType& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
switch (tag.numberOfPoints())
{
case 3:
return parametricDistance(Triangle{}, pcoords);
default:
return internal::findParametricDistance(pcoords, 2);
}
}
template<typename CoordType>
VTKC_EXEC inline bool cellInside(Polygon tag, const CoordType& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
using T = ComponentType<CoordType>;
switch (tag.numberOfPoints())
{
case 3:
return cellInside(Triangle{}, pcoords);
case 4:
return cellInside(Quad{}, pcoords);
default:
return (((component(pcoords, 0) - T(0.5f)) * (component(pcoords, 0) - T(0.5f))) +
((component(pcoords, 1) - T(0.5f)) * (component(pcoords, 1) - T(0.5f)))) <= T(0.25f);
}
}
namespace internal
{
template <typename CoordType>
VTKC_EXEC inline vtkc::ErrorCode polygonToSubTrianglePCoords(
Polygon tag,
const CoordType& polygonPC,
IdComponent& p0,
IdComponent& p1,
ComponentType<CoordType> trianglePC[2]) noexcept
{
using T = ComponentType<CoordType>;
constexpr T epsilon = std::is_same<T, float>::value ? T(1e-5f) : T(1e-9f);
// Find the sub-triangle containing pcoords
auto x = component(polygonPC, 1) - T(0.5f);
auto y = component(polygonPC, 0) - T(0.5f);
if (VTKC_MATH_CALL(abs, (x)) < (T(4) * epsilon) && VTKC_MATH_CALL(abs, (y)) < (T(4) * epsilon))
{
// we are at the center
p0 = 0;
p1 = 1;
trianglePC[0] = trianglePC[1] = T(0);
return ErrorCode::SUCCESS;
}
constexpr double two_pi = 2.0 * 3.14159265359;
T angle = VTKC_MATH_CALL(atan2, (x), (y));
if (angle < T(0))
{
angle += static_cast<T>(two_pi);
}
T deltaAngle = static_cast<T>(two_pi) / static_cast<T>(tag.numberOfPoints());
p0 = static_cast<IdComponent>(VTKC_MATH_CALL(floor, (angle / deltaAngle)));
p1 = (p0 + 1) % tag.numberOfPoints();
// Build triangle with polygon pcoords as its wcoords
T triPts[9] = { T(0), T(0), T(0), T(0), T(0), T(0), T(0), T(0), T(0) };
VTKC_RETURN_ON_ERROR(parametricCenter(tag, triPts))
VTKC_RETURN_ON_ERROR(parametricPoint(tag, p0, triPts + 3))
VTKC_RETURN_ON_ERROR(parametricPoint(tag, p1, triPts + 6))
// Find the parametric coord on the triangle
T triWC[3] = { component(polygonPC, 0), component(polygonPC, 1), T(0) };
VTKC_RETURN_ON_ERROR(worldToParametric(Triangle{}, makeFieldAccessorFlatSOAConst(triPts, 3), triWC, trianglePC))
return ErrorCode::SUCCESS;
}
template <typename Values>
VTKC_EXEC inline typename Values::ValueType polygonInterpolateComponentAtCenter(
Polygon tag, const Values& values, IdComponent comp) noexcept
{
using T = internal::ClosestFloatType<typename Values::ValueType>;
auto weight = T{1} / static_cast<T>(tag.numberOfPoints());
auto result = static_cast<T>(values.getValue(0, comp));
for (IdComponent i = 1; i < tag.numberOfPoints(); ++i)
{
result += static_cast<T>(values.getValue(i, comp));
}
result *= weight;
return static_cast<typename Values::ValueType>(result);
}
} // namespace internal
template <typename Values, typename CoordType, typename Result>
VTKC_EXEC vtkc::ErrorCode interpolate(
Polygon tag,
const Values& values,
const CoordType& pcoords,
Result&& result) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
switch (tag.numberOfPoints())
{
case 3:
return interpolate(Triangle{}, values, pcoords, std::forward<Result>(result));
case 4:
return interpolate(Quad{}, values, pcoords, std::forward<Result>(result));
default:
break;
}
using ResultCompType = ComponentType<Result>;
using ProcessingType = internal::ClosestFloatType<typename Values::ValueType>;
IdComponent p0, p1;
ComponentType<CoordType> triPc[2];
VTKC_RETURN_ON_ERROR(internal::polygonToSubTrianglePCoords(tag, pcoords, p0, p1, triPc))
// compute polygon interpolation from triangle weights
for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
{
ProcessingType triVals[3];
triVals[0] = static_cast<ProcessingType>(internal::polygonInterpolateComponentAtCenter(tag, values, c));
triVals[1] = static_cast<ProcessingType>(values.getValue(p0, c));
triVals[2] = static_cast<ProcessingType>(values.getValue(p1, c));
ResultCompType val = 0;
VTKC_RETURN_ON_ERROR(interpolate(Triangle{}, makeFieldAccessorNestedSOA(triVals), triPc, &val))
component(result, c) = val;
}
return ErrorCode::SUCCESS;
}
namespace internal
{
// To find the gradient in a polygon (of 5 or more points), we will extract a small triangle near
// the desired parameteric coordinates (pcoords). We return the field values (outField) and world
// coordinates (outWCoords) for this triangle, which is all that is needed to find the gradient
// in a triangle.
//
// The trangle will be "pointing" away from the center of the polygon, and pcoords will be placed
// at the apex of the triangle. This is because if pcoords is at or near the edge of the polygon,
// we do not want to push any of the points over the edge, and it is not trivial to determine
// exactly where the edge of the polygon is.
template <typename CoordType>
VTKC_EXEC inline void polygonGetTriangleAroundPCoords(
const CoordType& pcoords, ComponentType<CoordType> pc1[2], ComponentType<CoordType> pc2[2]) noexcept
{
using T = ComponentType<CoordType>;
// Find the unit vector pointing from the center of the polygon to pcoords
Vector<T, 2> radialVector(component(pcoords, 0) - 0.5f, component(pcoords, 1) - 0.5f);
auto magSqr = dot(radialVector, radialVector);
if (magSqr > 8.0f * 1e-4f)
{
radialVector /= VTKC_MATH_CALL(sqrt, (magSqr));
}
else
{
// pcoords is in the center of the polygon. Just point in an arbitrary direction
radialVector[0] = T(1);
radialVector[1] = T(0);
}
// We want the two points away from pcoords to be back toward the center but moved at 45 degrees
// off the radius. Simple geometry shows us that the (not quite unit) vectors of those two
// directions are (-radialVector[1] - radialVector[0], radialVector[0] - radialVector[1]) and
// (radialVector[1] - radialVector[0], -radialVector[0] - radialVector[1]).
//
// *\ (-radialVector[1], radialVector[0]) //
// | \ //
// | \ (-radialVector[1] - radialVector[0], radialVector[0] - radialVector[1]) //
// | \ //
// +-------* radialVector //
// | / //
// | / (radialVector[1] - radialVector[0], -radialVector[0] - radialVector[1]) //
// | / //
// */ (radialVector[1], -radialVector[0]) //
// This scaling value is somewhat arbitrary. It is small enough to be "close" to the selected
// point and small enough to be guaranteed to be inside the polygon, but large enough to
// get an accurate gradient.
static constexpr T scale = 0.05f;
pc1[0] = pcoords[0] + scale * (-radialVector[1] - radialVector[0]);
pc1[1] = pcoords[1] + scale * (radialVector[0] - radialVector[1]);
pc2[0] = pcoords[0] + scale * (radialVector[1] - radialVector[0]);
pc2[1] = pcoords[1] + scale * (-radialVector[0] - radialVector[1]);
}
} // namespace internal
template <typename Points, typename Values, typename CoordType, typename Result>
VTKC_EXEC inline vtkc::ErrorCode derivative(
Polygon tag,
const Points& points,
const Values& values,
const CoordType& pcoords,
Result&& dx,
Result&& dy,
Result&& dz) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
switch (tag.numberOfPoints())
{
case 3:
return derivative(Triangle{},
points,
values,
pcoords,
std::forward<Result>(dx),
std::forward<Result>(dy),
std::forward<Result>(dz));
case 4:
return derivative(Quad{},
points,
values,
pcoords,
std::forward<Result>(dx),
std::forward<Result>(dy),
std::forward<Result>(dz));
default:
break;
}
using ResultCompType = ComponentType<Result>;
using ProcessingType = internal::ClosestFloatType<typename Values::ValueType>;
// Get the parametric coordinates of a small triangle, with pcoords as one of the vertices
ComponentType<CoordType> ptPc1[2], ptPc2[2];
internal::polygonGetTriangleAroundPCoords(pcoords, ptPc1, ptPc2);
// Compute world coordinates of the points of the triangle
internal::Vector<ProcessingType, 3> triPts[3];
VTKC_RETURN_ON_ERROR(interpolate(tag, points, pcoords, triPts[0]))
VTKC_RETURN_ON_ERROR(interpolate(tag, points, ptPc1, triPts[1]))
VTKC_RETURN_ON_ERROR(interpolate(tag, points, ptPc2, triPts[2]))
// Compute the derivative on the triangle
//----------------------------------------
// 2-D coordinate system on the triangle's plane
internal::Space2D<ProcessingType> triSpace(triPts[0], triPts[1], triPts[2]);
internal::Vector<ProcessingType, 2> pts2d[3];
for (int i = 0; i < 3; ++i)
{
pts2d[i] = triSpace.to2DPoint(triPts[i]);
}
// pre-compute once
internal::Matrix<ProcessingType, 2, 2> jacobian;
internal::jacobian2D(Triangle{}, makeFieldAccessorNestedSOA(pts2d, 2), nullptr, jacobian);
internal::Matrix<ProcessingType, 2, 2> invJacobian;
VTKC_RETURN_ON_ERROR(internal::matrixInverse(jacobian, invJacobian))
// Compute sub-triangle information of the three vertices of the derivation triangle to
// reduce the amount of redundant computations in the loop.
IdComponent subP1P2[3][2];
ComponentType<CoordType> pcs[3][2];
internal::polygonToSubTrianglePCoords(tag, pcoords, subP1P2[0][0], subP1P2[0][1], pcs[0]);
internal::polygonToSubTrianglePCoords(tag, ptPc1, subP1P2[1][0], subP1P2[1][1], pcs[1]);
internal::polygonToSubTrianglePCoords(tag, ptPc2, subP1P2[2][0], subP1P2[2][1], pcs[2]);
for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
{
// Interpolate component values at the vertices of the derivation triangle.
auto vCenter = static_cast<ProcessingType>(internal::polygonInterpolateComponentAtCenter(tag, values, c));
ProcessingType triVals[3];
for (int i = 0; i < 3; ++i)
{
ProcessingType field[3] = {vCenter,
static_cast<ProcessingType>(values.getValue(subP1P2[i][0], c)),
static_cast<ProcessingType>(values.getValue(subP1P2[i][1], c))};
VTKC_RETURN_ON_ERROR(interpolate(Triangle{}, makeFieldAccessorNestedSOA(field), pcs[i], triVals + i))
}
// Compute derivative in the triangle
internal::Vector<ProcessingType, 2> dvdp;
parametricDerivative(Triangle{}, makeFieldAccessorNestedSOA(triVals), 0, nullptr, dvdp);
auto d2D = matrixMultiply(dvdp, invJacobian);
auto d3D = triSpace.to3DVec(d2D);
component(dx, c) = static_cast<ResultCompType>(d3D[0]);
component(dy, c) = static_cast<ResultCompType>(d3D[1]);
component(dz, c) = static_cast<ResultCompType>(d3D[2]);
}
return ErrorCode::SUCCESS;
}
template <typename Points, typename PCoordType, typename WCoordType>
VTKC_EXEC inline vtkc::ErrorCode parametricToWorld(
Polygon tag,
const Points& points,
const PCoordType& pcoords,
WCoordType&& wcoords) noexcept
{
return interpolate(tag, points, pcoords, std::forward<WCoordType>(wcoords));
}
template <typename Points, typename WCoordType, typename PCoordType>
VTKC_EXEC inline vtkc::ErrorCode worldToParametric(
Polygon tag,
const Points& points,
const WCoordType& wcoords,
PCoordType&& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(PCoordType);
switch (tag.numberOfPoints())
{
case 3:
return worldToParametric(Triangle{}, points, wcoords, pcoords);
case 4:
return worldToParametric(Quad{}, points, wcoords, pcoords);
default:
break;
}
using T = ComponentType<WCoordType>;
auto numPoints = tag.numberOfPoints();
// Find the position of the center point.
internal::Vector<T, 3> wcoordCenter{T(0), T(0), T(0)};
for (IdComponent pointIndex = 0; pointIndex < numPoints; ++pointIndex)
{
wcoordCenter[0] += points.getValue(pointIndex, 0);
wcoordCenter[1] += points.getValue(pointIndex, 1);
wcoordCenter[2] += points.getValue(pointIndex, 2);
}
wcoordCenter /= static_cast<T>(numPoints);
// Find the normal vector to the polygon. If the polygon is planar, convex,
// and in general position, any three points will give a normal in the same
// direction. Although not perfectly robust, we can reduce the effect of
// non-planar, non-convex, or degenerate polygons by picking three points
// topologically far from each other. Note that we do not care about the
// length of the normal in this case.
internal::Vector<T, 3> polygonNormal;
{
internal::Vector<T, 3> v1p1, v1p2;
points.getTuple(0, v1p1);
points.getTuple(numPoints / 3, v1p2);
internal::Vector<T, 3> v2p1, v2p2;
points.getTuple(1, v2p1);
points.getTuple(2 * numPoints / 3, v2p2);
polygonNormal = internal::cross(v1p2 - v1p1, v2p2 - v2p1);
}
// Find which triangle wcoords is located in. We do this by defining the
// equations for the planes through the radial edges and perpendicular to the
// polygon. The point is in the triangle if it is on the correct side of both
// planes.
internal::Vector<T, 3> wc{ component(wcoords, 0), component(wcoords, 1), component(wcoords, 2) };
IdComponent firstPointIndex;
IdComponent secondPointIndex = 0;
internal::Vector<T, 3> firstPoint, secondPoint;
bool foundTriangle = false;
for (firstPointIndex = 0; firstPointIndex < numPoints - 1; ++firstPointIndex)
{
points.getTuple(firstPointIndex, firstPoint);
auto vecInPlane = firstPoint - wcoordCenter;
auto planeNormal = internal::cross(polygonNormal, vecInPlane);
auto planeOffset = internal::dot(planeNormal, wcoordCenter);
if (internal::dot(planeNormal, wc) < planeOffset)
{
// wcoords on wrong side of plane, thus outside of triangle
continue;
}
secondPointIndex = firstPointIndex + 1;
points.getTuple(secondPointIndex, secondPoint);
vecInPlane = secondPoint - wcoordCenter;
planeNormal = internal::cross(polygonNormal, vecInPlane);
planeOffset = internal::dot(planeNormal, wcoordCenter);
if (internal::dot(planeNormal, wc) > planeOffset)
{
// wcoords on wrong side of plane, thus outside of triangle
continue;
}
foundTriangle = true;
break;
}
if (!foundTriangle)
{
// wcoord was outside of all triangles we checked. It must be inside the
// one triangle we did not check (the one between the first and last
// polygon points).
firstPointIndex = numPoints - 1;
points.getTuple(firstPointIndex, firstPoint);
secondPointIndex = 0;
points.getTuple(secondPointIndex, secondPoint);
}
// Build a structure containing the points of the triangle wcoords is in and
// use the triangle version of this function to find the parametric
// coordinates.
internal::Vector<T, 3> triangleWCoords[3] = { wcoordCenter, firstPoint, secondPoint };
internal::Vector<T, 3> trianglePCoords;
VTKC_RETURN_ON_ERROR(worldToParametric(
Triangle{}, makeFieldAccessorNestedSOA(triangleWCoords, 3), wc, trianglePCoords))
// trianglePCoords is in the triangle's parameter space rather than the
// polygon's parameter space. We can find the polygon's parameter space by
// repurposing parametricToWorld by using the
// polygon parametric coordinates as a proxy for world coordinates.
VTKC_RETURN_ON_ERROR(parametricCenter(tag, triangleWCoords[0]))
VTKC_RETURN_ON_ERROR(parametricPoint(tag, firstPointIndex, triangleWCoords[1]))
VTKC_RETURN_ON_ERROR(parametricPoint(tag, secondPointIndex, triangleWCoords[2]))
triangleWCoords[0][2] = triangleWCoords[1][2] = triangleWCoords[2][2] = T(0);
VTKC_RETURN_ON_ERROR(
parametricToWorld(Triangle{}, makeFieldAccessorNestedSOA(triangleWCoords, 3), trianglePCoords, wc))
component(pcoords, 0) = static_cast<ComponentType<PCoordType>>(wc[0]);
component(pcoords, 1) = static_cast<ComponentType<PCoordType>>(wc[1]);
return ErrorCode::SUCCESS;
}
} //namespace vtkc
#endif //vtk_c_Polygon_h

293
vtkm/thirdparty/vtkc/vtkmvtkc/vtkc/Pyramid.h vendored Normal file
View File

@ -0,0 +1,293 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.md 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_c_Pyramid_h
#define vtk_c_Pyramid_h
#include <vtkc/ErrorCode.h>
#include <vtkc/Shapes.h>
#include <vtkc/internal/Common.h>
namespace vtkc
{
class Pyramid : public Cell
{
public:
constexpr VTKC_EXEC Pyramid() : Cell(ShapeId::PYRAMID, 5) {}
constexpr VTKC_EXEC explicit Pyramid(const Cell& cell) : Cell(cell) {}
};
VTKC_EXEC inline vtkc::ErrorCode validate(Pyramid tag) noexcept
{
if (tag.shape() != ShapeId::PYRAMID)
{
return ErrorCode::WRONG_SHAPE_ID_FOR_TAG_TYPE;
}
if (tag.numberOfPoints() != 5)
{
return ErrorCode::INVALID_NUMBER_OF_POINTS;
}
return ErrorCode::SUCCESS;
}
template<typename CoordType>
VTKC_EXEC inline vtkc::ErrorCode parametricCenter(Pyramid, CoordType&& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
component(pcoords, 0) = 0.5f;
component(pcoords, 1) = 0.5f;
component(pcoords, 2) = 0.2f;
return ErrorCode::SUCCESS;
}
template<typename CoordType>
VTKC_EXEC inline vtkc::ErrorCode parametricPoint(
Pyramid, IdComponent pointId, CoordType&& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
switch (pointId)
{
case 0:
component(pcoords, 0) = 0.0f;
component(pcoords, 1) = 0.0f;
component(pcoords, 2) = 0.0f;
break;
case 1:
component(pcoords, 0) = 1.0f;
component(pcoords, 1) = 0.0f;
component(pcoords, 2) = 0.0f;
break;
case 2:
component(pcoords, 0) = 1.0f;
component(pcoords, 1) = 1.0f;
component(pcoords, 2) = 0.0f;
break;
case 3:
component(pcoords, 0) = 0.0f;
component(pcoords, 1) = 1.0f;
component(pcoords, 2) = 0.0f;
break;
case 4:
component(pcoords, 0) = 0.5f;
component(pcoords, 1) = 0.5f;
component(pcoords, 2) = 1.0f;
break;
default:
return ErrorCode::INVALID_POINT_ID;
}
return ErrorCode::SUCCESS;
}
template<typename CoordType>
VTKC_EXEC inline ComponentType<CoordType> parametricDistance(Pyramid, const CoordType& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
return internal::findParametricDistance(pcoords, 3);
}
template<typename CoordType>
VTKC_EXEC inline bool cellInside(Pyramid, const CoordType& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
using T = ComponentType<CoordType>;
return component(pcoords, 0) >= T{0} && component(pcoords, 0) <= T{1} &&
component(pcoords, 1) >= T{0} && component(pcoords, 1) <= T{1} &&
component(pcoords, 2) >= T{0} && component(pcoords, 2) <= T{1};
}
template <typename Values, typename CoordType, typename Result>
VTKC_EXEC inline vtkc::ErrorCode interpolate(
Pyramid,
const Values& values,
const CoordType& pcoords,
Result&& result) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
using T = internal::ClosestFloatType<typename Values::ValueType>;
for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
{
auto baseV0 = internal::lerp(static_cast<T>(values.getValue(0, c)),
static_cast<T>(values.getValue(1, c)),
static_cast<T>(component(pcoords, 0)));
auto baseV1 = internal::lerp(static_cast<T>(values.getValue(3, c)),
static_cast<T>(values.getValue(2, c)),
static_cast<T>(component(pcoords, 0)));
auto baseV = internal::lerp(baseV0, baseV1, static_cast<T>(component(pcoords, 1)));
auto v = internal::lerp(baseV,
static_cast<T>(values.getValue(4, c)),
static_cast<T>(component(pcoords, 2)));
component(result, c) = static_cast<ComponentType<Result>>(v);
}
return ErrorCode::SUCCESS;
}
namespace internal
{
template <typename Values, typename CoordType, typename Result>
VTKC_EXEC inline void parametricDerivative(Pyramid,
const Values& values,
IdComponent comp,
const CoordType& pcoords,
Result&& result) noexcept
{
using T = internal::ClosestFloatType<typename Values::ValueType>;
T p0 = static_cast<T>(component(pcoords, 0));
T p1 = static_cast<T>(component(pcoords, 1));
T p2 = static_cast<T>(component(pcoords, 2));
T rm = T{1} - p0;
T sm = T{1} - p1;
T tm = T{1} - p2;
T dr = (static_cast<T>(values.getValue(0, comp)) * -sm * tm) +
(static_cast<T>(values.getValue(1, comp)) * sm * tm) +
(static_cast<T>(values.getValue(2, comp)) * p1 * tm) +
(static_cast<T>(values.getValue(3, comp)) * -p1 * tm);
T ds = (static_cast<T>(values.getValue(0, comp)) * -rm * tm) +
(static_cast<T>(values.getValue(1, comp)) * -p0 * tm) +
(static_cast<T>(values.getValue(2, comp)) * p0 * tm) +
(static_cast<T>(values.getValue(3, comp)) * rm * tm);
T dt = (static_cast<T>(values.getValue(0, comp)) * -rm * sm) +
(static_cast<T>(values.getValue(1, comp)) * -p0 * sm) +
(static_cast<T>(values.getValue(2, comp)) * -p0 * p1) +
(static_cast<T>(values.getValue(3, comp)) * -rm * p1) +
(static_cast<T>(values.getValue(4, comp)));
component(result, 0) = static_cast<ComponentType<Result>>(dr);
component(result, 1) = static_cast<ComponentType<Result>>(ds);
component(result, 2) = static_cast<ComponentType<Result>>(dt);
}
} // internal
template <typename Points, typename Values, typename CoordType, typename Result>
VTKC_EXEC inline vtkc::ErrorCode derivative(
Pyramid,
const Points& points,
const Values& values,
const CoordType& pcoords,
Result&& dx,
Result&& dy,
Result&& dz) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
using ProcessingType = internal::ClosestFloatType<typename Values::ValueType>;
using ResultCompType = ComponentType<Result>;
if (component(pcoords, 2) > ComponentType<CoordType>(.999f))
{
// If we are at the apex of the pyramid we need to do something special.
// As we approach the apex, the derivatives of the parametric shape
// functions in x and y go to 0 while the inverse of the Jacobian
// also goes to 0. This results in 0/0 but using l'Hopital's rule
// we could actually compute the value of the limit, if we had a
// functional expression to compute the gradient. We're on a computer
// so we don't but we can cheat and do a linear extrapolation of the
// derivatives which really ends up as the same thing.
internal::Matrix<ProcessingType, 3, 3> j, ij1, ij2;
const ComponentType<CoordType> pc1[3] = {0.5f, 0.5f, (2.0f * 0.998f) - component(pcoords, 2)};
internal::jacobian3D(Pyramid{}, points, pc1, j);
VTKC_RETURN_ON_ERROR(internal::matrixInverse(j, ij1))
const ComponentType<CoordType> pc2[3] = {0.5f, 0.5f, 0.998f};
internal::jacobian3D(Pyramid{}, points, pc2, j);
VTKC_RETURN_ON_ERROR(internal::matrixInverse(j, ij2))
for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
{
internal::Vector<ProcessingType, 3> dvdp;
internal::parametricDerivative(Pyramid{}, values, c, pc1, dvdp);
auto d1 = matrixMultiply(dvdp, ij1);
internal::parametricDerivative(Pyramid{}, values, c, pc2, dvdp);
auto d2 = matrixMultiply(dvdp, ij2);
component(dx, c) = static_cast<ResultCompType>((d2[0] * 2.0f) - d1[0]);
component(dy, c) = static_cast<ResultCompType>((d2[1] * 2.0f) - d1[1]);
component(dz, c) = static_cast<ResultCompType>((d2[2] * 2.0f) - d1[2]);
}
}
else
{
return internal::derivative3D(Pyramid{},
points,
values,
pcoords,
std::forward<Result>(dx),
std::forward<Result>(dy),
std::forward<Result>(dz));
}
return ErrorCode::SUCCESS;
}
template <typename Points, typename PCoordType, typename WCoordType>
VTKC_EXEC inline vtkc::ErrorCode parametricToWorld(
Pyramid,
const Points& points,
const PCoordType& pcoords,
WCoordType&& wcoords) noexcept
{
return interpolate(Pyramid{}, points, pcoords, std::forward<WCoordType>(wcoords));
}
template <typename Points, typename WCoordType, typename PCoordType>
VTKC_EXEC inline vtkc::ErrorCode worldToParametric(
Pyramid,
const Points& points,
const WCoordType& wcoords,
PCoordType&& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(PCoordType);
using TIn = typename Points::ValueType;
using TOut = ComponentType<PCoordType>;
internal::Vector<TIn, 3> wcVec{component(wcoords, 0), component(wcoords, 1), component(wcoords, 2)};
// Newton's method fails if the wcoord is too close to the apex. Just return the pcoords at the
// apex for those cases.
internal::Vector<TOut, 3> pcBaseCenter(0.5f, 0.5f, 0.0f);
internal::Vector<TIn, 3> apex, wcBaseCenter;
points.getTuple(4, apex);
VTKC_RETURN_ON_ERROR(parametricToWorld(Pyramid{}, points, pcBaseCenter, wcBaseCenter))
auto apexToBase = wcBaseCenter - apex;
auto apexToWc = wcVec - apex;
auto dist2ApexToBase = internal::dot(apexToBase, apexToBase);
auto dist2ApexToWC = internal::dot(apexToWc, apexToWc);
if (dist2ApexToWC <= (1e-6f * dist2ApexToBase))
{
return parametricPoint(Pyramid{}, 4, pcoords);
}
return internal::worldToParametric3D(
Pyramid{}, points, wcoords, std::forward<PCoordType>(pcoords));
}
} // vtkc
#endif // vtk_c_Pyramid_h

195
vtkm/thirdparty/vtkc/vtkmvtkc/vtkc/Quad.h vendored Normal file
View File

@ -0,0 +1,195 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.md 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_c_Quad_h
#define vtk_c_Quad_h
#include <vtkc/ErrorCode.h>
#include <vtkc/Shapes.h>
#include <vtkc/internal/Common.h>
namespace vtkc
{
class Quad : public Cell
{
public:
constexpr VTKC_EXEC Quad() : Cell(ShapeId::QUAD, 4) {}
constexpr VTKC_EXEC explicit Quad(const Cell& cell) : Cell(cell) {}
};
VTKC_EXEC inline vtkc::ErrorCode validate(Quad tag) noexcept
{
if (tag.shape() != ShapeId::QUAD && tag.shape() != ShapeId::PIXEL)
{
return ErrorCode::WRONG_SHAPE_ID_FOR_TAG_TYPE;
}
if (tag.numberOfPoints() != 4)
{
return ErrorCode::INVALID_NUMBER_OF_POINTS;
}
return ErrorCode::SUCCESS;
}
template<typename CoordType>
VTKC_EXEC inline vtkc::ErrorCode parametricCenter(Quad, CoordType&& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
component(pcoords, 0) = 0.5f;
component(pcoords, 1) = 0.5f;
return ErrorCode::SUCCESS;
}
template<typename CoordType>
VTKC_EXEC inline vtkc::ErrorCode parametricPoint(
Quad, IdComponent pointId, CoordType&& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
switch (pointId)
{
case 0:
component(pcoords, 0) = 0.0f;
component(pcoords, 1) = 0.0f;
break;
case 1:
component(pcoords, 0) = 1.0f;
component(pcoords, 1) = 0.0f;
break;
case 2:
component(pcoords, 0) = 1.0f;
component(pcoords, 1) = 1.0f;
break;
case 3:
component(pcoords, 0) = 0.0f;
component(pcoords, 1) = 1.0f;
break;
default:
return ErrorCode::INVALID_POINT_ID;
}
return ErrorCode::SUCCESS;
}
template<typename CoordType>
VTKC_EXEC inline ComponentType<CoordType> parametricDistance(Quad, const CoordType& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
return internal::findParametricDistance(pcoords, 2);
}
template<typename CoordType>
VTKC_EXEC inline bool cellInside(Quad, const CoordType& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
using T = ComponentType<CoordType>;
return component(pcoords, 0) >= T{0} && component(pcoords, 0) <= T{1} &&
component(pcoords, 1) >= T{0} && component(pcoords, 1) <= T{1};
}
template <typename Values, typename CoordType, typename Result>
VTKC_EXEC inline vtkc::ErrorCode interpolate(
Quad,
const Values& values,
const CoordType& pcoords,
Result&& result) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
using T = internal::ClosestFloatType<typename Values::ValueType>;
for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
{
auto v0 = internal::lerp(static_cast<T>(values.getValue(0, c)),
static_cast<T>(values.getValue(1, c)),
static_cast<T>(component(pcoords, 0)));
auto v1 = internal::lerp(static_cast<T>(values.getValue(3, c)),
static_cast<T>(values.getValue(2, c)),
static_cast<T>(component(pcoords, 0)));
auto v = internal::lerp(v0, v1, static_cast<T>(component(pcoords, 1)));
component(result, c) = static_cast<ComponentType<Result>>(v);
}
return ErrorCode::SUCCESS;
}
namespace internal
{
template <typename Values, typename CoordType, typename Result>
VTKC_EXEC inline void parametricDerivative(
Quad, const Values& values, IdComponent comp, const CoordType& pcoords, Result&& result) noexcept
{
using T = internal::ClosestFloatType<typename Values::ValueType>;
T p0 = static_cast<T>(component(pcoords, 0));
T p1 = static_cast<T>(component(pcoords, 1));
T rm = T{1} - p0;
T sm = T{1} - p1;
T dr = (static_cast<T>(values.getValue(0, comp)) * -sm) +
(static_cast<T>(values.getValue(1, comp)) * sm) +
(static_cast<T>(values.getValue(2, comp)) * p1) +
(static_cast<T>(values.getValue(3, comp)) * -p1);
T ds = (static_cast<T>(values.getValue(0, comp)) * -rm) +
(static_cast<T>(values.getValue(1, comp)) * -p0) +
(static_cast<T>(values.getValue(2, comp)) * p0) +
(static_cast<T>(values.getValue(3, comp)) * rm);
component(result, 0) = static_cast<ComponentType<Result>>(dr);
component(result, 1) = static_cast<ComponentType<Result>>(ds);
}
} // internal
template <typename Points, typename Values, typename CoordType, typename Result>
VTKC_EXEC inline vtkc::ErrorCode derivative(
Quad,
const Points& points,
const Values& values,
const CoordType& pcoords,
Result&& dx,
Result&& dy,
Result&& dz) noexcept
{
return internal::derivative2D(Quad{},
points,
values,
pcoords,
std::forward<Result>(dx),
std::forward<Result>(dy),
std::forward<Result>(dz));
}
template <typename Points, typename PCoordType, typename WCoordType>
VTKC_EXEC inline vtkc::ErrorCode parametricToWorld(
Quad,
const Points& points,
const PCoordType& pcoords,
WCoordType&& wcoords) noexcept
{
return interpolate(Quad{}, points, pcoords, std::forward<WCoordType>(wcoords));
}
template <typename Points, typename WCoordType, typename PCoordType>
VTKC_EXEC inline vtkc::ErrorCode worldToParametric(
Quad,
const Points& points,
const WCoordType& wcoords,
PCoordType&& pcoords) noexcept
{
return internal::worldToParametric2D(Quad{}, points, wcoords, std::forward<PCoordType>(pcoords));
}
} //namespace vtkc
#endif //vtk_c_Quad_h

142
vtkm/thirdparty/vtkc/vtkmvtkc/vtkc/Shapes.h vendored Normal file
View File

@ -0,0 +1,142 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.md 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_c_Shapes_h
#define vtk_c_Shapes_h
#include <vtkc/internal/Config.h>
#include <cstdint>
namespace vtkc
{
enum ShapeId : IdShape
{
// Linear cells
EMPTY = 0,
VERTEX = 1,
//POLY_VERTEX = 2,
LINE = 3,
//POLY_LINE = 4,
TRIANGLE = 5,
//TRIANGLE_STRIP = 6,
POLYGON = 7,
PIXEL = 8,
QUAD = 9,
TETRA = 10,
VOXEL = 11,
HEXAHEDRON = 12,
WEDGE = 13,
PYRAMID = 14,
NUMBER_OF_CELL_SHAPES
};
class Cell
{
public:
constexpr VTKC_EXEC Cell() : Shape(ShapeId::EMPTY), NumberOfPoints(0) {}
constexpr VTKC_EXEC Cell(IdShape shape, IdComponent numberOfPoints)
: Shape(shape), NumberOfPoints(numberOfPoints)
{
}
constexpr VTKC_EXEC IdShape shape() const noexcept { return this->Shape; }
constexpr VTKC_EXEC IdComponent numberOfPoints() const noexcept { return this->NumberOfPoints; }
protected:
IdShape Shape;
IdComponent NumberOfPoints;
};
/// \brief Check if a shape id is valid
/// \param[in] shapeId The id to check.
/// \return true if id is a valid shape id.
///
constexpr inline VTKC_EXEC bool isValidShape(IdShape shapeId)
{
return (shapeId >= ShapeId::EMPTY) && (shapeId < ShapeId::NUMBER_OF_CELL_SHAPES);
}
/// \brief Returns the dimensionality of a cell
/// \param[in] shapeId The shape id of the cell.
/// \return The dimensionality of the cell, -1 for invalid shapes.
///
inline VTKC_EXEC int dimension(IdShape shapeId)
{
switch (shapeId)
{
case VERTEX:
return 0;
case LINE:
return 1;
case TRIANGLE:
case POLYGON:
case PIXEL:
case QUAD:
return 2;
case TETRA:
case VOXEL:
case HEXAHEDRON:
case WEDGE:
case PYRAMID:
return 3;
case EMPTY:
default:
return -1;
}
}
/// \brief Returns the dimensionality of a cell
/// \param[in] cell The cell.
/// \return The dimensionality of the cell, -1 for invalid shapes.
///
inline VTKC_EXEC int dimension(Cell cell)
{
return dimension(cell.shape());
}
// forward declare cell tags
class Vertex;
class Line;
class Triangle;
class Polygon;
class Pixel;
class Quad;
class Tetra;
class Voxel;
class Hexahedron;
class Wedge;
class Pyramid;
} //namespace vtkc
#define vtkcGenericCellShapeMacroCase(cellId, cell, call) \
case cellId: \
{ \
using CellTag = cell; \
call; \
} \
break
#define vtkcGenericCellShapeMacro(call) \
vtkcGenericCellShapeMacroCase(vtkc::ShapeId::VERTEX, vtkc::Vertex, call); \
vtkcGenericCellShapeMacroCase(vtkc::ShapeId::LINE, vtkc::Line, call); \
vtkcGenericCellShapeMacroCase(vtkc::ShapeId::TRIANGLE, vtkc::Triangle, call); \
vtkcGenericCellShapeMacroCase(vtkc::ShapeId::POLYGON, vtkc::Polygon, call); \
vtkcGenericCellShapeMacroCase(vtkc::ShapeId::PIXEL, vtkc::Pixel, call); \
vtkcGenericCellShapeMacroCase(vtkc::ShapeId::QUAD, vtkc::Quad, call); \
vtkcGenericCellShapeMacroCase(vtkc::ShapeId::TETRA, vtkc::Tetra, call); \
vtkcGenericCellShapeMacroCase(vtkc::ShapeId::VOXEL, vtkc::Voxel, call); \
vtkcGenericCellShapeMacroCase(vtkc::ShapeId::HEXAHEDRON, vtkc::Hexahedron, call); \
vtkcGenericCellShapeMacroCase(vtkc::ShapeId::WEDGE, vtkc::Wedge, call); \
vtkcGenericCellShapeMacroCase(vtkc::ShapeId::PYRAMID, vtkc::Pyramid, call)
#endif //vtk_c_Shapes_h

208
vtkm/thirdparty/vtkc/vtkmvtkc/vtkc/Tetra.h vendored Normal file
View File

@ -0,0 +1,208 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.md 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_c_Tetra_h
#define vtk_c_Tetra_h
#include <vtkc/ErrorCode.h>
#include <vtkc/Shapes.h>
#include <vtkc/internal/Common.h>
namespace vtkc
{
class Tetra : public Cell
{
public:
constexpr VTKC_EXEC Tetra() : Cell(ShapeId::TETRA, 4) {}
constexpr VTKC_EXEC explicit Tetra(const Cell& cell) : Cell(cell) {}
};
VTKC_EXEC inline vtkc::ErrorCode validate(Tetra tag) noexcept
{
if (tag.shape() != ShapeId::TETRA)
{
return ErrorCode::WRONG_SHAPE_ID_FOR_TAG_TYPE;
}
if (tag.numberOfPoints() != 4)
{
return ErrorCode::INVALID_NUMBER_OF_POINTS;
}
return ErrorCode::SUCCESS;
}
template<typename CoordType>
VTKC_EXEC inline vtkc::ErrorCode parametricCenter(Tetra, CoordType&& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
component(pcoords, 0) = 0.25f;
component(pcoords, 1) = 0.25f;
component(pcoords, 2) = 0.25f;
return ErrorCode::SUCCESS;
}
template<typename CoordType>
VTKC_EXEC inline vtkc::ErrorCode parametricPoint(
Tetra, IdComponent pointId, CoordType&& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
switch (pointId)
{
case 0:
component(pcoords, 0) = 0.0f;
component(pcoords, 1) = 0.0f;
component(pcoords, 2) = 0.0f;
break;
case 1:
component(pcoords, 0) = 1.0f;
component(pcoords, 1) = 0.0f;
component(pcoords, 2) = 0.0f;
break;
case 2:
component(pcoords, 0) = 0.0f;
component(pcoords, 1) = 1.0f;
component(pcoords, 2) = 0.0f;
break;
case 3:
component(pcoords, 0) = 0.0f;
component(pcoords, 1) = 0.0f;
component(pcoords, 2) = 1.0f;
break;
default:
return ErrorCode::INVALID_POINT_ID;
}
return ErrorCode::SUCCESS;
}
template<typename CoordType>
VTKC_EXEC inline ComponentType<CoordType> parametricDistance(Tetra, const CoordType& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
ComponentType<CoordType> weights[4];
weights[0] = ComponentType<CoordType>{1} - component(pcoords, 0) - component(pcoords, 1) - component(pcoords, 2);
weights[1] = component(pcoords, 0);
weights[2] = component(pcoords, 1);
weights[3] = component(pcoords, 2);
interpolationWeights(Tetra{}, pcoords, weights);
return internal::findParametricDistance(weights, 4);
}
template<typename CoordType>
VTKC_EXEC inline bool cellInside(Tetra, const CoordType& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
using T = ComponentType<CoordType>;
return component(pcoords, 0) >= T{0} &&
component(pcoords, 1) >= T{0} &&
component(pcoords, 2) >= T{0} &&
(component(pcoords, 0) + component(pcoords, 1) + component(pcoords, 2)) <= T{1};
}
template <typename CoordType, typename T>
VTKC_EXEC inline void interpolationWeights(Tetra, const CoordType& pcoords, T weights[4]) noexcept
{
weights[1] = static_cast<T>(component(pcoords, 0));
weights[2] = static_cast<T>(component(pcoords, 1));
weights[3] = static_cast<T>(component(pcoords, 2));
weights[0] = T(1) - weights[1] - weights[2] - weights[3];
}
template <typename Values, typename CoordType, typename Result>
VTKC_EXEC inline vtkc::ErrorCode interpolate(
Tetra,
const Values& values,
const CoordType& pcoords,
Result&& result) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
using T = internal::ClosestFloatType<typename Values::ValueType>;
auto w0 = T(1) - static_cast<T>(component(pcoords, 0) + component(pcoords, 1) + component(pcoords, 2));
auto w1 = static_cast<T>(component(pcoords, 0));
auto w2 = static_cast<T>(component(pcoords, 1));
auto w3 = static_cast<T>(component(pcoords, 2));
for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
{
auto v = static_cast<T>(values.getValue(0, c)) * w0 +
static_cast<T>(values.getValue(1, c)) * w1 +
static_cast<T>(values.getValue(2, c)) * w2 +
static_cast<T>(values.getValue(3, c)) * w3;
component(result, c) = static_cast<ComponentType<Result>>(v);
}
return ErrorCode::SUCCESS;
}
namespace internal
{
template <typename Values, typename CoordType, typename Result>
VTKC_EXEC inline void parametricDerivative(
Tetra, const Values& values, IdComponent comp, const CoordType&, Result&& result) noexcept
{
component(result, 0) = static_cast<ComponentType<Result>>(values.getValue(1, comp) -
values.getValue(0, comp));
component(result, 1) = static_cast<ComponentType<Result>>(values.getValue(2, comp) -
values.getValue(0, comp));
component(result, 2) = static_cast<ComponentType<Result>>(values.getValue(3, comp) -
values.getValue(0, comp));
}
} // internal
template <typename Points, typename Values, typename CoordType, typename Result>
VTKC_EXEC inline vtkc::ErrorCode derivative(
Tetra,
const Points& points,
const Values& values,
const CoordType& pcoords,
Result&& dx,
Result&& dy,
Result&& dz) noexcept
{
return internal::derivative3D(Tetra{},
points,
values,
pcoords,
std::forward<Result>(dx),
std::forward<Result>(dy),
std::forward<Result>(dz));
}
template <typename Points, typename PCoordType, typename WCoordType>
VTKC_EXEC inline vtkc::ErrorCode parametricToWorld(
Tetra,
const Points& points,
const PCoordType& pcoords,
WCoordType&& wcoords) noexcept
{
return interpolate(Tetra{}, points, pcoords, std::forward<WCoordType>(wcoords));
}
template <typename Points, typename WCoordType, typename PCoordType>
VTKC_EXEC inline vtkc::ErrorCode worldToParametric(
Tetra,
const Points& points,
const WCoordType& wcoords,
PCoordType&& pcoords) noexcept
{
return internal::worldToParametric3D(Tetra{}, points, wcoords, std::forward<PCoordType>(pcoords));
}
} // vtkc
#endif // vtk_c_Tetra_h

271
vtkm/thirdparty/vtkc/vtkmvtkc/vtkc/Triangle.h vendored Normal file
View File

@ -0,0 +1,271 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.md 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_c_Triangle_h
#define vtk_c_Triangle_h
#include <vtkc/ErrorCode.h>
#include <vtkc/Shapes.h>
#include <vtkc/internal/Common.h>
namespace vtkc
{
class Triangle : public Cell
{
public:
constexpr VTKC_EXEC Triangle() : Cell(ShapeId::TRIANGLE, 3) {}
constexpr VTKC_EXEC explicit Triangle(const Cell& cell) : Cell(cell) {}
};
VTKC_EXEC inline vtkc::ErrorCode validate(Triangle tag) noexcept
{
if (tag.shape() != ShapeId::TRIANGLE)
{
return ErrorCode::WRONG_SHAPE_ID_FOR_TAG_TYPE;
}
if (tag.numberOfPoints() != 3)
{
return ErrorCode::INVALID_NUMBER_OF_POINTS;
}
return ErrorCode::SUCCESS;
}
template<typename CoordType>
VTKC_EXEC inline vtkc::ErrorCode parametricCenter(Triangle, CoordType&& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
using T = ComponentType<CoordType>;
component(pcoords, 0) = T(1)/T(3);
component(pcoords, 1) = T(1)/T(3);
component(pcoords, 2) = T(0);
return ErrorCode::SUCCESS;
}
template<typename CoordType>
VTKC_EXEC inline vtkc::ErrorCode parametricPoint(
Triangle, IdComponent pointId, CoordType&& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
component(pcoords, 2) = 0.0f;
switch (pointId)
{
case 0:
component(pcoords, 0) = 0.0f;
component(pcoords, 1) = 0.0f;
break;
case 1:
component(pcoords, 0) = 1.0f;
component(pcoords, 1) = 0.0f;
break;
case 2:
component(pcoords, 0) = 0.0f;
component(pcoords, 1) = 1.0f;
break;
default:
return ErrorCode::INVALID_POINT_ID;
}
return ErrorCode::SUCCESS;
}
template<typename CoordType>
VTKC_EXEC inline ComponentType<CoordType> parametricDistance(Triangle, const CoordType& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
ComponentType<CoordType> weights[3];
weights[0] = ComponentType<CoordType>{1} - component(pcoords, 0) - component(pcoords, 1);
weights[1] = component(pcoords, 0);
weights[2] = component(pcoords, 1);
return internal::findParametricDistance(weights, 3);
}
template<typename CoordType>
VTKC_EXEC inline bool cellInside(Triangle, const CoordType& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
using T = ComponentType<CoordType>;
return component(pcoords, 0) >= T{0} && component(pcoords, 1) >= T{0} &&
(component(pcoords, 0) + component(pcoords, 1)) <= T{1};
}
template <typename Values, typename CoordType, typename Result>
VTKC_EXEC inline vtkc::ErrorCode interpolate(
Triangle,
const Values& values,
const CoordType& pcoords,
Result&& result) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
using T = internal::ClosestFloatType<typename Values::ValueType>;
auto w0 = T(1) - static_cast<T>(component(pcoords, 0) + component(pcoords, 1));
auto w1 = static_cast<T>(component(pcoords, 0));
auto w2 = static_cast<T>(component(pcoords, 1));
for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
{
auto v = static_cast<T>(values.getValue(0, c)) * w0 +
static_cast<T>(values.getValue(1, c)) * w1 +
static_cast<T>(values.getValue(2, c)) * w2;
component(result, c) = static_cast<ComponentType<Result>>(v);
}
return ErrorCode::SUCCESS;
}
namespace internal
{
template <typename Values, typename CoordType, typename Result>
VTKC_EXEC inline void parametricDerivative(
Triangle, const Values& values, IdComponent comp, const CoordType&, Result&& result) noexcept
{
component(result, 0) = static_cast<ComponentType<Result>>(values.getValue(1, comp) -
values.getValue(0, comp));
component(result, 1) = static_cast<ComponentType<Result>>(values.getValue(2, comp) -
values.getValue(0, comp));
}
} // internal
template <typename Points, typename Values, typename CoordType, typename Result>
VTKC_EXEC inline vtkc::ErrorCode derivative(
Triangle,
const Points& points,
const Values& values,
const CoordType& pcoords,
Result&& dx,
Result&& dy,
Result&& dz) noexcept
{
return internal::derivative2D(Triangle{},
points,
values,
pcoords,
std::forward<Result>(dx),
std::forward<Result>(dy),
std::forward<Result>(dz));
}
template <typename Points, typename PCoordType, typename WCoordType>
VTKC_EXEC inline vtkc::ErrorCode parametricToWorld(
Triangle,
const Points& points,
const PCoordType& pcoords,
WCoordType&& wcoords) noexcept
{
return interpolate(Triangle{}, points, pcoords, std::forward<WCoordType>(wcoords));
}
//-----------------------------------------------------------------------------
// The following implementation is lifted from VTK-m
//-----------------------------------------------------------------------------
// 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.
//
template <typename Points, typename WCoordType, typename PCoordType>
VTKC_EXEC inline vtkc::ErrorCode worldToParametric(
Triangle,
const Points& points,
const WCoordType& wcoords,
PCoordType&& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(PCoordType);
using TIn = typename Points::ValueType;
using TOut = ComponentType<PCoordType>;
internal::Vector<TIn, 3> pts[3];
for (int i = 0; i < 3; ++i)
{
points.getTuple(i, pts[i]);
}
internal::Vector<TIn, 3> wc(component(wcoords, 0), component(wcoords, 1), component(wcoords, 2));
auto triangleNormal = internal::cross(pts[1] - pts[0], pts[2] - pts[0]);
for (IdComponent i = 0; i < 2; ++i)
{
auto& p0 = pts[0];
auto& p1 = pts[i + 1];
auto& p2 = pts[2 - i];
auto planeNormal = internal::cross(triangleNormal, p2 - p0);
component(pcoords, i) = static_cast<TOut>(internal::dot(wc - p0, planeNormal)) /
static_cast<TOut>(internal::dot(p1 - p0, planeNormal));
}
return ErrorCode::SUCCESS;
}
} //namespace vtkc
#endif //vtk_c_Triangle_h

131
vtkm/thirdparty/vtkc/vtkmvtkc/vtkc/Vertex.h vendored Normal file
View File

@ -0,0 +1,131 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.md 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_c_Vertex_h
#define vtk_c_Vertex_h
#include <vtkc/ErrorCode.h>
#include <vtkc/Shapes.h>
namespace vtkc
{
class Vertex : public Cell
{
public:
constexpr VTKC_EXEC Vertex() : Cell(ShapeId::VERTEX, 1) {}
constexpr VTKC_EXEC explicit Vertex(const Cell& cell) : Cell(cell) {}
};
VTKC_EXEC inline vtkc::ErrorCode validate(Vertex tag) noexcept
{
if (tag.shape() != ShapeId::VERTEX)
{
return ErrorCode::WRONG_SHAPE_ID_FOR_TAG_TYPE;
}
if (tag.numberOfPoints() != 1)
{
return ErrorCode::INVALID_NUMBER_OF_POINTS;
}
return ErrorCode::SUCCESS;
}
template<typename CoordType>
VTKC_EXEC inline vtkc::ErrorCode parametricCenter(Vertex, CoordType&&) noexcept
{
return ErrorCode::SUCCESS;
}
template<typename CoordType>
VTKC_EXEC inline vtkc::ErrorCode parametricPoint(Vertex, IdComponent pointId, CoordType&&) noexcept
{
if (pointId == 0)
{
return ErrorCode::SUCCESS;
}
else
{
return ErrorCode::INVALID_POINT_ID;
}
}
template<typename CoordType>
VTKC_EXEC inline ComponentType<CoordType> parametricDistance(Vertex, const CoordType&) noexcept
{
return ComponentType<CoordType>{1};
}
template<typename CoordType>
VTKC_EXEC inline bool cellInside(Vertex, const CoordType&) noexcept
{
return false;
}
template <typename Values, typename CoordType, typename Result>
VTKC_EXEC inline vtkc::ErrorCode interpolate(
Vertex,
const Values& values,
const CoordType&,
Result&& result) noexcept
{
using T = ComponentType<Result>;
for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
{
component(result, c) = static_cast<T>(values.getValue(0, c));
}
return ErrorCode::SUCCESS;
}
template <typename Points, typename Values, typename CoordType, typename Result>
VTKC_EXEC inline vtkc::ErrorCode derivative(
Vertex,
const Points&,
const Values& values,
const CoordType&,
Result&& dx,
Result&& dy,
Result&& dz) noexcept
{
using T = ComponentType<Result>;
for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
{
component(dx, c) = component(dy, c) = component(dz, c) = T{0};
}
return ErrorCode::SUCCESS;
}
template <typename Points, typename PCoordType, typename WCoordType>
VTKC_EXEC inline vtkc::ErrorCode parametricToWorld(
Vertex,
const Points& points,
const PCoordType&,
WCoordType&& wcoords) noexcept
{
using T = ComponentType<WCoordType>;
component(wcoords, 0) = static_cast<T>(points.getValue(0, 0));
component(wcoords, 1) = static_cast<T>(points.getValue(0, 1));
component(wcoords, 2) = static_cast<T>(points.getValue(0, 2));
return ErrorCode::SUCCESS;
}
template <typename Points, typename WCoordType, typename PCoordType>
VTKC_EXEC inline vtkc::ErrorCode worldToParametric(
Vertex,
const Points&,
const WCoordType&,
PCoordType&&) noexcept
{
return ErrorCode::SUCCESS;
}
} //namespace vtkc
#endif //vtk_c_Vertex_h

136
vtkm/thirdparty/vtkc/vtkmvtkc/vtkc/Voxel.h vendored Normal file
View File

@ -0,0 +1,136 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.md 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_c_Voxel_h
#define vtk_c_Voxel_h
#include <vtkc/ErrorCode.h>
#include <vtkc/Hexahedron.h>
#include <vtkc/internal/Common.h>
namespace vtkc
{
class Voxel : public Hexahedron
{
public:
constexpr VTKC_EXEC Voxel() : Hexahedron(Cell(ShapeId::VOXEL, 8)) {}
constexpr VTKC_EXEC explicit Voxel(const Cell& cell) : Hexahedron(cell) {}
};
namespace internal
{
template <typename Points, typename T>
VTKC_EXEC inline int getVoxelSpacing(const Points& points, T spacing[3])
{
int zeros = 0;
for (int i = 0; i < 3; ++i)
{
spacing[i] = static_cast<T>(points.getValue(6, i) - points.getValue(0, i));
if (spacing[i] == T(0))
{
zeros |= 1 << i;
}
}
return zeros;
}
} // internal
template <typename Points, typename Values, typename CoordType, typename Result>
VTKC_EXEC inline vtkc::ErrorCode derivative(
Voxel,
const Points& points,
const Values& values,
const CoordType& pcoords,
Result&& dx,
Result&& dy,
Result&& dz) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
using ProcessingType = internal::ClosestFloatType<typename Values::ValueType>;
using ResultCompType = ComponentType<Result>;
ProcessingType spacing[3];
if (internal::getVoxelSpacing(points, spacing) != 0)
{
return ErrorCode::DEGENERATE_CELL_DETECTED;
}
for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
{
ProcessingType dvdp[3];
internal::parametricDerivative(Hexahedron{}, values, c, pcoords, dvdp);
component(dx, c) = static_cast<ResultCompType>(dvdp[0] / spacing[0]);
component(dy, c) = static_cast<ResultCompType>(dvdp[1] / spacing[1]);
component(dz, c) = static_cast<ResultCompType>(dvdp[2] / spacing[2]);
}
return ErrorCode::SUCCESS;
}
template <typename Points, typename PCoordType, typename WCoordType>
VTKC_EXEC vtkc::ErrorCode parametricToWorld(
Voxel,
const Points& points,
const PCoordType& pcoords,
WCoordType&& wcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(PCoordType);
using T = ComponentType<WCoordType>;
T spacing[3];
if (internal::getVoxelSpacing(points, spacing) != 0)
{
return ErrorCode::DEGENERATE_CELL_DETECTED;
}
component(wcoords, 0) = points.getValue(0, 0) +
(spacing[0] * static_cast<T>(component(pcoords, 0)));
component(wcoords, 1) = points.getValue(0, 1) +
(spacing[1] * static_cast<T>(component(pcoords, 1)));
component(wcoords, 2) = points.getValue(0, 2) +
(spacing[2] * static_cast<T>(component(pcoords, 2)));
return ErrorCode::SUCCESS;
}
template <typename Points, typename WCoordType, typename PCoordType>
VTKC_EXEC vtkc::ErrorCode worldToParametric(
Voxel,
const Points& points,
const WCoordType& wcoords,
PCoordType&& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(PCoordType);
using T = ComponentType<PCoordType>;
T spacing[3];
if (internal::getVoxelSpacing(points, spacing) != 0)
{
return ErrorCode::DEGENERATE_CELL_DETECTED;
}
component(pcoords, 0) = static_cast<T>(component(wcoords, 0) - points.getValue(0, 0)) /
spacing[0];
component(pcoords, 1) = static_cast<T>(component(wcoords, 1) - points.getValue(0, 1)) /
spacing[1];
component(pcoords, 2) = static_cast<T>(component(wcoords, 2) - points.getValue(0, 2)) /
spacing[2];
return ErrorCode::SUCCESS;
}
} //namespace vtkc
#endif //vtk_c_Voxel_h

228
vtkm/thirdparty/vtkc/vtkmvtkc/vtkc/Wedge.h vendored Normal file
View File

@ -0,0 +1,228 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.md 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_c_Wedge_h
#define vtk_c_Wedge_h
#include <vtkc/ErrorCode.h>
#include <vtkc/Shapes.h>
#include <vtkc/internal/Common.h>
namespace vtkc
{
class Wedge : public Cell
{
public:
constexpr VTKC_EXEC Wedge() : Cell(ShapeId::WEDGE, 6) {}
constexpr VTKC_EXEC explicit Wedge(const Cell& cell) : Cell(cell) {}
};
VTKC_EXEC inline vtkc::ErrorCode validate(Wedge tag) noexcept
{
if (tag.shape() != ShapeId::WEDGE)
{
return ErrorCode::WRONG_SHAPE_ID_FOR_TAG_TYPE;
}
if (tag.numberOfPoints() != 6)
{
return ErrorCode::INVALID_NUMBER_OF_POINTS;
}
return ErrorCode::SUCCESS;
}
template<typename CoordType>
VTKC_EXEC inline vtkc::ErrorCode parametricCenter(Wedge, CoordType&& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
using T = ComponentType<CoordType>;
component(pcoords, 0) = T(1)/T(3);
component(pcoords, 1) = T(1)/T(3);
component(pcoords, 2) = 0.5f;
return ErrorCode::SUCCESS;
}
template<typename CoordType>
VTKC_EXEC inline vtkc::ErrorCode parametricPoint(
Wedge, IdComponent pointId, CoordType&& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
switch (pointId)
{
case 0:
component(pcoords, 0) = 0.0f;
component(pcoords, 1) = 0.0f;
component(pcoords, 2) = 0.0f;
break;
case 1:
component(pcoords, 0) = 1.0f;
component(pcoords, 1) = 0.0f;
component(pcoords, 2) = 0.0f;
break;
case 2:
component(pcoords, 0) = 0.0f;
component(pcoords, 1) = 1.0f;
component(pcoords, 2) = 0.0f;
break;
case 3:
component(pcoords, 0) = 0.0f;
component(pcoords, 1) = 0.0f;
component(pcoords, 2) = 1.0f;
break;
case 4:
component(pcoords, 0) = 1.0f;
component(pcoords, 1) = 0.0f;
component(pcoords, 2) = 1.0f;
break;
case 5:
component(pcoords, 0) = 0.0f;
component(pcoords, 1) = 1.0f;
component(pcoords, 2) = 1.0f;
break;
default:
return ErrorCode::INVALID_POINT_ID;
}
return ErrorCode::SUCCESS;
}
template<typename CoordType>
VTKC_EXEC inline ComponentType<CoordType> parametricDistance(Wedge, const CoordType& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
return internal::findParametricDistance(pcoords, 3);
}
template<typename CoordType>
VTKC_EXEC inline bool cellInside(Wedge, const CoordType& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
using T = ComponentType<CoordType>;
return component(pcoords, 0) >= T{0} &&
component(pcoords, 1) >= T{0} &&
component(pcoords, 2) >= T{0} &&
(component(pcoords, 0) + component(pcoords, 1)) <= T{1} &&
component(pcoords, 2) <= T{1};
}
template <typename Values, typename CoordType, typename Result>
VTKC_EXEC vtkc::ErrorCode interpolate(
Wedge,
const Values& values,
const CoordType& pcoords,
Result&& result) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
using T = internal::ClosestFloatType<typename Values::ValueType>;
auto p0 = static_cast<T>(component(pcoords, 0));
auto p1 = static_cast<T>(component(pcoords, 1));
auto p2 = static_cast<T>(component(pcoords, 2));
auto sm = T{1} - p0 - p1;
for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
{
T t1 = static_cast<T>(values.getValue(0, c)) * sm +
static_cast<T>(values.getValue(1, c)) * p0 +
static_cast<T>(values.getValue(2, c)) * p1;
T t2 = static_cast<T>(values.getValue(3, c)) * sm +
static_cast<T>(values.getValue(4, c)) * p0 +
static_cast<T>(values.getValue(5, c)) * p1;
component(result, c) = static_cast<ComponentType<Result>>(internal::lerp(t1, t2, p2));
}
return ErrorCode::SUCCESS;
}
namespace internal
{
template <typename Values, typename CoordType, typename Result>
VTKC_EXEC inline void parametricDerivative(
Wedge, const Values& values, IdComponent comp, const CoordType& pcoords, Result&& result) noexcept
{
using T = internal::ClosestFloatType<typename Values::ValueType>;
auto p0 = static_cast<T>(component(pcoords, 0));
auto p1 = static_cast<T>(component(pcoords, 1));
auto p2 = static_cast<T>(component(pcoords, 2));
auto rm = T{1} - p2;
auto sm = T{1} - p0 - p1;
T dr = (static_cast<T>(values.getValue(0, comp)) * -rm) +
(static_cast<T>(values.getValue(1, comp)) * rm) +
(static_cast<T>(values.getValue(3, comp)) * -p2) +
(static_cast<T>(values.getValue(4, comp)) * p2);
T ds = (static_cast<T>(values.getValue(0, comp)) * -rm) +
(static_cast<T>(values.getValue(2, comp)) * rm) +
(static_cast<T>(values.getValue(3, comp)) * -p2) +
(static_cast<T>(values.getValue(5, comp)) * p2);
T dt = (static_cast<T>(values.getValue(0, comp)) * -sm) +
(static_cast<T>(values.getValue(1, comp)) * -p0) +
(static_cast<T>(values.getValue(2, comp)) * -p1) +
(static_cast<T>(values.getValue(3, comp)) * sm) +
(static_cast<T>(values.getValue(4, comp)) * p0) +
(static_cast<T>(values.getValue(5, comp)) * p1);
component(result, 0) = static_cast<ComponentType<Result>>(dr);
component(result, 1) = static_cast<ComponentType<Result>>(ds);
component(result, 2) = static_cast<ComponentType<Result>>(dt);
}
} // internal
template <typename Points, typename Values, typename CoordType, typename Result>
VTKC_EXEC inline vtkc::ErrorCode derivative(
Wedge,
const Points& points,
const Values& values,
const CoordType& pcoords,
Result&& dx,
Result&& dy,
Result&& dz) noexcept
{
return internal::derivative3D(Wedge{},
points,
values,
pcoords,
std::forward<Result>(dx),
std::forward<Result>(dy),
std::forward<Result>(dz));
}
template <typename Points, typename PCoordType, typename WCoordType>
VTKC_EXEC inline vtkc::ErrorCode parametricToWorld(
Wedge,
const Points& points,
const PCoordType& pcoords,
WCoordType&& wcoords) noexcept
{
return interpolate(Wedge{}, points, pcoords, std::forward<WCoordType>(wcoords));
}
template <typename Points, typename WCoordType, typename PCoordType>
VTKC_EXEC inline vtkc::ErrorCode worldToParametric(
Wedge,
const Points& points,
const WCoordType& wcoords,
PCoordType&& pcoords) noexcept
{
return internal::worldToParametric3D(Wedge{}, points, wcoords, std::forward<PCoordType>(pcoords));
}
} // vtkc
#endif // vtk_c_Wedge_h

310
vtkm/thirdparty/vtkc/vtkmvtkc/vtkc/internal/Common.h vendored Normal file
View File

@ -0,0 +1,310 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.md 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_c_internal_Common_h
#define vtk_c_internal_Common_h
#include <vtkc/FieldAccessor.h>
#include <vtkc/internal/Math.h>
#include <cstdint>
#include <type_traits>
#define VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(PCoordType) \
static_assert(std::is_floating_point<ComponentType<PCoordType>>::value, \
"parametric coordinates should be of floating point type");
namespace vtkc
{
namespace internal
{
///=========================================================================
template<typename CoordType>
VTKC_EXEC inline ComponentType<CoordType> findParametricDistance(
const CoordType& pvals, IdComponent numVals) noexcept
{
using T = ComponentType<CoordType>;
T pDistMax = 0.0f;
for (IdComponent i = 0; i < numVals; ++i)
{
ComponentType<CoordType> pDist = 0.0f;
if (component(pvals, i) < T(0))
{
pDist = -1.0f * component(pvals, i);
}
else if (component(pvals, i) > T(1))
{
pDist = component(pvals, i) - 1.0f;
}
if (pDist > pDistMax)
{
pDistMax = pDist;
}
}
return pDistMax;
}
///=========================================================================
/// Forward declaration
#define FORWARD_DECLAR_PARAMETRIC_DERIVATIVE(tag) \
template <typename Values, typename CoordType, typename Result> \
VTKC_EXEC inline void parametricDerivative( \
vtkc::tag, const Values& values, IdComponent comp, const CoordType&, Result&& result) noexcept
FORWARD_DECLAR_PARAMETRIC_DERIVATIVE(Triangle);
FORWARD_DECLAR_PARAMETRIC_DERIVATIVE(Quad);
FORWARD_DECLAR_PARAMETRIC_DERIVATIVE(Tetra);
FORWARD_DECLAR_PARAMETRIC_DERIVATIVE(Hexahedron);
FORWARD_DECLAR_PARAMETRIC_DERIVATIVE(Wedge);
FORWARD_DECLAR_PARAMETRIC_DERIVATIVE(Pyramid);
#undef FORWARD_DECLAR_PARAMETRIC_DERIVATIVE
///=========================================================================
template <typename T>
class Space2D
{
public:
explicit VTKC_EXEC Space2D(const Vector<T, 3>& origin, const Vector<T, 3>& p1, const Vector<T, 3>& p2) noexcept
{
this->Origin = origin;
this->XAxis = p1 - origin;
auto normal = internal::cross(this->XAxis, p2 - origin);
this->YAxis = internal::cross(normal, this->XAxis);
internal::normalize(this->XAxis);
internal::normalize(this->YAxis);
}
VTKC_EXEC Vector<T, 2> to2DPoint(Vector<T, 3> pt) const noexcept
{
pt -= this->Origin;
return Vector<T, 2>{ internal::dot(pt, this->XAxis), internal::dot(pt, this->YAxis) };
}
VTKC_EXEC Vector<T, 3> to3DVec(const Vector<T, 2>& vec) const noexcept
{
return (this->XAxis * vec[0]) + (this->YAxis * vec[1]);
}
private:
Vector<T, 3> Origin;
Vector<T, 3> XAxis, YAxis;
};
template <typename CellTag, typename Points, typename PCoords, typename T>
VTKC_EXEC inline void jacobian2D(
CellTag tag, const Points& points, const PCoords& pcoords, Matrix<T, 2, 2>& jacobian) noexcept
{
T pd[2];
parametricDerivative(tag, points, 0, pcoords, pd);
jacobian(0, 0) = pd[0];
jacobian(0, 1) = pd[1];
parametricDerivative(tag, points, 1, pcoords, pd);
jacobian(1, 0) = pd[0];
jacobian(1, 1) = pd[1];
}
template <typename CellTag, typename Points, typename Values, typename CoordType, typename Result>
VTKC_EXEC inline vtkc::ErrorCode derivative2D(
CellTag tag,
const Points& points,
const Values& values,
const CoordType& pcoords,
Result&& dx,
Result&& dy,
Result&& dz) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
using ProcessingType = internal::ClosestFloatType<typename Values::ValueType>;
using ResultCompType = ComponentType<Result>;
constexpr IdComponent numPoints = CellTag{}.numberOfPoints();
Vector<ProcessingType, 3> pts[numPoints];
for (int i = 0; i < numPoints; ++i)
{
points.getTuple(i, pts[i]);
}
// 2-D coordinate system on the cell's plane
Space2D<ProcessingType> planeSpace(pts[0], pts[1], pts[numPoints - 1]);
Vector<ProcessingType, 2> pts2d[numPoints];
for (int i = 0; i < numPoints; ++i)
{
pts2d[i] = planeSpace.to2DPoint(pts[i]);
}
Matrix<ProcessingType, 2, 2> jacobian;
jacobian2D(tag, makeFieldAccessorNestedSOA(pts2d, 2), pcoords, jacobian);
Matrix<ProcessingType, 2, 2> invJacobian;
VTKC_RETURN_ON_ERROR(matrixInverse(jacobian, invJacobian))
for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
{
Vector<ProcessingType, 2> dvdp;
parametricDerivative(tag, values, c, pcoords, dvdp);
auto d2D = matrixMultiply(dvdp, invJacobian);
auto d3D = planeSpace.to3DVec(d2D);
component(dx, c) = static_cast<ResultCompType>(d3D[0]);
component(dy, c) = static_cast<ResultCompType>(d3D[1]);
component(dz, c) = static_cast<ResultCompType>(d3D[2]);
}
return ErrorCode::SUCCESS;
}
template <typename CellTag, typename Points, typename WCoordType, typename PCoordType>
VTKC_EXEC inline vtkc::ErrorCode worldToParametric2D(
CellTag tag,
const Points& points,
const WCoordType& wcoords,
PCoordType&& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(PCoordType);
using TIn = typename Points::ValueType;
using TOut = ComponentType<PCoordType>;
constexpr IdComponent numPoints = CellTag{}.numberOfPoints();
Vector<TIn, 3> pts[numPoints];
for (int i = 0; i < numPoints; ++i)
{
points.getTuple(i, pts[i]);
}
// 2-D coordinate system on the cell's plane
Space2D<TIn> planeSpace(pts[0], pts[1], pts[numPoints - 1]);
Vector<TIn, 2> pts2d[numPoints];
for (int i = 0; i < numPoints; ++i)
{
pts2d[i] = planeSpace.to2DPoint(pts[i]);
}
auto jacobianEvaluator = [&pts2d](const Vector<TOut, 2>& pc, Matrix<TIn, 2, 2>& jacobian) {
jacobian2D(CellTag{}, makeFieldAccessorNestedSOA(pts2d, 2), pc, jacobian);
return ErrorCode::SUCCESS;
};
auto functionEvaluator = [&points, &planeSpace](const Vector<TOut, 2>& pc, Vector<TIn, 2>& wc) {
Vector<TIn, 3> wc3(0);
VTKC_RETURN_ON_ERROR(parametricToWorld(CellTag{}, points, pc, wc3))
wc = planeSpace.to2DPoint(wc3);
return ErrorCode::SUCCESS;
};
Vector<TIn, 3> wcVec{component(wcoords, 0), component(wcoords, 1), component(wcoords, 2)};
Vector<TOut, 2> pcVec;
VTKC_RETURN_ON_ERROR(parametricCenter(tag, pcVec))
VTKC_RETURN_ON_ERROR(newtonsMethod(
jacobianEvaluator, functionEvaluator, planeSpace.to2DPoint(wcVec), pcVec))
component(pcoords, 0) = pcVec[0];
component(pcoords, 1) = pcVec[1];
return ErrorCode::SUCCESS;
}
///=========================================================================
template <typename CellTag, typename Points, typename PCoords, typename T>
VTKC_EXEC inline void jacobian3D(
CellTag tag, const Points& points, const PCoords& pcoords, Matrix<T, 3, 3>& jacobian) noexcept
{
T pd[3];
parametricDerivative(tag, points, 0, pcoords, pd);
jacobian(0, 0) = pd[0];
jacobian(0, 1) = pd[1];
jacobian(0, 2) = pd[2];
parametricDerivative(tag, points, 1, pcoords, pd);
jacobian(1, 0) = pd[0];
jacobian(1, 1) = pd[1];
jacobian(1, 2) = pd[2];
parametricDerivative(tag, points, 2, pcoords, pd);
jacobian(2, 0) = pd[0];
jacobian(2, 1) = pd[1];
jacobian(2, 2) = pd[2];
}
template <typename CellTag, typename Points, typename Values, typename CoordType, typename Result>
VTKC_EXEC inline vtkc::ErrorCode derivative3D(
CellTag tag,
const Points& points,
const Values& values,
const CoordType& pcoords,
Result&& dx,
Result&& dy,
Result&& dz) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
using ProcessingType = internal::ClosestFloatType<typename Values::ValueType>;
using ResultCompType = ComponentType<Result>;
Matrix<ProcessingType, 3, 3> jacobian;
jacobian3D(tag, points, pcoords, jacobian);
Matrix<ProcessingType, 3, 3> invJacobian;
VTKC_RETURN_ON_ERROR(matrixInverse(jacobian, invJacobian))
for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
{
Vector<ProcessingType, 3> dvdp;
parametricDerivative(tag, values, c, pcoords, dvdp);
auto deriv = matrixMultiply(dvdp, invJacobian);
component(dx, c) = static_cast<ResultCompType>(deriv[0]);
component(dy, c) = static_cast<ResultCompType>(deriv[1]);
component(dz, c) = static_cast<ResultCompType>(deriv[2]);
}
return ErrorCode::SUCCESS;
}
template <typename CellTag, typename Points, typename WCoordType, typename PCoordType>
VTKC_EXEC inline vtkc::ErrorCode worldToParametric3D(
CellTag tag,
const Points& points,
const WCoordType& wcoords,
PCoordType&& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(PCoordType);
using TIn = typename Points::ValueType;
using TOut = ComponentType<PCoordType>;
auto jacobianEvaluator = [tag, &points](const Vector<TOut, 3>& pc, Matrix<TIn, 3, 3>& jacobian) {
jacobian3D(tag, points, pc, jacobian);
return ErrorCode::SUCCESS;
};
auto functionEvaluator = [tag, &points](const Vector<TOut, 3>& pc, Vector<TIn, 3>& wc) {
return parametricToWorld(tag, points, pc, wc);
};
internal::Vector<TIn, 3> wcVec{component(wcoords, 0), component(wcoords, 1), component(wcoords, 2)};
internal::Vector<TOut, 3> pcVec;
VTKC_RETURN_ON_ERROR(parametricCenter(tag, pcVec))
VTKC_RETURN_ON_ERROR(newtonsMethod(jacobianEvaluator, functionEvaluator, wcVec, pcVec))
component(pcoords, 0) = pcVec[0];
component(pcoords, 1) = pcVec[1];
component(pcoords, 2) = pcVec[2];
return ErrorCode::SUCCESS;
}
}
} // vtkc::internal
#endif //vtk_c_internal_Common_h

38
vtkm/thirdparty/vtkc/vtkmvtkc/vtkc/internal/Config.h vendored Normal file
View File

@ -0,0 +1,38 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.md 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_c_internal_Config_h
#define vtk_c_internal_Config_h
#include <cstdint>
#include <type_traits>
#ifdef __CUDACC__
# define VTKC_EXEC __device__ __host__
#else
# define VTKC_EXEC
#endif
namespace vtkc
{
namespace internal
{
template <typename T>
using ClosestFloatType =
typename std::enable_if<std::is_arithmetic<T>::value,
typename std::conditional<sizeof(T) <= 4, float, double>::type>::type;
}
using IdShape = std::int8_t;
using IdComponent = std::int32_t;
} // namespace vtkc
#endif // vtk_c_internal_Config_h

712
vtkm/thirdparty/vtkc/vtkmvtkc/vtkc/internal/Math.h vendored Normal file
View File

@ -0,0 +1,712 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.md 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_c_internal_Math_h
#define vtk_c_internal_Math_h
#include <cmath>
#ifdef __CUDA_ARCH__
# define VTKC_MATH_CALL(func, ...) func(__VA_ARGS__)
#else
# define VTKC_MATH_CALL(func, ...) std::func(__VA_ARGS__)
#endif
namespace vtkc
{
namespace internal
{
///=========================================================================
/// Basic vector class and functions for internal use
template <typename T, IdComponent Dim>
class Vector
{
public:
#if defined(_MSC_VER) && (_MSC_VER == 1900)
// workaround vs2015 bug producing the following error:
// error C2476: ‘constexpr’ constructor does not initialize all members
VTKC_EXEC Vector() noexcept {};
#else
VTKC_EXEC
constexpr Vector() noexcept {};
#endif
VTKC_EXEC
explicit Vector(const T& val) noexcept
{
for (auto& c : Data)
{
c = val;
}
}
template <typename... Ts>
VTKC_EXEC
constexpr explicit Vector(const T& c1, const Ts&... cs) noexcept
: Data{c1, cs...}
{
static_assert(sizeof...(Ts) == (Dim - 1), "Invalid number of components passed");
}
static constexpr IdComponent getNumberOfComponents() noexcept
{
return Dim;
}
VTKC_EXEC
T& operator[](IdComponent c) noexcept
{
return this->Data[c];
}
VTKC_EXEC
constexpr const T& operator[](IdComponent c) const noexcept
{
return this->Data[c];
}
VTKC_EXEC
T* data() noexcept
{
return this->Data;
}
VTKC_EXEC
constexpr const T* data() const noexcept
{
return this->Data;
}
private:
T Data[Dim];
};
//---------------------------------------------------------------------------
template <typename T, IdComponent Dim>
VTKC_EXEC
Vector<T, Dim>& operator+=(Vector<T, Dim>& v1, const Vector<T, Dim>& v2) noexcept
{
for (IdComponent i = 0; i < Dim; ++i)
{
v1[i] += v2[i];
}
return v1;
}
template <typename T, IdComponent Dim>
VTKC_EXEC
Vector<T, Dim>& operator-=(Vector<T, Dim>& v1, const Vector<T, Dim>& v2) noexcept
{
for (IdComponent i = 0; i < Dim; ++i)
{
v1[i] -= v2[i];
}
return v1;
}
template <typename T, IdComponent Dim>
VTKC_EXEC
Vector<T, Dim>& operator*=(Vector<T, Dim>& v1, const Vector<T, Dim>& v2) noexcept
{
for (IdComponent i = 0; i < Dim; ++i)
{
v1[i] *= v2[i];
}
return v1;
}
template <typename T, IdComponent Dim>
VTKC_EXEC
Vector<T, Dim>& operator*=(Vector<T, Dim>& v, const T& s) noexcept
{
for (IdComponent i = 0; i < Dim; ++i)
{
v[i] *= s;
}
return v;
}
template <typename T, IdComponent Dim>
VTKC_EXEC
Vector<T, Dim>& operator/=(Vector<T, Dim>& v1, const Vector<T, Dim>& v2) noexcept
{
for (IdComponent i = 0; i < Dim; ++i)
{
v1[i] /= v2[i];
}
return v1;
}
template <typename T, IdComponent Dim>
VTKC_EXEC
Vector<T, Dim>& operator/=(Vector<T, Dim>& v, const T& s) noexcept
{
for (IdComponent i = 0; i < Dim; ++i)
{
v[i] /= s;
}
return v;
}
template <typename T, IdComponent Dim>
VTKC_EXEC
Vector<T, Dim> operator+(Vector<T, Dim> v1, const Vector<T, Dim>& v2) noexcept
{
return v1 += v2;
}
template <typename T, IdComponent Dim>
VTKC_EXEC
Vector<T, Dim> operator-(Vector<T, Dim> v1, const Vector<T, Dim>& v2) noexcept
{
return v1 -= v2;
}
template <typename T, IdComponent Dim>
VTKC_EXEC
Vector<T, Dim> operator*(Vector<T, Dim> v1, const Vector<T, Dim>& v2) noexcept
{
return v1 *= v2;
}
template <typename T, IdComponent Dim>
VTKC_EXEC
Vector<T, Dim> operator*(Vector<T, Dim> v, const T& s) noexcept
{
return v *= s;
}
template <typename T, IdComponent Dim>
VTKC_EXEC
Vector<T, Dim> operator*(const T& s, Vector<T, Dim> v) noexcept
{
return v *= s;
}
template <typename T, IdComponent Dim>
VTKC_EXEC
Vector<T, Dim> operator/(Vector<T, Dim> v1, const Vector<T, Dim>& v2) noexcept
{
return v1 /= v2;
}
template <typename T, IdComponent Dim>
VTKC_EXEC
Vector<T, Dim> operator/(Vector<T, Dim> v, const T& s) noexcept
{
return v /= s;
}
//---------------------------------------------------------------------------
template <typename T, IdComponent Dim>
VTKC_EXEC
T dot(const Vector<T, Dim>& v1, const Vector<T, Dim>& v2) noexcept
{
T result{};
for (IdComponent i = 0; i < Dim; ++i)
{
result += v1[i] * v2[i];
}
return result;
}
template <typename T>
VTKC_EXEC
Vector<T, 3> cross(const Vector<T, 3>& v1, const Vector<T, 3>& v2) noexcept
{
return Vector<T, 3>(v1[1] * v2[2] - v1[2] * v2[1],
v1[2] * v2[0] - v1[0] * v2[2],
v1[0] * v2[1] - v1[1] * v2[0]);
}
template <typename T, IdComponent Dim>
VTKC_EXEC
T magnitude(const Vector<T, Dim>& v)
{
return VTKC_MATH_CALL(sqrt, (dot(v, v)));
}
template <typename T, IdComponent Dim>
VTKC_EXEC
void normalize(Vector<T, Dim>& v)
{
v /= magnitude(v);
}
template <typename T, IdComponent Dim>
VTKC_EXEC
Vector<T, Dim> normal(Vector<T, Dim> v)
{
normalize(v);
return v;
}
///=========================================================================
/// Basic matrix class and functions for internal use
template <typename T, IdComponent NumberOfRows, IdComponent NumberOfColumns>
class Matrix
{
public:
#if defined(_MSC_VER) && (_MSC_VER == 1900)
// workaround vs2015 bug producing the following error:
// error C2476: ‘constexpr’ constructor does not initialize all members
VTKC_EXEC Matrix() noexcept {};
#else
VTKC_EXEC
constexpr Matrix() noexcept {};
#endif
VTKC_EXEC
T& operator()(IdComponent r, IdComponent c) noexcept
{
return this->Columns[c][r];
}
VTKC_EXEC
constexpr const T& operator()(IdComponent r, IdComponent c) const noexcept
{
return this->Columns[c][r];
}
VTKC_EXEC
static constexpr IdComponent getNumberOfRows() noexcept
{
return NumberOfRows;
}
VTKC_EXEC
static constexpr IdComponent getNumberOfColumns() noexcept
{
return NumberOfColumns;
}
VTKC_EXEC
constexpr const Vector<T, NumberOfRows>& getColumn(IdComponent c) const
{
return this->Columns[c];
}
VTKC_EXEC
Vector<T, NumberOfColumns> getRow(IdComponent r) const
{
Vector<T, NumberOfColumns> row;
for (IdComponent i = 0; i < NumberOfColumns; ++i)
{
row[i] = this->Columns[i][r];
}
return row;
}
VTKC_EXEC
Matrix& operator+=(const Matrix& m) noexcept
{
for (IdComponent i = 0; i < NumberOfColumns; ++i)
{
this->Columns[i] += m.Columns[i];
}
return *this;
}
VTKC_EXEC
Matrix& operator-=(const Matrix& m) noexcept
{
for (IdComponent i = 0; i < NumberOfColumns; ++i)
{
this->Columns[i] -= m.Columns[i];
}
return *this;
}
VTKC_EXEC
Matrix& operator*=(const T& s) noexcept
{
for (auto& c : this->Columns)
{
c *= s;
}
return *this;
}
VTKC_EXEC
Matrix& operator/=(const T& s) noexcept
{
for (auto& c : this->Columns)
{
c /= s;
}
return *this;
}
private:
Vector<T, NumberOfRows> Columns[NumberOfColumns];
};
//---------------------------------------------------------------------------
template <typename T, IdComponent NRows, IdComponent NCols>
VTKC_EXEC
Matrix<T, NRows, NCols> operator+(Matrix<T, NRows, NCols> m1,
const Matrix<T, NRows, NCols>& m2) noexcept
{
return m1 += m2;
}
template <typename T, IdComponent NRows, IdComponent NCols>
VTKC_EXEC
Matrix<T, NRows, NCols> operator-(Matrix<T, NRows, NCols> m1,
const Matrix<T, NRows, NCols>& m2) noexcept
{
return m1 -= m2;
}
template <typename T, IdComponent NRows, IdComponent NCols>
VTKC_EXEC
Matrix<T, NRows, NCols> operator*(Matrix<T, NRows, NCols> m, const T& s) noexcept
{
return m *= s;
}
template <typename T, IdComponent NRows, IdComponent NCols>
VTKC_EXEC
Matrix<T, NRows, NCols> operator/(Matrix<T, NRows, NCols> m, const T& s) noexcept
{
return m /= s;
}
//---------------------------------------------------------------------------
template <typename T,
IdComponent NumRow,
IdComponent NumCol,
IdComponent NumInternal>
VTKC_EXEC
Matrix<T, NumRow, NumCol> matrixMultiply(
const Matrix<T, NumRow, NumInternal>& leftFactor,
const Matrix<T, NumInternal, NumCol>& rightFactor) noexcept
{
Matrix<T, NumRow, NumCol> result;
for (IdComponent rowIndex = 0; rowIndex < NumRow; ++rowIndex)
{
for (IdComponent colIndex = 0; colIndex < NumCol; ++colIndex)
{
T sum = T(leftFactor(rowIndex, 0) * rightFactor(0, colIndex));
for (IdComponent internalIndex = 1; internalIndex < NumInternal; internalIndex++)
{
sum = T(sum + (leftFactor(rowIndex, internalIndex) * rightFactor(internalIndex, colIndex)));
}
result(rowIndex, colIndex) = sum;
}
}
return result;
}
template <typename T, IdComponent NumRow, IdComponent NumCol>
VTKC_EXEC
Vector<T, NumRow> matrixMultiply(
const Matrix<T, NumRow, NumCol>& leftFactor,
const Vector<T, NumCol>& rightFactor) noexcept
{
Vector<T, NumRow> result;
for (IdComponent rowIndex = 0; rowIndex < NumRow; ++rowIndex)
{
result[rowIndex] = dot(leftFactor.getRow(rowIndex), rightFactor);
}
return result;
}
template <typename T, IdComponent NumRow, IdComponent NumCol>
VTKC_EXEC
Vector<T, NumCol> matrixMultiply(
const Vector<T, NumRow>& leftFactor,
const Matrix<T, NumRow, NumCol>& rightFactor) noexcept
{
Vector<T, NumCol> result;
for (IdComponent colIndex = 0; colIndex < NumCol; ++colIndex)
{
result[colIndex] = dot(leftFactor, rightFactor.getColumn(colIndex));
}
return result;
}
//---------------------------------------------------------------------------
// The following code is copied from VTK-m
namespace detail
{
template <typename T, IdComponent Size>
VTKC_EXEC
vtkc::ErrorCode matrixLUPFactorFindPivot(Matrix<T, Size, Size>& m,
Vector<IdComponent, Size>& permutation,
IdComponent topCornerIndex,
T& inversionParity) noexcept
{
constexpr T epsilon = std::is_floating_point<T>::value ?
(std::is_same<T, float>::value ? T(1e-5f) : T(1e-9)) :
T{};
IdComponent maxRowIndex = topCornerIndex;
T maxValue = VTKC_MATH_CALL(abs, (m(maxRowIndex, topCornerIndex)));
for (IdComponent rowIndex = topCornerIndex + 1; rowIndex < Size; rowIndex++)
{
T compareValue = VTKC_MATH_CALL(abs, (m(rowIndex, topCornerIndex)));
if (maxValue < compareValue)
{
maxValue = compareValue;
maxRowIndex = rowIndex;
}
}
if (maxValue < epsilon)
{
return vtkc::ErrorCode::MATRIX_LUP_FACTORIZATION_FAILED;
}
if (maxRowIndex != topCornerIndex)
{
// Swap rows in matrix.
for (IdComponent i = 0; i < Size; ++i)
{
auto temp = m(topCornerIndex, i);
m(topCornerIndex, i) = m(maxRowIndex, i);
m(maxRowIndex, i) = temp;
}
// Record change in permutation matrix.
IdComponent maxOriginalRowIndex = permutation[maxRowIndex];
permutation[maxRowIndex] = permutation[topCornerIndex];
permutation[topCornerIndex] = maxOriginalRowIndex;
// Keep track of inversion parity.
inversionParity = -inversionParity;
}
return vtkc::ErrorCode::SUCCESS;
}
// Used with MatrixLUPFactor
template <typename T, IdComponent Size>
VTKC_EXEC
void matrixLUPFactorFindUpperTriangleElements(
Matrix<T, Size, Size>& m, IdComponent topCornerIndex) noexcept
{
// Compute values for upper triangle on row topCornerIndex
for (IdComponent colIndex = topCornerIndex + 1; colIndex < Size; colIndex++)
{
m(topCornerIndex, colIndex) /= m(topCornerIndex, topCornerIndex);
}
// Update the rest of the matrix for calculations on subsequent rows
for (IdComponent rowIndex = topCornerIndex + 1; rowIndex < Size; rowIndex++)
{
for (IdComponent colIndex = topCornerIndex + 1; colIndex < Size; colIndex++)
{
m(rowIndex, colIndex) -= m(rowIndex, topCornerIndex) * m(topCornerIndex, colIndex);
}
}
}
/// Performs an LUP-factorization on the given matrix using Crout's method. The
/// LU-factorization takes a matrix A and decomposes it into a lower triangular
/// matrix L and upper triangular matrix U such that A = LU. The
/// LUP-factorization also allows permutation of A, which makes the
/// decomposition always possible so long as A is not singular. In addition to
/// matrices L and U, LUP also finds permutation matrix P containing all zeros
/// except one 1 per row and column such that PA = LU.
///
/// The result is done in place such that the lower triangular matrix, L, is
/// stored in the lower-left triangle of A including the diagonal. The upper
/// triangular matrix, U, is stored in the upper-right triangle of L not
/// including the diagonal. The diagonal of U in Crout's method is all 1's (and
/// therefore not explicitly stored).
///
/// The permutation matrix P is represented by the permutation vector. If
/// permutation[i] = j then row j in the original matrix A has been moved to
/// row i in the resulting matrices. The permutation matrix P can be
/// represented by a matrix with p_i,j = 1 if permutation[i] = j and 0
/// otherwise. If using LUP-factorization to compute a determinant, you also
/// need to know the parity (whether there is an odd or even amount) of
/// inversions. An inversion is an instance of a smaller number appearing after
/// a larger number in permutation. Although you can compute the inversion
/// parity after the fact, this function keeps track of it with much less
/// compute resources. The parameter inversionParity is set to 1.0 for even
/// parity and -1.0 for odd parity.
///
/// Not all matrices (specifically singular matrices) have an
/// LUP-factorization. If the LUP-factorization succeeds, valid is set to true.
/// Otherwise, valid is set to false and the result is indeterminant.
///
template <typename T, IdComponent Size>
VTKC_EXEC
vtkc::ErrorCode matrixLUPFactor(Matrix<T, Size, Size>& m,
Vector<IdComponent, Size>& permutation,
T& inversionParity) noexcept
{
// Initialize permutation.
for (IdComponent index = 0; index < Size; index++)
{
permutation[index] = index;
}
inversionParity = T(1);
for (IdComponent rowIndex = 0; rowIndex < Size; rowIndex++)
{
VTKC_RETURN_ON_ERROR(matrixLUPFactorFindPivot(m, permutation, rowIndex, inversionParity))
matrixLUPFactorFindUpperTriangleElements(m, rowIndex);
}
return vtkc::ErrorCode::SUCCESS;
}
/// Use a previous factorization done with MatrixLUPFactor to solve the
/// system Ax = b. Instead of A, this method takes in the LU and P
/// matrices calculated by MatrixLUPFactor from A. The x matrix is returned.
///
template <typename T, IdComponent Size>
VTKC_EXEC
Vector<T, Size> matrixLUPSolve(
const Matrix<T, Size, Size>& LU,
const Vector<IdComponent, Size>& permutation,
const Vector<T, Size>& b) noexcept
{
// The LUP-factorization gives us PA = LU or equivalently A = inv(P)LU.
// Substituting into Ax = b gives us inv(P)LUx = b or LUx = Pb.
// Now consider the intermediate vector y = Ux.
// Substituting in the previous two equations yields Ly = Pb.
// Solving Ly = Pb is easy because L is triangular and P is just a
// permutation.
Vector<T, Size> y;
for (IdComponent rowIndex = 0; rowIndex < Size; rowIndex++)
{
y[rowIndex] = b[permutation[rowIndex]];
// Recall that L is stored in the lower triangle of LU including diagonal.
for (IdComponent colIndex = 0; colIndex < rowIndex; colIndex++)
{
y[rowIndex] -= LU(rowIndex, colIndex) * y[colIndex];
}
y[rowIndex] /= LU(rowIndex, rowIndex);
}
// Now that we have y, we can easily solve Ux = y for x.
Vector<T, Size> x;
for (IdComponent rowIndex = Size - 1; rowIndex >= 0; rowIndex--)
{
// Recall that U is stored in the upper triangle of LU with the diagonal
// implicitly all 1's.
x[rowIndex] = y[rowIndex];
for (IdComponent colIndex = rowIndex + 1; colIndex < Size; colIndex++)
{
x[rowIndex] -= LU(rowIndex, colIndex) * x[colIndex];
}
}
return x;
}
} // namespace detail
template <typename T, IdComponent Size>
VTKC_EXEC
vtkc::ErrorCode solveLinearSystem(const Matrix<T, Size, Size>& A,
const Vector<T, Size>& b,
Vector<T, Size>& x) noexcept
{
// First, we will make an LUP-factorization to help us.
Matrix<T, Size, Size> LU = A;
Vector<IdComponent, Size> permutation;
T inversionParity; // Unused.
VTKC_RETURN_ON_ERROR(detail::matrixLUPFactor(LU, permutation, inversionParity))
// Next, use the decomposition to solve the system.
x = detail::matrixLUPSolve(LU, permutation, b);
return vtkc::ErrorCode::SUCCESS;
}
/// Find and return the inverse of the given matrix. If the matrix is singular,
/// the inverse will not be correct and valid will be set to false.
///
template <typename T, IdComponent Size>
VTKC_EXEC
vtkc::ErrorCode matrixInverse(const Matrix<T, Size, Size>& A, Matrix<T, Size, Size>& invA) noexcept
{
// First, we will make an LUP-factorization to help us.
Matrix<T, Size, Size> LU = A;
Vector<IdComponent, Size> permutation;
T inversionParity; // Unused
VTKC_RETURN_ON_ERROR(detail::matrixLUPFactor(LU, permutation, inversionParity))
// We will use the decomposition to solve AX = I for X where X is
// clearly the inverse of A. Our solve method only works for vectors,
// so we solve for one column of invA at a time.
Vector<T, Size> ICol(T(0));
for (IdComponent colIndex = 0; colIndex < Size; colIndex++)
{
ICol[colIndex] = 1;
Vector<T, Size> invACol = detail::matrixLUPSolve(LU, permutation, ICol);
ICol[colIndex] = 0;
for (IdComponent i = 0; i < Size; ++i)
{
invA(i, colIndex) = invACol[i];
}
}
return vtkc::ErrorCode::SUCCESS;
}
///=========================================================================
template <typename JacobianFunctor, typename FunctionFunctor, typename T, IdComponent Size>
VTKC_EXEC
inline vtkc::ErrorCode newtonsMethod(
const JacobianFunctor& jacobianEvaluator,
const FunctionFunctor& functionEvaluator,
const Vector<T, Size>& rhs,
Vector<T, Size>& result,
T convergeDifference = 1e-3f,
int maxIterations = 10) noexcept
{
auto x = result; // pass initial guess in result
bool converged = false;
for (int i = 0; !converged && i < maxIterations; ++i)
{
Matrix<T, Size, Size> jacobian;
Vector<T, Size> fx(0);
VTKC_RETURN_ON_ERROR(jacobianEvaluator(x, jacobian))
VTKC_RETURN_ON_ERROR(functionEvaluator(x, fx))
Vector<T, Size> deltax;
VTKC_RETURN_ON_ERROR(solveLinearSystem(jacobian, fx - rhs, deltax))
x -= deltax;
converged = true;
for (int c = 0; c < Size; ++c)
{
converged &= (VTKC_MATH_CALL(abs, (deltax[c])) < convergeDifference);
}
}
result = x;
return converged ? ErrorCode::SUCCESS : ErrorCode::SOLUTION_DID_NOT_CONVERGE;
}
///=========================================================================
template <typename T>
VTKC_EXEC inline T lerp(T v0, T v1, T t)
{
static_assert(std::is_floating_point<T>::value,
"lerp requires floating point parameters");
return VTKC_MATH_CALL(fma, (t), (v1), (VTKC_MATH_CALL(fma, (-t), (v0), (v0))));
}
}
} // vtkc::internal
#endif // vtk_c_internal_Math_h

251
vtkm/thirdparty/vtkc/vtkmvtkc/vtkc/vtkc.h vendored Normal file
View File

@ -0,0 +1,251 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.md 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_c_vtkc_h
#define vtk_c_vtkc_h
#include <vtkc/Hexahedron.h>
#include <vtkc/Line.h>
#include <vtkc/Pixel.h>
#include <vtkc/Polygon.h>
#include <vtkc/Pyramid.h>
#include <vtkc/Quad.h>
#include <vtkc/Tetra.h>
#include <vtkc/Triangle.h>
#include <vtkc/Vertex.h>
#include <vtkc/Voxel.h>
#include <vtkc/Wedge.h>
#include <utility>
namespace vtkc
{
/// \brief Perform basic checks to validate cell's state.
/// \param[in] tag The cell tag to validate.
/// \return vtkc::ErrorCode::SUCCESS if valid.
///
VTKC_EXEC inline vtkc::ErrorCode validate(Cell tag) noexcept
{
ErrorCode status = ErrorCode::SUCCESS;
switch (tag.shape())
{
vtkcGenericCellShapeMacro(status = validate(CellTag{tag}));
default:
status = ErrorCode::INVALID_SHAPE_ID;
break;
}
return status;
}
/// \brief Return center of a cell in parametric coordinates.
/// \remark Note that the parametric center is not always located at (0.5,0.5,0.5).
/// \param[in] tag The cell tag.
/// \param[out] pcoords The center of the cell in parametric coordinates.
/// \return A vtkc::ErrorCode value indicating the status of the operation.
///
template<typename CoordType>
VTKC_EXEC inline vtkc::ErrorCode parametricCenter(Cell tag, CoordType&& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
ErrorCode status = ErrorCode::SUCCESS;
switch (tag.shape())
{
vtkcGenericCellShapeMacro(status = parametricCenter(CellTag{tag}, std::forward<CoordType>(pcoords)));
default:
status = ErrorCode::INVALID_SHAPE_ID;
break;
}
return status;
}
/// \brief Return the parametric coordinates of a cell's point.
/// \param[in] tag The cell tag.
/// \param[in] pointId The point number.
/// \param[out] pcoords The parametric coordinates of a cell's point.
/// \return A vtkc::ErrorCode value indicating the status of the operation.
///
template<typename CoordType>
VTKC_EXEC inline vtkc::ErrorCode parametricPoint(
Cell tag, IdComponent pointId, CoordType&& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
ErrorCode status = ErrorCode::SUCCESS;
switch (tag.shape())
{
vtkcGenericCellShapeMacro(status = parametricPoint(CellTag{tag}, pointId, std::forward<CoordType>(pcoords)));
default:
status = ErrorCode::INVALID_SHAPE_ID;
break;
}
return status;
}
/// \brief Return the parametric distance of a parametric coordinate to a cell.
/// \param[in] tag The cell tag.
/// \param[in] pcoords The parametric coordinates of the point.
/// \return The parametric distance of the point to the cell.
/// If point is inside the cell, 0 is returned.
/// \pre tag should be a valid cell, otherwise the result is undefined.
///
template<typename CoordType>
VTKC_EXEC inline ComponentType<CoordType> parametricDistance(Cell tag, const CoordType& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
ComponentType<CoordType> dist{0};
switch (tag.shape())
{
vtkcGenericCellShapeMacro(dist = parametricDistance(CellTag{tag}, pcoords));
default:
break;
}
return dist;
}
/// \brief Check if the given parametric point lies inside a cell.
/// \param[in] tag The cell tag.
/// \param[in] pcoords The parametric coordinates of the point.
/// \return true if inside, false otherwise.
/// \pre tag should be a valid cell, otherwise the result is undefined.
///
template<typename CoordType>
VTKC_EXEC inline bool cellInside(Cell tag, const CoordType& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
bool inside = false;
switch (tag.shape())
{
vtkcGenericCellShapeMacro(inside = cellInside(CellTag{tag}, pcoords));
default:
break;
}
return inside;
}
/// \brief Interpolate \c values at the paramteric coordinates \c pcoords
/// \param[in] tag The cell tag.
/// \param[in] values A \c FieldAccessor for values to interpolate
/// \param[in] pcoords The parametric coordinates.
/// \param[out] result The interpolation result
/// \return A vtkc::ErrorCode value indicating the status of the operation.
///
template <typename Values, typename CoordType, typename Result>
VTKC_EXEC inline vtkc::ErrorCode interpolate(
Cell tag,
const Values& values,
const CoordType& pcoords,
Result&& result) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
ErrorCode status = ErrorCode::SUCCESS;
switch (tag.shape())
{
vtkcGenericCellShapeMacro(status = interpolate(CellTag{tag}, values, pcoords, std::forward<Result>(result)));
default:
status = ErrorCode::INVALID_SHAPE_ID;
}
return status;
}
/// \brief Compute derivative of \c values at the paramteric coordinates \c pcoords
/// \param[in] tag The cell tag.
/// \param[in] points A \c FieldAccessor for points of the cell
/// \param[in] values A \c FieldAccessor for the values to compute derivative of.
/// \param[in] pcoords The parametric coordinates.
/// \param[out] dx The derivative along X
/// \param[out] dy The derivative along Y
/// \param[out] dz The derivative along Z
/// \return A vtkc::ErrorCode value indicating the status of the operation.
///
template <typename Points, typename Values, typename CoordType, typename Result>
VTKC_EXEC inline vtkc::ErrorCode derivative(
Cell tag,
const Points& points,
const Values& values,
const CoordType& pcoords,
Result&& dx,
Result&& dy,
Result&& dz) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
ErrorCode status = ErrorCode::SUCCESS;
switch (tag.shape())
{
vtkcGenericCellShapeMacro(status = derivative(CellTag{tag},
points,
values,
pcoords,
std::forward<Result>(dx),
std::forward<Result>(dy),
std::forward<Result>(dz)));
default:
status = ErrorCode::INVALID_SHAPE_ID;
}
return status;
}
/// \brief Compute world coordinates from parametric coordinates
/// \param[in] tag The cell tag.
/// \param[in] points A \c FieldAccessor for points of the cell
/// \param[in] pcoords The parametric coordinates.
/// \param[out] wcoords The world coordinates.
///
template <typename Points, typename PCoordType, typename WCoordType>
VTKC_EXEC inline vtkc::ErrorCode parametricToWorld(
Cell tag,
const Points& points,
const PCoordType& pcoords,
WCoordType&& wcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(PCoordType);
ErrorCode status = ErrorCode::SUCCESS;
switch (tag.shape())
{
vtkcGenericCellShapeMacro(status = parametricToWorld(CellTag{tag}, points, pcoords, std::forward<WCoordType>(wcoords)));
default:
status = ErrorCode::INVALID_SHAPE_ID;
}
return status;
}
/// \brief Compute parametric coordinates from world coordinates
/// \param[in] tag The cell tag.
/// \param[in] points A \c FieldAccessor for points of the cell
/// \param[in] wcoords The world coordinates.
/// \param[out] pcoords The parametric coordinates.
///
template <typename Points, typename WCoordType, typename PCoordType>
VTKC_EXEC inline vtkc::ErrorCode worldToParametric(
Cell tag,
const Points& points,
const WCoordType& wcoords,
PCoordType&& pcoords) noexcept
{
VTKC_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(PCoordType);
ErrorCode status = ErrorCode::SUCCESS;
switch (tag.shape())
{
vtkcGenericCellShapeMacro(status = worldToParametric(CellTag{tag}, points, wcoords, std::forward<PCoordType>(pcoords)));
default:
status = ErrorCode::INVALID_SHAPE_ID;
}
return status;
}
} //namespace vtkc
#endif //vtk_c_vtkc_h

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();