fdaccc22db
Change the VTKM_CONT_EXPORT to VTKM_CONT. (Likewise for EXEC and EXEC_CONT.) Remove the inline from these macros so that they can be applied to everything, including implementations in a library. Because inline is not declared in these modifies, you have to add the keyword to functions and methods where the implementation is not inlined in the class.
481 lines
18 KiB
C++
481 lines
18 KiB
C++
//============================================================================
|
|
// 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.
|
|
//
|
|
// Copyright 2015 Sandia Corporation.
|
|
// Copyright 2015 UT-Battelle, LLC.
|
|
// Copyright 2015 Los Alamos National Security.
|
|
//
|
|
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
|
|
// 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.
|
|
//============================================================================
|
|
#ifndef vtk_m_exec_Interpolate_h
|
|
#define vtk_m_exec_Interpolate_h
|
|
|
|
#include <vtkm/Assert.h>
|
|
#include <vtkm/CellShape.h>
|
|
#include <vtkm/Math.h>
|
|
#include <vtkm/VecRectilinearPointCoordinates.h>
|
|
#include <vtkm/VectorAnalysis.h>
|
|
#include <vtkm/exec/FunctorBase.h>
|
|
|
|
namespace vtkm {
|
|
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)
|
|
{
|
|
VTKM_ASSERT(pointWCoords.GetNumberOfComponents() == 3);
|
|
|
|
// 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 coordiante. (The v coordinate follows the other edge
|
|
// accordingly.) Thus, if we can find the interesection 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 interesection 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 coordiante is simply d. The v coordinate follows
|
|
// similarly.
|
|
//
|
|
|
|
typedef typename WorldCoordVector::ComponentType Vector3;
|
|
typedef typename Vector3::ComponentType T;
|
|
|
|
Vector3 pcoords(T(0));
|
|
Vector3 triangleNormal =
|
|
vtkm::TriangleNormal(pointWCoords[0], pointWCoords[1], pointWCoords[2]);
|
|
|
|
for (vtkm::IdComponent dimension = 0; dimension < 2; dimension++)
|
|
{
|
|
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;
|
|
}
|
|
|
|
return pcoords;
|
|
}
|
|
|
|
}
|
|
|
|
/// \brief Interpolate a point field in a cell.
|
|
///
|
|
/// Given the point field values for each node and the parametric coordinates
|
|
/// of a point within the cell, interpolates the field to that point.
|
|
///
|
|
template<typename FieldVecType,
|
|
typename ParametricCoordType>
|
|
VTKM_EXEC
|
|
typename FieldVecType::ComponentType
|
|
CellInterpolate(const FieldVecType &pointFieldValues,
|
|
const vtkm::Vec<ParametricCoordType,3> ¶metricCoords,
|
|
vtkm::CellShapeTagGeneric shape,
|
|
const vtkm::exec::FunctorBase &worklet)
|
|
{
|
|
typename FieldVecType::ComponentType result;
|
|
switch (shape.Id)
|
|
{
|
|
vtkmGenericCellShapeMacro(
|
|
result = CellInterpolate(pointFieldValues,
|
|
parametricCoords,
|
|
CellShapeTag(),
|
|
worklet));
|
|
default:
|
|
worklet.RaiseError("Unknown cell shape sent to interpolate.");
|
|
return typename FieldVecType::ComponentType();
|
|
}
|
|
return result;
|
|
}
|
|
|
|
//-----------------------------------------------------------------------------
|
|
template<typename FieldVecType,
|
|
typename ParametricCoordType>
|
|
VTKM_EXEC
|
|
typename FieldVecType::ComponentType
|
|
CellInterpolate(const FieldVecType &,
|
|
const vtkm::Vec<ParametricCoordType,3> &,
|
|
vtkm::CellShapeTagEmpty,
|
|
const vtkm::exec::FunctorBase &worklet)
|
|
{
|
|
worklet.RaiseError("Attempted to interpolate an empty cell.");
|
|
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> ¶metricCoords,
|
|
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::Vec<vtkm::FloatDefault,3>
|
|
CellInterpolate(const vtkm::VecRectilinearPointCoordinates<1> &field,
|
|
const vtkm::Vec<ParametricCoordType,3> &pcoords,
|
|
vtkm::CellShapeTagLine,
|
|
const vtkm::exec::FunctorBase &)
|
|
{
|
|
typedef vtkm::Vec<vtkm::FloatDefault,3> T;
|
|
|
|
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);
|
|
typedef typename FieldVecType::ComponentType T;
|
|
return static_cast<T>( (field[0] * (1 - pcoords[0] - pcoords[1]))
|
|
+ (field[1] * pcoords[0])
|
|
+ (field[2] * pcoords[1]));
|
|
}
|
|
|
|
//-----------------------------------------------------------------------------
|
|
template<typename FieldVecType,
|
|
typename ParametricCoordType>
|
|
VTKM_EXEC
|
|
typename FieldVecType::ComponentType
|
|
CellInterpolate(const FieldVecType &field,
|
|
const vtkm::Vec<ParametricCoordType,3> &pcoords,
|
|
vtkm::CellShapeTagPolygon,
|
|
const vtkm::exec::FunctorBase &worklet)
|
|
{
|
|
const vtkm::IdComponent numPoints = field.GetNumberOfComponents();
|
|
VTKM_ASSERT(numPoints > 0);
|
|
switch (numPoints)
|
|
{
|
|
case 1:
|
|
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);
|
|
}
|
|
|
|
// 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.
|
|
|
|
typedef typename FieldVecType::ComponentType FieldType;
|
|
|
|
// 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 += static_cast<ParametricCoordType>(2*vtkm::Pi());
|
|
}
|
|
const ParametricCoordType deltaAngle =
|
|
static_cast<ParametricCoordType>(2*vtkm::Pi()/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);
|
|
|
|
typedef typename FieldVecType::ComponentType T;
|
|
|
|
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::Vec<vtkm::FloatDefault,3>
|
|
CellInterpolate(const vtkm::VecRectilinearPointCoordinates<2> &field,
|
|
const vtkm::Vec<ParametricCoordType,3> &pcoords,
|
|
vtkm::CellShapeTagQuad,
|
|
const vtkm::exec::FunctorBase &)
|
|
{
|
|
typedef vtkm::Vec<vtkm::FloatDefault,3> T;
|
|
|
|
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]);
|
|
}
|
|
|
|
//-----------------------------------------------------------------------------
|
|
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);
|
|
typedef typename FieldVecType::ComponentType T;
|
|
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);
|
|
|
|
typedef typename FieldVecType::ComponentType T;
|
|
|
|
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::Vec<vtkm::FloatDefault,3>
|
|
CellInterpolate(const vtkm::VecRectilinearPointCoordinates<3> &field,
|
|
const vtkm::Vec<ParametricCoordType,3> &pcoords,
|
|
vtkm::CellShapeTagHexahedron,
|
|
const vtkm::exec::FunctorBase &)
|
|
{
|
|
vtkm::Vec<vtkm::FloatDefault,3> 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);
|
|
|
|
typedef typename FieldVecType::ComponentType T;
|
|
|
|
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);
|
|
|
|
typedef typename FieldVecType::ComponentType T;
|
|
|
|
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]);
|
|
}
|
|
|
|
}
|
|
} // namespace vtkm::exec
|
|
|
|
#endif //vtk_m_exec_Interpolate_h
|