75a060a243
For polygon cell shapes (that are not triangles or quadrilaterals), interpolations are done by finding the center point and creating a triangle fan around that point. Previously, the gradient was computed in the same way as interpolation: identifying the correct triangle and computing the gradient for that triangle. The problem with that approach is that makes the gradient discontinuous at the boundaries of this implicit triangle fan. To make things worse, this discontinuity happens right at each vertex where gradient calculations happen frequently. This means that when you ask for the gradient at the vertex, you might get wildly different answers based on floating point imprecision. Get around this problem by creating a small triangle around the point in question, interpolating values to that triangle, and use that for the gradient. This makes for a smoother gradient transition around these internal boundaries.
1199 lines
49 KiB
C++
1199 lines
49 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.
|
|
//============================================================================
|
|
#ifndef vtk_m_exec_Derivative_h
|
|
#define vtk_m_exec_Derivative_h
|
|
|
|
#include <vtkm/Assert.h>
|
|
#include <vtkm/BaseComponent.h>
|
|
#include <vtkm/CellShape.h>
|
|
#include <vtkm/Math.h>
|
|
#include <vtkm/Matrix.h>
|
|
#include <vtkm/VecAxisAlignedPointCoordinates.h>
|
|
#include <vtkm/VectorAnalysis.h>
|
|
|
|
#include <vtkm/exec/CellInterpolate.h>
|
|
#include <vtkm/exec/FunctorBase.h>
|
|
#include <vtkm/exec/Jacobian.h>
|
|
|
|
namespace vtkm
|
|
{
|
|
namespace exec
|
|
{
|
|
|
|
// The derivative for a 2D polygon in 3D space is underdetermined since there
|
|
// is no information in the direction perpendicular to the polygon. To compute
|
|
// derivatives for general polygons, we build a 2D space for the polygon's
|
|
// plane and solve the derivative there.
|
|
|
|
namespace
|
|
{
|
|
#define VTKM_ACCUM_PARAMETRIC_DERIVATIVE_3D(pointIndex, weight0, weight1, weight2) \
|
|
parametricDerivative[0] += field[pointIndex] * (weight0); \
|
|
parametricDerivative[1] += field[pointIndex] * (weight1); \
|
|
parametricDerivative[2] += field[pointIndex] * (weight2)
|
|
|
|
// Find the derivative of a field in parametric space. That is, find the
|
|
// vector [ds/du, ds/dv, ds/dw].
|
|
template <typename FieldVecType, typename ParametricCoordType>
|
|
VTKM_EXEC vtkm::Vec<typename FieldVecType::ComponentType, 3> ParametricDerivative(
|
|
const FieldVecType& field,
|
|
const vtkm::Vec<ParametricCoordType, 3>& pcoords,
|
|
vtkm::CellShapeTagHexahedron)
|
|
{
|
|
using FieldType = typename FieldVecType::ComponentType;
|
|
using GradientType = vtkm::Vec<FieldType, 3>;
|
|
|
|
GradientType pc(pcoords);
|
|
GradientType rc = GradientType(FieldType(1)) - pc;
|
|
|
|
GradientType parametricDerivative(FieldType(0));
|
|
VTKM_DERIVATIVE_WEIGHTS_HEXAHEDRON(pc, rc, VTKM_ACCUM_PARAMETRIC_DERIVATIVE_3D);
|
|
return parametricDerivative;
|
|
}
|
|
|
|
template <typename FieldVecType, typename ParametricCoordType>
|
|
VTKM_EXEC vtkm::Vec<typename FieldVecType::ComponentType, 3> ParametricDerivative(
|
|
const FieldVecType& field,
|
|
const vtkm::Vec<ParametricCoordType, 3>& pcoords,
|
|
vtkm::CellShapeTagWedge)
|
|
{
|
|
using FieldType = typename FieldVecType::ComponentType;
|
|
using GradientType = vtkm::Vec<FieldType, 3>;
|
|
|
|
GradientType pc(pcoords);
|
|
GradientType rc = GradientType(FieldType(1)) - pc;
|
|
|
|
GradientType parametricDerivative(FieldType(0));
|
|
VTKM_DERIVATIVE_WEIGHTS_WEDGE(pc, rc, VTKM_ACCUM_PARAMETRIC_DERIVATIVE_3D);
|
|
|
|
return parametricDerivative;
|
|
}
|
|
|
|
template <typename FieldVecType, typename ParametricCoordType>
|
|
VTKM_EXEC vtkm::Vec<typename FieldVecType::ComponentType, 3> ParametricDerivative(
|
|
const FieldVecType& field,
|
|
const vtkm::Vec<ParametricCoordType, 3>& pcoords,
|
|
vtkm::CellShapeTagPyramid)
|
|
{
|
|
using FieldType = typename FieldVecType::ComponentType;
|
|
using GradientType = vtkm::Vec<FieldType, 3>;
|
|
|
|
GradientType pc(pcoords);
|
|
GradientType rc = GradientType(FieldType(1)) - pc;
|
|
|
|
GradientType parametricDerivative(FieldType(0));
|
|
VTKM_DERIVATIVE_WEIGHTS_PYRAMID(pc, rc, VTKM_ACCUM_PARAMETRIC_DERIVATIVE_3D);
|
|
|
|
return parametricDerivative;
|
|
}
|
|
|
|
#undef VTKM_ACCUM_PARAMETRIC_DERIVATIVE_3D
|
|
|
|
#define VTKM_ACCUM_PARAMETRIC_DERIVATIVE_2D(pointIndex, weight0, weight1) \
|
|
parametricDerivative[0] += field[pointIndex] * (weight0); \
|
|
parametricDerivative[1] += field[pointIndex] * (weight1)
|
|
|
|
template <typename FieldVecType, typename ParametricCoordType>
|
|
VTKM_EXEC vtkm::Vec<typename FieldVecType::ComponentType, 2> ParametricDerivative(
|
|
const FieldVecType& field,
|
|
const vtkm::Vec<ParametricCoordType, 3>& pcoords,
|
|
vtkm::CellShapeTagQuad)
|
|
{
|
|
using FieldType = typename FieldVecType::ComponentType;
|
|
using GradientType = vtkm::Vec<FieldType, 2>;
|
|
|
|
GradientType pc(static_cast<FieldType>(pcoords[0]), static_cast<FieldType>(pcoords[1]));
|
|
GradientType rc = GradientType(FieldType(1)) - pc;
|
|
|
|
GradientType parametricDerivative(FieldType(0));
|
|
VTKM_DERIVATIVE_WEIGHTS_QUAD(pc, rc, VTKM_ACCUM_PARAMETRIC_DERIVATIVE_2D);
|
|
|
|
return parametricDerivative;
|
|
}
|
|
|
|
#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 FieldVecType,
|
|
typename ParametricCoordType>
|
|
VTKM_EXEC
|
|
vtkm::Vec<typename FieldVecType::ComponentType,2>
|
|
ParametricDerivative(const FieldVecType &field,
|
|
const vtkm::Vec<ParametricCoordType,3> &pcoords,
|
|
vtkm::CellShapeTagPolygon)
|
|
{
|
|
using FieldType = typename FieldVecType::ComponentType;
|
|
using GradientType = vtkm::Vec<FieldType,2>;
|
|
|
|
const vtkm::IdComponent numPoints = field.GetNumberOfComponents();
|
|
FieldType deltaAngle = static_cast<FieldType>(2*vtkm::Pi()/numPoints);
|
|
|
|
GradientType pc(pcoords[0], pcoords[1]);
|
|
|
|
GradientType parametricDerivative(0);
|
|
for (vtkm::IdComponent pointIndex = 0; pointIndex < numPoints; pointIndex++)
|
|
{
|
|
FieldType angle = pointIndex*deltaAngle;
|
|
vtkm::Vec<FieldType,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<FieldType,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.
|
|
FieldType magSqr = vtkm::MagnitudeSquared(pvec);
|
|
FieldType invMag = vtkm::RSqrt(magSqr);
|
|
FieldType scale = invMag*invMag*invMag;
|
|
vtkm::Vec<FieldType,2> weight = scale*pvec;
|
|
|
|
parametricDerivative[0] += field[pointIndex] * weight[0];
|
|
parametricDerivative[1] += field[pointIndex] * weight[1];
|
|
}
|
|
|
|
return parametricDerivative;
|
|
}
|
|
#endif
|
|
|
|
#undef VTKM_ACCUM_PARAMETRIC_DERIVATIVE_2D
|
|
|
|
} // namespace unnamed
|
|
|
|
namespace detail
|
|
{
|
|
|
|
template <typename FieldType,
|
|
int NumPoints,
|
|
typename WorldCoordType,
|
|
typename ParametricCoordType,
|
|
typename CellShapeTag>
|
|
VTKM_EXEC vtkm::Vec<FieldType, 3> CellDerivativeFor3DCell(
|
|
const vtkm::Vec<FieldType, NumPoints>& field,
|
|
const WorldCoordType& wCoords,
|
|
const vtkm::Vec<ParametricCoordType, 3>& pcoords,
|
|
CellShapeTag tag)
|
|
{
|
|
using GradientType = vtkm::Vec<FieldType, 3>;
|
|
|
|
// For reasons that should become apparent in a moment, we actually want
|
|
// the transpose of the Jacobian.
|
|
vtkm::Matrix<FieldType, 3, 3> jacobianTranspose;
|
|
vtkm::exec::JacobianFor3DCell(wCoords, pcoords, jacobianTranspose, tag);
|
|
jacobianTranspose = vtkm::MatrixTranspose(jacobianTranspose);
|
|
|
|
GradientType parametricDerivative = ParametricDerivative(field, pcoords, tag);
|
|
|
|
// If we write out the matrices below, it should become clear that the
|
|
// Jacobian transpose times the field derivative in world space equals
|
|
// the field derivative in parametric space.
|
|
//
|
|
// | | | | | |
|
|
// | dx/du dy/du dz/du | | ds/dx | | ds/du |
|
|
// | | | | | |
|
|
// | dx/dv dy/dv dz/dv | | ds/dy | = | ds/dv |
|
|
// | | | | | |
|
|
// | dx/dw dy/dw dz/dw | | ds/dz | | ds/dw |
|
|
// | | | | | |
|
|
//
|
|
// Now we just need to solve this linear system to find the derivative in
|
|
// world space.
|
|
|
|
bool valid; // Ignored.
|
|
return vtkm::SolveLinearSystem(jacobianTranspose, parametricDerivative, valid);
|
|
}
|
|
|
|
template <typename T,
|
|
int N,
|
|
int NumPoints,
|
|
typename WorldCoordType,
|
|
typename ParametricCoordType,
|
|
typename CellShapeTag>
|
|
VTKM_EXEC vtkm::Vec<vtkm::Vec<T, N>, 3> CellDerivativeFor3DCell(
|
|
const vtkm::Vec<vtkm::Vec<T, N>, NumPoints>& field,
|
|
const WorldCoordType& wCoords,
|
|
const vtkm::Vec<ParametricCoordType, 3>& pcoords,
|
|
CellShapeTag tag)
|
|
{
|
|
//We have been given a vector field so we need to solve for each
|
|
//component of the vector. For explanation of the logic used see the
|
|
//scalar version of CellDerivativeFor3DCell.
|
|
vtkm::Matrix<T, 3, 3> perComponentJacobianTranspose;
|
|
vtkm::exec::JacobianFor3DCell(wCoords, pcoords, perComponentJacobianTranspose, tag);
|
|
perComponentJacobianTranspose = vtkm::MatrixTranspose(perComponentJacobianTranspose);
|
|
|
|
bool valid; // Ignored.
|
|
vtkm::Vec<vtkm::Vec<T, N>, 3> result(vtkm::Vec<T, N>(0.0f));
|
|
vtkm::Vec<T, NumPoints> perPoint;
|
|
for (vtkm::IdComponent i = 0; i < N; ++i)
|
|
{
|
|
for (vtkm::IdComponent c = 0; c < NumPoints; ++c)
|
|
{
|
|
perPoint[c] = field[c][i];
|
|
}
|
|
vtkm::Vec<T, 3> p = ParametricDerivative(perPoint, pcoords, tag);
|
|
vtkm::Vec<T, 3> grad = vtkm::SolveLinearSystem(perComponentJacobianTranspose, p, valid);
|
|
|
|
result[0][i] += grad[0];
|
|
result[1][i] += grad[1];
|
|
result[2][i] += grad[2];
|
|
}
|
|
return result;
|
|
}
|
|
|
|
template <typename FieldType, typename LUType, typename ParametricCoordType, typename CellShapeTag>
|
|
VTKM_EXEC vtkm::Vec<FieldType, 3> CellDerivativeFor2DCellFinish(
|
|
const vtkm::Vec<FieldType, 4>& field,
|
|
const vtkm::Matrix<LUType, 2, 2>& LUFactorization,
|
|
const vtkm::IdComponent2& LUPermutation,
|
|
const vtkm::Vec<ParametricCoordType, 3>& pcoords,
|
|
const vtkm::exec::internal::Space2D<LUType>& space,
|
|
CellShapeTag,
|
|
vtkm::TypeTraitsScalarTag)
|
|
{
|
|
// Finish solving linear equation. See CellDerivativeFor2DCell implementation
|
|
// for more detail.
|
|
vtkm::Vec<FieldType, 2> parametricDerivative =
|
|
ParametricDerivative(field, pcoords, CellShapeTag());
|
|
|
|
vtkm::Vec<FieldType, 2> gradient2D =
|
|
vtkm::detail::MatrixLUPSolve(LUFactorization, LUPermutation, parametricDerivative);
|
|
|
|
return space.ConvertVecFromSpace(gradient2D);
|
|
}
|
|
|
|
template <typename FieldType, typename LUType, typename ParametricCoordType, typename CellShapeTag>
|
|
VTKM_EXEC vtkm::Vec<FieldType, 3> CellDerivativeFor2DCellFinish(
|
|
const vtkm::Vec<FieldType, 4>& field,
|
|
const vtkm::Matrix<LUType, 2, 2>& LUFactorization,
|
|
const vtkm::IdComponent2& LUPermutation,
|
|
const vtkm::Vec<ParametricCoordType, 3>& pcoords,
|
|
const vtkm::exec::internal::Space2D<LUType>& space,
|
|
CellShapeTag,
|
|
vtkm::TypeTraitsVectorTag)
|
|
{
|
|
using FieldTraits = vtkm::VecTraits<FieldType>;
|
|
using FieldComponentType = typename FieldTraits::ComponentType;
|
|
|
|
vtkm::Vec<FieldType, 3> gradient(vtkm::TypeTraits<FieldType>::ZeroInitialization());
|
|
|
|
for (vtkm::IdComponent fieldComponent = 0;
|
|
fieldComponent < FieldTraits::GetNumberOfComponents(field[0]);
|
|
fieldComponent++)
|
|
{
|
|
vtkm::Vec<FieldComponentType, 4> subField(FieldTraits::GetComponent(field[0], fieldComponent),
|
|
FieldTraits::GetComponent(field[1], fieldComponent),
|
|
FieldTraits::GetComponent(field[2], fieldComponent),
|
|
FieldTraits::GetComponent(field[3], fieldComponent));
|
|
|
|
vtkm::Vec<FieldComponentType, 3> subGradient = CellDerivativeFor2DCellFinish(
|
|
subField,
|
|
LUFactorization,
|
|
LUPermutation,
|
|
pcoords,
|
|
space,
|
|
CellShapeTag(),
|
|
typename vtkm::TypeTraits<FieldComponentType>::DimensionalityTag());
|
|
|
|
FieldTraits::SetComponent(gradient[0], fieldComponent, subGradient[0]);
|
|
FieldTraits::SetComponent(gradient[1], fieldComponent, subGradient[1]);
|
|
FieldTraits::SetComponent(gradient[2], fieldComponent, subGradient[2]);
|
|
}
|
|
|
|
return gradient;
|
|
}
|
|
|
|
template <typename FieldType, typename LUType, typename ParametricCoordType, typename CellShapeTag>
|
|
VTKM_EXEC vtkm::Vec<FieldType, 3> CellDerivativeFor2DCellFinish(
|
|
const vtkm::Vec<FieldType, 4>& field,
|
|
const vtkm::Matrix<LUType, 2, 2>& LUFactorization,
|
|
const vtkm::IdComponent2& LUPermutation,
|
|
const vtkm::Vec<ParametricCoordType, 3>& pcoords,
|
|
const vtkm::exec::internal::Space2D<LUType>& space,
|
|
CellShapeTag,
|
|
vtkm::TypeTraitsMatrixTag)
|
|
{
|
|
return CellDerivativeFor2DCellFinish(field,
|
|
LUFactorization,
|
|
LUPermutation,
|
|
pcoords,
|
|
space,
|
|
CellShapeTag(),
|
|
vtkm::TypeTraitsVectorTag());
|
|
}
|
|
|
|
template <typename FieldVecType,
|
|
typename WorldCoordType,
|
|
typename ParametricCoordType,
|
|
typename CellShapeTag>
|
|
VTKM_EXEC vtkm::Vec<typename FieldVecType::ComponentType, 3> CellDerivativeFor2DCell(
|
|
const FieldVecType& field,
|
|
const WorldCoordType& wCoords,
|
|
const vtkm::Vec<ParametricCoordType, 3>& pcoords,
|
|
CellShapeTag)
|
|
{
|
|
using FieldType = typename FieldVecType::ComponentType;
|
|
using BaseFieldType = typename BaseComponent<FieldType>::Type;
|
|
|
|
// We have an underdetermined system in 3D, so create a 2D space in the
|
|
// plane that the polygon sits.
|
|
vtkm::exec::internal::Space2D<BaseFieldType> space(
|
|
wCoords[0], wCoords[1], wCoords[wCoords.GetNumberOfComponents() - 1]);
|
|
|
|
// For reasons that should become apparent in a moment, we actually want
|
|
// the transpose of the Jacobian.
|
|
vtkm::Matrix<BaseFieldType, 2, 2> jacobianTranspose;
|
|
vtkm::exec::JacobianFor2DCell(wCoords, pcoords, space, jacobianTranspose, CellShapeTag());
|
|
jacobianTranspose = vtkm::MatrixTranspose(jacobianTranspose);
|
|
|
|
// Find the derivative of the field in parametric coordinate space. That is,
|
|
// find the vector [ds/du, ds/dv].
|
|
// Commented because this is actually done in CellDerivativeFor2DCellFinish
|
|
// to handle vector fields.
|
|
// vtkm::Vec<BaseFieldType,2> parametricDerivative =
|
|
// ParametricDerivative(field,pcoords,CellShapeTag());
|
|
|
|
// If we write out the matrices below, it should become clear that the
|
|
// Jacobian transpose times the field derivative in world space equals
|
|
// the field derivative in parametric space.
|
|
//
|
|
// | | | | | |
|
|
// | db0/du db1/du | | ds/db0 | | ds/du |
|
|
// | | | | = | |
|
|
// | db0/dv db1/dv | | ds/db1 | | ds/dv |
|
|
// | | | | | |
|
|
//
|
|
// Now we just need to solve this linear system to find the derivative in
|
|
// world space.
|
|
|
|
bool valid; // Ignored.
|
|
// If you look at the implementation of vtkm::SolveLinearSystem, you will see
|
|
// that it is done in two parts. First, it does an LU factorization and then
|
|
// uses that result to complete the solve. The factorization part talkes the
|
|
// longest amount of time, and if we are performing the gradient on a vector
|
|
// field, the factorization can be reused for each component of the vector
|
|
// field. Thus, we are going to call the internals of SolveLinearSystem
|
|
// ourselves to do the factorization and then apply it to all components.
|
|
vtkm::IdComponent2 permutation;
|
|
BaseFieldType inversionParity; // Unused
|
|
vtkm::detail::MatrixLUPFactor(jacobianTranspose, permutation, inversionParity, valid);
|
|
// MatrixLUPFactor does in place factorization. jacobianTranspose is now the
|
|
// LU factorization.
|
|
return CellDerivativeFor2DCellFinish(field,
|
|
jacobianTranspose,
|
|
permutation,
|
|
pcoords,
|
|
space,
|
|
CellShapeTag(),
|
|
typename vtkm::TypeTraits<FieldType>::DimensionalityTag());
|
|
}
|
|
|
|
} // namespace detail
|
|
|
|
//-----------------------------------------------------------------------------
|
|
/// \brief Take the derivative (get the gradient) of a point field in a cell.
|
|
///
|
|
/// Given the point field values for each node and the parametric coordinates
|
|
/// of a point within the cell, finds the derivative with respect to each
|
|
/// coordinate (i.e. the gradient) at that point. The derivative is not always
|
|
/// constant in some "linear" cells.
|
|
///
|
|
template <typename FieldVecType, typename WorldCoordType, typename ParametricCoordType>
|
|
VTKM_EXEC vtkm::Vec<typename FieldVecType::ComponentType, 3> CellDerivative(
|
|
const FieldVecType& pointFieldValues,
|
|
const WorldCoordType& worldCoordinateValues,
|
|
const vtkm::Vec<ParametricCoordType, 3>& parametricCoords,
|
|
vtkm::CellShapeTagGeneric shape,
|
|
const vtkm::exec::FunctorBase& worklet)
|
|
{
|
|
vtkm::Vec<typename FieldVecType::ComponentType, 3> result;
|
|
switch (shape.Id)
|
|
{
|
|
vtkmGenericCellShapeMacro(
|
|
result = CellDerivative(
|
|
pointFieldValues, worldCoordinateValues, parametricCoords, CellShapeTag(), worklet));
|
|
default:
|
|
worklet.RaiseError("Unknown cell shape sent to derivative.");
|
|
return vtkm::Vec<typename FieldVecType::ComponentType, 3>();
|
|
}
|
|
return result;
|
|
}
|
|
|
|
//-----------------------------------------------------------------------------
|
|
template <typename FieldVecType, typename WorldCoordType, typename ParametricCoordType>
|
|
VTKM_EXEC vtkm::Vec<typename FieldVecType::ComponentType, 3> CellDerivative(
|
|
const FieldVecType&,
|
|
const WorldCoordType&,
|
|
const vtkm::Vec<ParametricCoordType, 3>&,
|
|
vtkm::CellShapeTagEmpty,
|
|
const vtkm::exec::FunctorBase& worklet)
|
|
{
|
|
worklet.RaiseError("Attempted to take derivative in empty cell.");
|
|
return vtkm::Vec<typename FieldVecType::ComponentType, 3>();
|
|
}
|
|
|
|
//-----------------------------------------------------------------------------
|
|
template <typename FieldVecType, typename WorldCoordType, typename ParametricCoordType>
|
|
VTKM_EXEC vtkm::Vec<typename FieldVecType::ComponentType, 3> CellDerivative(
|
|
const FieldVecType& field,
|
|
const WorldCoordType& wCoords,
|
|
const vtkm::Vec<ParametricCoordType, 3>&,
|
|
vtkm::CellShapeTagVertex,
|
|
const vtkm::exec::FunctorBase& vtkmNotUsed(worklet))
|
|
{
|
|
(void)field;
|
|
(void)wCoords;
|
|
VTKM_ASSERT(field.GetNumberOfComponents() == 1);
|
|
VTKM_ASSERT(wCoords.GetNumberOfComponents() == 1);
|
|
|
|
using GradientType = vtkm::Vec<typename FieldVecType::ComponentType, 3>;
|
|
return vtkm::TypeTraits<GradientType>::ZeroInitialization();
|
|
}
|
|
|
|
//-----------------------------------------------------------------------------
|
|
namespace detail
|
|
{
|
|
|
|
template <typename FieldType, typename WorldCoordType>
|
|
VTKM_EXEC vtkm::Vec<FieldType, 3> CellDerivativeLineImpl(
|
|
const FieldType& deltaField,
|
|
const WorldCoordType& vec, // direction of line
|
|
const typename WorldCoordType::ComponentType& vecMagSqr,
|
|
vtkm::TypeTraitsScalarTag)
|
|
{
|
|
using GradientType = vtkm::Vec<FieldType, 3>;
|
|
|
|
// The derivative of a line is in the direction of the line. Its length is
|
|
// equal to the difference of the scalar field divided by the length of the
|
|
// line segment. Thus, the derivative is characterized by
|
|
// (deltaField*vec)/mag(vec)^2.
|
|
|
|
return (deltaField / static_cast<FieldType>(vecMagSqr)) * GradientType(vec);
|
|
}
|
|
|
|
template <typename FieldType, typename WorldCoordType, typename VectorTypeTraitsTag>
|
|
VTKM_EXEC vtkm::Vec<FieldType, 3> CellDerivativeLineImpl(
|
|
const FieldType& deltaField,
|
|
const WorldCoordType& vec, // direction of line
|
|
const typename WorldCoordType::ComponentType& vecMagSqr,
|
|
VectorTypeTraitsTag)
|
|
{
|
|
using FieldTraits = vtkm::VecTraits<FieldType>;
|
|
using FieldComponentType = typename FieldTraits::ComponentType;
|
|
using GradientType = vtkm::Vec<FieldType, 3>;
|
|
|
|
GradientType gradient(vtkm::TypeTraits<FieldType>::ZeroInitialization());
|
|
for (vtkm::IdComponent fieldComponent = 0;
|
|
fieldComponent < FieldTraits::GetNumberOfComponents(deltaField);
|
|
fieldComponent++)
|
|
{
|
|
using SubGradientType = vtkm::Vec<FieldComponentType, 3>;
|
|
SubGradientType subGradient =
|
|
CellDerivativeLineImpl(FieldTraits::GetComponent(deltaField, fieldComponent),
|
|
vec,
|
|
vecMagSqr,
|
|
typename vtkm::TypeTraits<FieldComponentType>::DimensionalityTag());
|
|
FieldTraits::SetComponent(gradient[0], fieldComponent, subGradient[0]);
|
|
FieldTraits::SetComponent(gradient[1], fieldComponent, subGradient[1]);
|
|
FieldTraits::SetComponent(gradient[2], fieldComponent, subGradient[2]);
|
|
}
|
|
|
|
return gradient;
|
|
}
|
|
|
|
} // namespace detail
|
|
|
|
template <typename FieldVecType, typename WorldCoordType, typename ParametricCoordType>
|
|
VTKM_EXEC vtkm::Vec<typename FieldVecType::ComponentType, 3> CellDerivative(
|
|
const FieldVecType& field,
|
|
const WorldCoordType& wCoords,
|
|
const vtkm::Vec<ParametricCoordType, 3>& vtkmNotUsed(pcoords),
|
|
vtkm::CellShapeTagLine,
|
|
const vtkm::exec::FunctorBase& vtkmNotUsed(worklet))
|
|
{
|
|
VTKM_ASSERT(field.GetNumberOfComponents() == 2);
|
|
VTKM_ASSERT(wCoords.GetNumberOfComponents() == 2);
|
|
|
|
using FieldType = typename FieldVecType::ComponentType;
|
|
using BaseComponentType = typename BaseComponent<FieldType>::Type;
|
|
|
|
FieldType deltaField(field[1] - field[0]);
|
|
vtkm::Vec<BaseComponentType, 3> vec(wCoords[1] - wCoords[0]);
|
|
|
|
return detail::CellDerivativeLineImpl(deltaField,
|
|
vec,
|
|
vtkm::MagnitudeSquared(vec),
|
|
typename vtkm::TypeTraits<FieldType>::DimensionalityTag());
|
|
}
|
|
|
|
template <typename FieldVecType, typename ParametricCoordType>
|
|
VTKM_EXEC vtkm::Vec<typename FieldVecType::ComponentType, 3> CellDerivative(
|
|
const FieldVecType& field,
|
|
const vtkm::VecAxisAlignedPointCoordinates<1>& wCoords,
|
|
const vtkm::Vec<ParametricCoordType, 3>& vtkmNotUsed(pcoords),
|
|
vtkm::CellShapeTagLine,
|
|
const vtkm::exec::FunctorBase& vtkmNotUsed(worklet))
|
|
{
|
|
VTKM_ASSERT(field.GetNumberOfComponents() == 2);
|
|
|
|
using T = typename FieldVecType::ComponentType;
|
|
|
|
return vtkm::Vec<T, 3>(
|
|
static_cast<T>((field[1] - field[0]) / wCoords.GetSpacing()[0]), T(0), T(0));
|
|
}
|
|
|
|
template <typename FieldVecType, typename WorldCoordType, typename ParametricCoordType>
|
|
VTKM_EXEC vtkm::Vec<typename FieldVecType::ComponentType, 3> CellDerivative(
|
|
const FieldVecType& field,
|
|
const WorldCoordType& wCoords,
|
|
const vtkm::Vec<ParametricCoordType, 3>& pcoords,
|
|
vtkm::CellShapeTagPolyLine,
|
|
const vtkm::exec::FunctorBase& worklet)
|
|
{
|
|
vtkm::IdComponent numPoints = field.GetNumberOfComponents();
|
|
VTKM_ASSERT(numPoints >= 1);
|
|
VTKM_ASSERT(numPoints == wCoords.GetNumberOfComponents());
|
|
|
|
switch (numPoints)
|
|
{
|
|
case 1:
|
|
return CellDerivative(field, wCoords, pcoords, vtkm::CellShapeTagVertex(), worklet);
|
|
case 2:
|
|
return CellDerivative(field, wCoords, pcoords, vtkm::CellShapeTagLine(), worklet);
|
|
}
|
|
|
|
using FieldType = typename FieldVecType::ComponentType;
|
|
using BaseComponentType = typename BaseComponent<FieldType>::Type;
|
|
|
|
ParametricCoordType dt;
|
|
dt = static_cast<ParametricCoordType>(1) / static_cast<ParametricCoordType>(numPoints - 1);
|
|
vtkm::IdComponent idx = static_cast<vtkm::IdComponent>(vtkm::Ceil(pcoords[0] / dt));
|
|
if (idx == 0)
|
|
idx = 1;
|
|
if (idx > numPoints - 1)
|
|
idx = numPoints - 1;
|
|
|
|
FieldType deltaField(field[idx] - field[idx - 1]);
|
|
vtkm::Vec<BaseComponentType, 3> vec(wCoords[idx] - wCoords[idx - 1]);
|
|
|
|
return detail::CellDerivativeLineImpl(deltaField,
|
|
vec,
|
|
vtkm::MagnitudeSquared(vec),
|
|
typename vtkm::TypeTraits<FieldType>::DimensionalityTag());
|
|
}
|
|
|
|
template <typename FieldVecType, typename ParametricCoordType>
|
|
VTKM_EXEC vtkm::Vec<typename FieldVecType::ComponentType, 3> CellDerivative(
|
|
const FieldVecType& field,
|
|
const vtkm::VecAxisAlignedPointCoordinates<1>& wCoords,
|
|
const vtkm::Vec<ParametricCoordType, 3>& pcoords,
|
|
vtkm::CellShapeTagPolyLine,
|
|
const vtkm::exec::FunctorBase& worklet)
|
|
{
|
|
vtkm::IdComponent numPoints = field.GetNumberOfComponents();
|
|
VTKM_ASSERT(numPoints >= 1);
|
|
|
|
switch (numPoints)
|
|
{
|
|
case 1:
|
|
return CellDerivative(field, wCoords, pcoords, vtkm::CellShapeTagVertex(), worklet);
|
|
case 2:
|
|
return CellDerivative(field, wCoords, pcoords, vtkm::CellShapeTagLine(), worklet);
|
|
}
|
|
|
|
ParametricCoordType dt;
|
|
dt = static_cast<ParametricCoordType>(1) / static_cast<ParametricCoordType>(numPoints - 1);
|
|
vtkm::IdComponent idx = static_cast<vtkm::IdComponent>(vtkm::Ceil(pcoords[0] / dt));
|
|
if (idx == 0)
|
|
idx = 1;
|
|
if (idx > numPoints - 1)
|
|
idx = numPoints - 1;
|
|
|
|
using T = typename FieldVecType::ComponentType;
|
|
|
|
return vtkm::Vec<T, 3>(
|
|
static_cast<T>((field[idx] - field[idx - 1]) / wCoords.GetSpacing()[0]), T(0), T(0));
|
|
}
|
|
|
|
//-----------------------------------------------------------------------------
|
|
namespace detail
|
|
{
|
|
|
|
template <typename ValueType, typename LUType>
|
|
VTKM_EXEC vtkm::Vec<ValueType, 3> TriangleDerivativeFinish(
|
|
const vtkm::Vec<ValueType, 3>& field,
|
|
const vtkm::Matrix<LUType, 3, 3>& LUFactorization,
|
|
const vtkm::IdComponent3& LUPermutation,
|
|
vtkm::TypeTraitsScalarTag)
|
|
{
|
|
// Finish solving linear equation. See TriangleDerivative implementation
|
|
// for more detail.
|
|
vtkm::Vec<LUType, 3> b(field[1] - field[0], field[2] - field[0], 0);
|
|
|
|
return vtkm::detail::MatrixLUPSolve(LUFactorization, LUPermutation, b);
|
|
}
|
|
|
|
template <typename ValueType, typename LUType>
|
|
VTKM_EXEC vtkm::Vec<ValueType, 3> TriangleDerivativeFinish(
|
|
const vtkm::Vec<ValueType, 3>& field,
|
|
const vtkm::Matrix<LUType, 3, 3>& LUFactorization,
|
|
const vtkm::IdComponent3& LUPermutation,
|
|
vtkm::TypeTraitsVectorTag)
|
|
{
|
|
using FieldTraits = vtkm::VecTraits<ValueType>;
|
|
using FieldComponentType = typename FieldTraits::ComponentType;
|
|
|
|
vtkm::Vec<ValueType, 3> gradient(vtkm::TypeTraits<ValueType>::ZeroInitialization());
|
|
|
|
for (vtkm::IdComponent fieldComponent = 0;
|
|
fieldComponent < FieldTraits::GetNumberOfComponents(field[0]);
|
|
fieldComponent++)
|
|
{
|
|
vtkm::Vec<FieldComponentType, 3> subField(FieldTraits::GetComponent(field[0], fieldComponent),
|
|
FieldTraits::GetComponent(field[1], fieldComponent),
|
|
FieldTraits::GetComponent(field[2], fieldComponent));
|
|
|
|
vtkm::Vec<FieldComponentType, 3> subGradient =
|
|
TriangleDerivativeFinish(subField,
|
|
LUFactorization,
|
|
LUPermutation,
|
|
typename vtkm::TypeTraits<FieldComponentType>::DimensionalityTag());
|
|
|
|
FieldTraits::SetComponent(gradient[0], fieldComponent, subGradient[0]);
|
|
FieldTraits::SetComponent(gradient[1], fieldComponent, subGradient[1]);
|
|
FieldTraits::SetComponent(gradient[2], fieldComponent, subGradient[2]);
|
|
}
|
|
|
|
return gradient;
|
|
}
|
|
|
|
template <typename ValueType, typename LUType>
|
|
VTKM_EXEC vtkm::Vec<ValueType, 3> TriangleDerivativeFinish(
|
|
const vtkm::Vec<ValueType, 3>& field,
|
|
const vtkm::Matrix<LUType, 3, 3>& LUFactorization,
|
|
const vtkm::IdComponent3& LUPermutation,
|
|
vtkm::TypeTraitsMatrixTag)
|
|
{
|
|
return TriangleDerivativeFinish(
|
|
field, LUFactorization, LUPermutation, vtkm::TypeTraitsVectorTag());
|
|
}
|
|
|
|
template <typename ValueType, typename WCoordType>
|
|
VTKM_EXEC vtkm::Vec<ValueType, 3> TriangleDerivative(const vtkm::Vec<ValueType, 3>& field,
|
|
const vtkm::Vec<WCoordType, 3>& wCoords)
|
|
{
|
|
using BaseComponentType = typename BaseComponent<ValueType>::Type;
|
|
|
|
// The scalar values of the three points in a triangle completely specify a
|
|
// linear field (with constant gradient) assuming the field is constant in
|
|
// the normal direction to the triangle. The field, defined by the 3-vector
|
|
// gradient g and scalar value s_origin, can be found with this set of 4
|
|
// equations and 4 unknowns.
|
|
//
|
|
// Dot(p0, g) + s_origin = s0
|
|
// Dot(p1, g) + s_origin = s1
|
|
// Dot(p2, g) + s_origin = s2
|
|
// Dot(n, g) = 0
|
|
//
|
|
// Where the p's are point coordinates and n is the normal vector. But we
|
|
// don't really care about s_origin. We just want to find the gradient g.
|
|
// With some simple elimination we, we can get rid of s_origin and be left
|
|
// with 3 equations and 3 unknowns.
|
|
//
|
|
// Dot(p1-p0, g) = s1 - s0
|
|
// Dot(p2-p0, g) = s2 - s0
|
|
// Dot(n, g) = 0
|
|
//
|
|
// We'll solve this by putting this in matrix form Ax = b where the rows of A
|
|
// are the differences in points and normal, b has the scalar differences,
|
|
// and x is really the gradient g.
|
|
|
|
vtkm::Vec<BaseComponentType, 3> v0 = wCoords[1] - wCoords[0];
|
|
vtkm::Vec<BaseComponentType, 3> v1 = wCoords[2] - wCoords[0];
|
|
vtkm::Vec<BaseComponentType, 3> n = vtkm::Cross(v0, v1);
|
|
|
|
vtkm::Matrix<BaseComponentType, 3, 3> A;
|
|
vtkm::MatrixSetRow(A, 0, v0);
|
|
vtkm::MatrixSetRow(A, 1, v1);
|
|
vtkm::MatrixSetRow(A, 2, n);
|
|
|
|
// If the triangle is degenerate, then valid will be false. For now we are
|
|
// ignoring it. We could detect it if we determine we need to although I have
|
|
// seen singular matrices missed due to floating point error.
|
|
//
|
|
bool valid;
|
|
|
|
// If you look at the implementation of vtkm::SolveLinearSystem, you will see
|
|
// that it is done in two parts. First, it does an LU factorization and then
|
|
// uses that result to complete the solve. The factorization part talkes the
|
|
// longest amount of time, and if we are performing the gradient on a vector
|
|
// field, the factorization can be reused for each component of the vector
|
|
// field. Thus, we are going to call the internals of SolveLinearSystem
|
|
// ourselves to do the factorization and then apply it to all components.
|
|
vtkm::IdComponent3 permutation;
|
|
BaseComponentType inversionParity; // Unused
|
|
vtkm::detail::MatrixLUPFactor(A, permutation, inversionParity, valid);
|
|
// MatrixLUPFactor does in place factorization. A is now the LU factorization.
|
|
return TriangleDerivativeFinish(
|
|
field, A, permutation, typename vtkm::TypeTraits<ValueType>::DimensionalityTag());
|
|
}
|
|
|
|
} // namespace detail
|
|
|
|
//-----------------------------------------------------------------------------
|
|
template <typename FieldVecType, typename WorldCoordType, typename ParametricCoordType>
|
|
VTKM_EXEC vtkm::Vec<typename FieldVecType::ComponentType, 3> CellDerivative(
|
|
const FieldVecType& inputField,
|
|
const WorldCoordType& wCoords,
|
|
const vtkm::Vec<ParametricCoordType, 3>& vtkmNotUsed(pcoords),
|
|
vtkm::CellShapeTagTriangle,
|
|
const vtkm::exec::FunctorBase& vtkmNotUsed(worklet))
|
|
{
|
|
VTKM_ASSERT(inputField.GetNumberOfComponents() == 3);
|
|
VTKM_ASSERT(wCoords.GetNumberOfComponents() == 3);
|
|
|
|
using ValueType = typename FieldVecType::ComponentType;
|
|
using WCoordType = typename WorldCoordType::ComponentType;
|
|
vtkm::Vec<ValueType, 3> field;
|
|
inputField.CopyInto(field);
|
|
vtkm::Vec<WCoordType, 3> wpoints;
|
|
wCoords.CopyInto(wpoints);
|
|
return detail::TriangleDerivative(field, wpoints);
|
|
}
|
|
|
|
namespace detail
|
|
{
|
|
|
|
//-----------------------------------------------------------------------------
|
|
template <typename FieldVecType, typename WorldCoordType, typename ParametricCoordType>
|
|
VTKM_EXEC void PolygonSubTriangle(const FieldVecType& inField,
|
|
const WorldCoordType& inWCoords,
|
|
const vtkm::Vec<ParametricCoordType, 3>& pcoords,
|
|
vtkm::Vec<typename FieldVecType::ComponentType, 3>& outField,
|
|
vtkm::Vec<typename WorldCoordType::ComponentType, 3>& outWCoords,
|
|
const vtkm::exec::FunctorBase& worklet)
|
|
{
|
|
// 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.
|
|
|
|
// First point is right at pcoords
|
|
outField[0] = vtkm::exec::CellInterpolate(inField, pcoords, vtkm::CellShapeTagPolygon{}, worklet);
|
|
outWCoords[0] =
|
|
vtkm::exec::CellInterpolate(inWCoords, pcoords, vtkm::CellShapeTagPolygon{}, worklet);
|
|
|
|
// Find the unit vector pointing from the center of the polygon to pcoords
|
|
vtkm::Vec<ParametricCoordType, 2> radialVector = { pcoords[0] - 0.5, pcoords[1] - 0.5 };
|
|
ParametricCoordType magnitudeSquared = vtkm::MagnitudeSquared(radialVector);
|
|
if (magnitudeSquared > 8 * vtkm::Epsilon<ParametricCoordType>())
|
|
{
|
|
radialVector = vtkm::RSqrt(magnitudeSquared) * radialVector;
|
|
}
|
|
else
|
|
{
|
|
// pcoords is in the center of the polygon. Just point in an arbitrary direction
|
|
radialVector = { 1.0, 0.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 to
|
|
// get an accurate gradient.
|
|
static constexpr ParametricCoordType scale = 0.05;
|
|
|
|
vtkm::Vec<ParametricCoordType, 3> backPcoord = {
|
|
pcoords[0] + scale * (-radialVector[1] - radialVector[0]),
|
|
pcoords[1] + scale * (radialVector[0] - radialVector[1]),
|
|
0
|
|
};
|
|
outField[1] =
|
|
vtkm::exec::CellInterpolate(inField, backPcoord, vtkm::CellShapeTagPolygon{}, worklet);
|
|
outWCoords[1] =
|
|
vtkm::exec::CellInterpolate(inWCoords, backPcoord, vtkm::CellShapeTagPolygon{}, worklet);
|
|
|
|
backPcoord = { pcoords[0] + scale * (radialVector[1] - radialVector[0]),
|
|
pcoords[1] + scale * (-radialVector[0] - radialVector[1]),
|
|
0 };
|
|
outField[2] =
|
|
vtkm::exec::CellInterpolate(inField, backPcoord, vtkm::CellShapeTagPolygon{}, worklet);
|
|
outWCoords[2] =
|
|
vtkm::exec::CellInterpolate(inWCoords, backPcoord, vtkm::CellShapeTagPolygon{}, worklet);
|
|
}
|
|
|
|
} // namespace detail
|
|
|
|
//-----------------------------------------------------------------------------
|
|
template <typename FieldVecType, typename WorldCoordType, typename ParametricCoordType>
|
|
VTKM_EXEC vtkm::Vec<typename FieldVecType::ComponentType, 3> CellDerivative(
|
|
const FieldVecType& field,
|
|
const WorldCoordType& wCoords,
|
|
const vtkm::Vec<ParametricCoordType, 3>& pcoords,
|
|
vtkm::CellShapeTagPolygon,
|
|
const vtkm::exec::FunctorBase& worklet)
|
|
{
|
|
VTKM_ASSERT(field.GetNumberOfComponents() == wCoords.GetNumberOfComponents());
|
|
|
|
const vtkm::IdComponent numPoints = field.GetNumberOfComponents();
|
|
VTKM_ASSERT(numPoints > 0);
|
|
|
|
switch (field.GetNumberOfComponents())
|
|
{
|
|
case 1:
|
|
return CellDerivative(field, wCoords, pcoords, vtkm::CellShapeTagVertex(), worklet);
|
|
case 2:
|
|
return CellDerivative(field, wCoords, pcoords, vtkm::CellShapeTagLine(), worklet);
|
|
case 3:
|
|
return CellDerivative(field, wCoords, pcoords, vtkm::CellShapeTagTriangle(), worklet);
|
|
case 4:
|
|
return CellDerivative(field, wCoords, pcoords, vtkm::CellShapeTagQuad(), worklet);
|
|
}
|
|
|
|
vtkm::Vec<typename FieldVecType::ComponentType, 3> triangleField;
|
|
vtkm::Vec<typename WorldCoordType::ComponentType, 3> triangleWCoords;
|
|
detail::PolygonSubTriangle(field, wCoords, pcoords, triangleField, triangleWCoords, worklet);
|
|
return detail::TriangleDerivative(triangleField, triangleWCoords);
|
|
}
|
|
|
|
//-----------------------------------------------------------------------------
|
|
template <typename FieldVecType, typename WorldCoordType, typename ParametricCoordType>
|
|
VTKM_EXEC vtkm::Vec<typename FieldVecType::ComponentType, 3> CellDerivative(
|
|
const FieldVecType& inputField,
|
|
const WorldCoordType& wCoords,
|
|
const vtkm::Vec<ParametricCoordType, 3>& pcoords,
|
|
vtkm::CellShapeTagQuad,
|
|
const vtkm::exec::FunctorBase& vtkmNotUsed(worklet))
|
|
{
|
|
VTKM_ASSERT(inputField.GetNumberOfComponents() == 4);
|
|
VTKM_ASSERT(wCoords.GetNumberOfComponents() == 4);
|
|
|
|
using ValueType = typename FieldVecType::ComponentType;
|
|
using WCoordType = typename WorldCoordType::ComponentType;
|
|
vtkm::Vec<ValueType, 4> field;
|
|
inputField.CopyInto(field);
|
|
vtkm::Vec<WCoordType, 4> wpoints;
|
|
wCoords.CopyInto(wpoints);
|
|
|
|
return detail::CellDerivativeFor2DCell(field, wpoints, pcoords, vtkm::CellShapeTagQuad());
|
|
}
|
|
|
|
template <typename FieldVecType, typename ParametricCoordType>
|
|
VTKM_EXEC vtkm::Vec<typename FieldVecType::ComponentType, 3> CellDerivative(
|
|
const FieldVecType& field,
|
|
const vtkm::VecAxisAlignedPointCoordinates<2>& wCoords,
|
|
const vtkm::Vec<ParametricCoordType, 3>& pcoords,
|
|
vtkm::CellShapeTagQuad,
|
|
const vtkm::exec::FunctorBase& vtkmNotUsed(worklet))
|
|
{
|
|
VTKM_ASSERT(field.GetNumberOfComponents() == 4);
|
|
|
|
using T = typename FieldVecType::ComponentType;
|
|
using VecT = vtkm::Vec<T, 2>;
|
|
|
|
VecT pc(static_cast<typename vtkm::VecTraits<T>::ComponentType>(pcoords[0]),
|
|
static_cast<typename vtkm::VecTraits<T>::ComponentType>(pcoords[1]));
|
|
VecT rc = VecT(T(1)) - pc;
|
|
|
|
VecT sum = field[0] * VecT(-rc[1], -rc[0]);
|
|
sum = sum + field[1] * VecT(rc[1], -pc[0]);
|
|
sum = sum + field[2] * VecT(pc[1], pc[0]);
|
|
sum = sum + field[3] * VecT(-pc[1], rc[0]);
|
|
|
|
return vtkm::Vec<T, 3>(static_cast<T>(sum[0] / wCoords.GetSpacing()[0]),
|
|
static_cast<T>(sum[1] / wCoords.GetSpacing()[1]),
|
|
T(0));
|
|
}
|
|
|
|
//-----------------------------------------------------------------------------
|
|
namespace detail
|
|
{
|
|
|
|
template <typename ValueType, typename LUType>
|
|
VTKM_EXEC vtkm::Vec<ValueType, 3> TetraDerivativeFinish(
|
|
const vtkm::Vec<ValueType, 4>& field,
|
|
const vtkm::Matrix<LUType, 3, 3>& LUFactorization,
|
|
const vtkm::IdComponent3& LUPermutation,
|
|
vtkm::TypeTraitsScalarTag)
|
|
{
|
|
// Finish solving linear equation. See TriangleDerivative implementation
|
|
// for more detail.
|
|
vtkm::Vec<LUType, 3> b(field[1] - field[0], field[2] - field[0], field[3] - field[0]);
|
|
|
|
return vtkm::detail::MatrixLUPSolve(LUFactorization, LUPermutation, b);
|
|
}
|
|
|
|
template <typename ValueType, typename LUType>
|
|
VTKM_EXEC vtkm::Vec<ValueType, 3> TetraDerivativeFinish(
|
|
const vtkm::Vec<ValueType, 4>& field,
|
|
const vtkm::Matrix<LUType, 3, 3>& LUFactorization,
|
|
const vtkm::IdComponent3& LUPermutation,
|
|
vtkm::TypeTraitsVectorTag)
|
|
{
|
|
using FieldTraits = vtkm::VecTraits<ValueType>;
|
|
using FieldComponentType = typename FieldTraits::ComponentType;
|
|
|
|
vtkm::Vec<ValueType, 3> gradient(vtkm::TypeTraits<ValueType>::ZeroInitialization());
|
|
|
|
for (vtkm::IdComponent fieldComponent = 0;
|
|
fieldComponent < FieldTraits::GetNumberOfComponents(field[0]);
|
|
fieldComponent++)
|
|
{
|
|
vtkm::Vec<FieldComponentType, 4> subField(FieldTraits::GetComponent(field[0], fieldComponent),
|
|
FieldTraits::GetComponent(field[1], fieldComponent),
|
|
FieldTraits::GetComponent(field[2], fieldComponent),
|
|
FieldTraits::GetComponent(field[3], fieldComponent));
|
|
|
|
vtkm::Vec<FieldComponentType, 3> subGradient =
|
|
TetraDerivativeFinish(subField,
|
|
LUFactorization,
|
|
LUPermutation,
|
|
typename vtkm::TypeTraits<FieldComponentType>::DimensionalityTag());
|
|
|
|
FieldTraits::SetComponent(gradient[0], fieldComponent, subGradient[0]);
|
|
FieldTraits::SetComponent(gradient[1], fieldComponent, subGradient[1]);
|
|
FieldTraits::SetComponent(gradient[2], fieldComponent, subGradient[2]);
|
|
}
|
|
|
|
return gradient;
|
|
}
|
|
|
|
template <typename ValueType, typename LUType>
|
|
VTKM_EXEC vtkm::Vec<ValueType, 3> TetraDerivativeFinish(
|
|
const vtkm::Vec<ValueType, 4>& field,
|
|
const vtkm::Matrix<LUType, 3, 3>& LUFactorization,
|
|
const vtkm::IdComponent3& LUPermutation,
|
|
vtkm::TypeTraitsMatrixTag)
|
|
{
|
|
return TetraDerivativeFinish(field, LUFactorization, LUPermutation, vtkm::TypeTraitsVectorTag());
|
|
}
|
|
|
|
template <typename ValueType, typename WorldCoordType>
|
|
VTKM_EXEC vtkm::Vec<ValueType, 3> TetraDerivative(const vtkm::Vec<ValueType, 4>& field,
|
|
const vtkm::Vec<WorldCoordType, 4>& wCoords)
|
|
{
|
|
using BaseComponentType = typename BaseComponent<ValueType>::Type;
|
|
|
|
// The scalar values of the four points in a tetrahedron completely specify a
|
|
// linear field (with constant gradient). The field, defined by the 3-vector
|
|
// gradient g and scalar value s_origin, can be found with this set of 4
|
|
// equations and 4 unknowns.
|
|
//
|
|
// Dot(p0, g) + s_origin = s0
|
|
// Dot(p1, g) + s_origin = s1
|
|
// Dot(p2, g) + s_origin = s2
|
|
// Dot(p3, g) + s_origin = s3
|
|
//
|
|
// Where the p's are point coordinates. But we don't really care about
|
|
// s_origin. We just want to find the gradient g. With some simple
|
|
// elimination we, we can get rid of s_origin and be left with 3 equations
|
|
// and 3 unknowns.
|
|
//
|
|
// Dot(p1-p0, g) = s1 - s0
|
|
// Dot(p2-p0, g) = s2 - s0
|
|
// Dot(p3-p0, g) = s3 - s0
|
|
//
|
|
// We'll solve this by putting this in matrix form Ax = b where the rows of A
|
|
// are the differences in points and normal, b has the scalar differences,
|
|
// and x is really the gradient g.
|
|
|
|
vtkm::Vec<BaseComponentType, 3> v0 = wCoords[1] - wCoords[0];
|
|
vtkm::Vec<BaseComponentType, 3> v1 = wCoords[2] - wCoords[0];
|
|
vtkm::Vec<BaseComponentType, 3> v2 = wCoords[3] - wCoords[0];
|
|
|
|
vtkm::Matrix<BaseComponentType, 3, 3> A;
|
|
vtkm::MatrixSetRow(A, 0, v0);
|
|
vtkm::MatrixSetRow(A, 1, v1);
|
|
vtkm::MatrixSetRow(A, 2, v2);
|
|
|
|
// If the tetrahedron is degenerate, then valid will be false. For now we are
|
|
// ignoring it. We could detect it if we determine we need to although I have
|
|
// seen singular matrices missed due to floating point error.
|
|
//
|
|
bool valid;
|
|
|
|
// If you look at the implementation of vtkm::SolveLinearSystem, you will see
|
|
// that it is done in two parts. First, it does an LU factorization and then
|
|
// uses that result to complete the solve. The factorization part talkes the
|
|
// longest amount of time, and if we are performing the gradient on a vector
|
|
// field, the factorization can be reused for each component of the vector
|
|
// field. Thus, we are going to call the internals of SolveLinearSystem
|
|
// ourselves to do the factorization and then apply it to all components.
|
|
vtkm::IdComponent3 permutation;
|
|
BaseComponentType inversionParity; // Unused
|
|
vtkm::detail::MatrixLUPFactor(A, permutation, inversionParity, valid);
|
|
// MatrixLUPFactor does in place factorization. A is now the LU factorization.
|
|
return TetraDerivativeFinish(
|
|
field, A, permutation, typename vtkm::TypeTraits<ValueType>::DimensionalityTag());
|
|
}
|
|
|
|
} // namespace detail
|
|
|
|
//-----------------------------------------------------------------------------
|
|
template <typename FieldVecType, typename WorldCoordType, typename ParametricCoordType>
|
|
VTKM_EXEC vtkm::Vec<typename FieldVecType::ComponentType, 3> CellDerivative(
|
|
const FieldVecType& inputField,
|
|
const WorldCoordType& wCoords,
|
|
const vtkm::Vec<ParametricCoordType, 3>& vtkmNotUsed(pcoords),
|
|
vtkm::CellShapeTagTetra,
|
|
const vtkm::exec::FunctorBase& vtkmNotUsed(worklet))
|
|
{
|
|
VTKM_ASSERT(inputField.GetNumberOfComponents() == 4);
|
|
VTKM_ASSERT(wCoords.GetNumberOfComponents() == 4);
|
|
|
|
using ValueType = typename FieldVecType::ComponentType;
|
|
using WCoordType = typename WorldCoordType::ComponentType;
|
|
vtkm::Vec<ValueType, 4> field;
|
|
inputField.CopyInto(field);
|
|
vtkm::Vec<WCoordType, 4> wpoints;
|
|
wCoords.CopyInto(wpoints);
|
|
return detail::TetraDerivative(field, wpoints);
|
|
}
|
|
|
|
//-----------------------------------------------------------------------------
|
|
template <typename FieldVecType, typename WorldCoordType, typename ParametricCoordType>
|
|
VTKM_EXEC vtkm::Vec<typename FieldVecType::ComponentType, 3> CellDerivative(
|
|
const FieldVecType& inputField,
|
|
const WorldCoordType& wCoords,
|
|
const vtkm::Vec<ParametricCoordType, 3>& pcoords,
|
|
vtkm::CellShapeTagHexahedron,
|
|
const vtkm::exec::FunctorBase& vtkmNotUsed(worklet))
|
|
{
|
|
VTKM_ASSERT(inputField.GetNumberOfComponents() == 8);
|
|
VTKM_ASSERT(wCoords.GetNumberOfComponents() == 8);
|
|
|
|
using ValueType = typename FieldVecType::ComponentType;
|
|
using WCoordType = typename WorldCoordType::ComponentType;
|
|
vtkm::Vec<ValueType, 8> field;
|
|
inputField.CopyInto(field);
|
|
vtkm::Vec<WCoordType, 8> wpoints;
|
|
wCoords.CopyInto(wpoints);
|
|
|
|
return detail::CellDerivativeFor3DCell(field, wpoints, pcoords, vtkm::CellShapeTagHexahedron());
|
|
}
|
|
|
|
template <typename FieldVecType, typename ParametricCoordType>
|
|
VTKM_EXEC vtkm::Vec<typename FieldVecType::ComponentType, 3> CellDerivative(
|
|
const FieldVecType& field,
|
|
const vtkm::VecAxisAlignedPointCoordinates<3>& wCoords,
|
|
const vtkm::Vec<ParametricCoordType, 3>& pcoords,
|
|
vtkm::CellShapeTagHexahedron,
|
|
const vtkm::exec::FunctorBase& vtkmNotUsed(worklet))
|
|
{
|
|
VTKM_ASSERT(field.GetNumberOfComponents() == 8);
|
|
|
|
using T = typename FieldVecType::ComponentType;
|
|
using VecT = vtkm::Vec<T, 3>;
|
|
|
|
VecT pc(static_cast<typename vtkm::VecTraits<T>::ComponentType>(pcoords[0]),
|
|
static_cast<typename vtkm::VecTraits<T>::ComponentType>(pcoords[1]),
|
|
static_cast<typename vtkm::VecTraits<T>::ComponentType>(pcoords[2]));
|
|
VecT rc = VecT(T(1)) - pc;
|
|
|
|
VecT sum = field[0] * VecT(-rc[1] * rc[2], -rc[0] * rc[2], -rc[0] * rc[1]);
|
|
sum = sum + field[1] * VecT(rc[1] * rc[2], -pc[0] * rc[2], -pc[0] * rc[1]);
|
|
sum = sum + field[2] * VecT(pc[1] * rc[2], pc[0] * rc[2], -pc[0] * pc[1]);
|
|
sum = sum + field[3] * VecT(-pc[1] * rc[2], rc[0] * rc[2], -rc[0] * pc[1]);
|
|
sum = sum + field[4] * VecT(-rc[1] * pc[2], -rc[0] * pc[2], rc[0] * rc[1]);
|
|
sum = sum + field[5] * VecT(rc[1] * pc[2], -pc[0] * pc[2], pc[0] * rc[1]);
|
|
sum = sum + field[6] * VecT(pc[1] * pc[2], pc[0] * pc[2], pc[0] * pc[1]);
|
|
sum = sum + field[7] * VecT(-pc[1] * pc[2], rc[0] * pc[2], rc[0] * pc[1]);
|
|
|
|
return sum / wCoords.GetSpacing();
|
|
}
|
|
|
|
//-----------------------------------------------------------------------------
|
|
template <typename FieldVecType, typename WorldCoordType, typename ParametricCoordType>
|
|
VTKM_EXEC vtkm::Vec<typename FieldVecType::ComponentType, 3> CellDerivative(
|
|
const FieldVecType& inputField,
|
|
const WorldCoordType& wCoords,
|
|
const vtkm::Vec<ParametricCoordType, 3>& pcoords,
|
|
vtkm::CellShapeTagWedge,
|
|
const vtkm::exec::FunctorBase& vtkmNotUsed(worklet))
|
|
{
|
|
VTKM_ASSERT(inputField.GetNumberOfComponents() == 6);
|
|
VTKM_ASSERT(wCoords.GetNumberOfComponents() == 6);
|
|
|
|
using ValueType = typename FieldVecType::ComponentType;
|
|
using WCoordType = typename WorldCoordType::ComponentType;
|
|
vtkm::Vec<ValueType, 6> field;
|
|
inputField.CopyInto(field);
|
|
vtkm::Vec<WCoordType, 6> wpoints;
|
|
wCoords.CopyInto(wpoints);
|
|
|
|
return detail::CellDerivativeFor3DCell(field, wpoints, pcoords, vtkm::CellShapeTagWedge());
|
|
}
|
|
|
|
//-----------------------------------------------------------------------------
|
|
template <typename FieldVecType, typename WorldCoordType, typename ParametricCoordType>
|
|
VTKM_EXEC vtkm::Vec<typename FieldVecType::ComponentType, 3> CellDerivative(
|
|
const FieldVecType& inputField,
|
|
const WorldCoordType& wCoords,
|
|
const vtkm::Vec<ParametricCoordType, 3>& pcoords,
|
|
vtkm::CellShapeTagPyramid,
|
|
const vtkm::exec::FunctorBase& vtkmNotUsed(worklet))
|
|
{
|
|
VTKM_ASSERT(inputField.GetNumberOfComponents() == 5);
|
|
VTKM_ASSERT(wCoords.GetNumberOfComponents() == 5);
|
|
|
|
using ValueType = typename FieldVecType::ComponentType;
|
|
using WCoordType = typename WorldCoordType::ComponentType;
|
|
vtkm::Vec<ValueType, 5> field;
|
|
inputField.CopyInto(field);
|
|
vtkm::Vec<WCoordType, 5> wpoints;
|
|
wCoords.CopyInto(wpoints);
|
|
|
|
if (pcoords[2] > static_cast<ParametricCoordType>(0.999))
|
|
{
|
|
// 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.
|
|
using PCoordType = vtkm::Vec<ParametricCoordType, 3>;
|
|
using ReturnType = vtkm::Vec<ValueType, 3>;
|
|
PCoordType pcoords1(0.5f, 0.5f, (2 * 0.998f) - pcoords[2]);
|
|
ReturnType derivative1 =
|
|
detail::CellDerivativeFor3DCell(field, wpoints, pcoords1, vtkm::CellShapeTagPyramid{});
|
|
PCoordType pcoords2(0.5f, 0.5f, 0.998f);
|
|
ReturnType derivative2 =
|
|
detail::CellDerivativeFor3DCell(field, wpoints, pcoords2, vtkm::CellShapeTagPyramid{});
|
|
return (ValueType(2) * derivative2) - derivative1;
|
|
}
|
|
|
|
return detail::CellDerivativeFor3DCell(field, wpoints, pcoords, vtkm::CellShapeTagPyramid());
|
|
}
|
|
}
|
|
} // namespace vtkm::exec
|
|
|
|
#endif //vtk_m_exec_Derivative_h
|