cc497e6a1b
These asserts are consolidated into the unified Assert.h. Also made some minor edits to add asserts where appropriate and a little bit of reconfiguring as found.
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_EXPORT
|
|
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_EXPORT
|
|
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_EXPORT
|
|
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_EXPORT
|
|
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_EXPORT
|
|
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_EXPORT
|
|
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_EXPORT
|
|
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_EXPORT
|
|
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_EXPORT
|
|
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_EXPORT
|
|
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_EXPORT
|
|
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_EXPORT
|
|
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_EXPORT
|
|
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_EXPORT
|
|
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_EXPORT
|
|
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
|