mirror of
https://gitlab.kitware.com/vtk/vtk-m
synced 2024-09-16 17:22:55 +00:00
1263 lines
47 KiB
C++
1263 lines
47 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_Derivative_h
|
|
#define vtk_m_exec_Derivative_h
|
|
|
|
#include <vtkm/Assert.h>
|
|
#include <vtkm/CellShape.h>
|
|
#include <vtkm/Math.h>
|
|
#include <vtkm/Matrix.h>
|
|
#include <vtkm/VecRectilinearPointCoordinates.h>
|
|
#include <vtkm/VectorAnalysis.h>
|
|
|
|
#include <vtkm/exec/CellInterpolate.h>
|
|
#include <vtkm/exec/Jacobian.h>
|
|
#include <vtkm/exec/FunctorBase.h>
|
|
|
|
namespace vtkm {
|
|
namespace exec {
|
|
|
|
namespace detail {
|
|
|
|
template<typename VecType, typename DimensionalityTag>
|
|
struct BaseComponentOfImpl;
|
|
|
|
template<typename VecType>
|
|
struct BaseComponentOfImpl<VecType, vtkm::TypeTraitsVectorTag>
|
|
{
|
|
private:
|
|
using ComponentType = typename vtkm::VecTraits<VecType>::ComponentType;
|
|
public:
|
|
using type =
|
|
typename BaseComponentOfImpl<
|
|
ComponentType,
|
|
typename vtkm::TypeTraits<ComponentType>::DimensionalityTag
|
|
>::type;
|
|
};
|
|
|
|
template<typename VecType>
|
|
struct BaseComponentOfImpl<VecType, vtkm::TypeTraitsMatrixTag>
|
|
: BaseComponentOfImpl<VecType, vtkm::TypeTraitsVectorTag>
|
|
{ };
|
|
|
|
template<typename ScalarType>
|
|
struct BaseComponentOfImpl<ScalarType, vtkm::TypeTraitsScalarTag>
|
|
{
|
|
using type = ScalarType;
|
|
};
|
|
|
|
} // namespace detail
|
|
|
|
/// Finds the base component type of a Vec. If you have a Vec of Vecs, it will
|
|
/// descend all Vecs until you get to the scalar type.
|
|
///
|
|
// If this becomes useful outside of CellDerivative, then this should probably
|
|
// be moved to a different header file.
|
|
template<typename VecType>
|
|
struct BaseComponentOf
|
|
{
|
|
using type =
|
|
typename detail::BaseComponentOfImpl<
|
|
VecType,
|
|
typename vtkm::TypeTraits<VecType>::DimensionalityTag>::type;
|
|
};
|
|
|
|
// 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)
|
|
{
|
|
typedef typename FieldVecType::ComponentType FieldType;
|
|
typedef vtkm::Vec<FieldType,3> GradientType;
|
|
|
|
GradientType pc(pcoords);
|
|
GradientType rc = GradientType(FieldType(1)) - pc;
|
|
|
|
GradientType parametricDerivative(FieldType(0));
|
|
VTKM_ACCUM_PARAMETRIC_DERIVATIVE_3D(0, -rc[1]*rc[2], -rc[0]*rc[2], -rc[0]*rc[1]);
|
|
VTKM_ACCUM_PARAMETRIC_DERIVATIVE_3D(1, rc[1]*rc[2], -pc[0]*rc[2], -pc[0]*rc[1]);
|
|
VTKM_ACCUM_PARAMETRIC_DERIVATIVE_3D(2, pc[1]*rc[2], pc[0]*rc[2], -pc[0]*pc[1]);
|
|
VTKM_ACCUM_PARAMETRIC_DERIVATIVE_3D(3, -pc[1]*rc[2], rc[0]*rc[2], -rc[0]*pc[1]);
|
|
VTKM_ACCUM_PARAMETRIC_DERIVATIVE_3D(4, -rc[1]*pc[2], -rc[0]*pc[2], rc[0]*rc[1]);
|
|
VTKM_ACCUM_PARAMETRIC_DERIVATIVE_3D(5, rc[1]*pc[2], -pc[0]*pc[2], pc[0]*rc[1]);
|
|
VTKM_ACCUM_PARAMETRIC_DERIVATIVE_3D(6, pc[1]*pc[2], pc[0]*pc[2], pc[0]*pc[1]);
|
|
VTKM_ACCUM_PARAMETRIC_DERIVATIVE_3D(7, -pc[1]*pc[2], rc[0]*pc[2], rc[0]*pc[1]);
|
|
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)
|
|
{
|
|
#if 0
|
|
// This is not working. Just leverage the hexahedron code that is working.
|
|
typedef typename FieldVecType::ComponentType FieldType;
|
|
typedef vtkm::Vec<FieldType,3> GradientType;
|
|
|
|
GradientType pc(pcoords);
|
|
GradientType rc = GradientType(1) - pc;
|
|
|
|
GradientType parametricDerivative(0);
|
|
VTKM_DERIVATIVE_WEIGHTS_WEDGE(pc, rc, VTKM_ACCUM_PARAMETRIC_DERIVATIVE_3D);
|
|
|
|
return parametricDerivative;
|
|
#else
|
|
return ParametricDerivative(vtkm::exec::internal::PermuteWedgeToHex(field),
|
|
pcoords,
|
|
vtkm::CellShapeTagHexahedron());
|
|
#endif
|
|
}
|
|
|
|
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)
|
|
{
|
|
#if 0
|
|
// This is not working. Just leverage the hexahedron code that is working.
|
|
typedef typename FieldVecType::ComponentType FieldType;
|
|
typedef vtkm::Vec<FieldType,3> GradientType;
|
|
|
|
GradientType pc(pcoords);
|
|
GradientType rc = GradientType(1) - pc;
|
|
|
|
GradientType parametricDerivative(0);
|
|
VTKM_DERIVATIVE_WEIGHTS_PYRAMID(pc, rc, VTKM_ACCUM_PARAMETRIC_DERIVATIVE_3D);
|
|
|
|
return parametricDerivative;
|
|
#else
|
|
return ParametricDerivative(vtkm::exec::internal::PermutePyramidToHex(field),
|
|
pcoords,
|
|
vtkm::CellShapeTagHexahedron());
|
|
#endif
|
|
}
|
|
|
|
#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)
|
|
{
|
|
typedef typename FieldVecType::ComponentType FieldType;
|
|
typedef vtkm::Vec<FieldType,2> GradientType;
|
|
|
|
GradientType pc(static_cast<FieldType>(pcoords[0]),
|
|
static_cast<FieldType>(pcoords[1]));
|
|
GradientType rc = GradientType(FieldType(1)) - pc;
|
|
|
|
GradientType parametricDerivative(FieldType(0));
|
|
VTKM_ACCUM_PARAMETRIC_DERIVATIVE_2D(0, -rc[1], -rc[0]);
|
|
VTKM_ACCUM_PARAMETRIC_DERIVATIVE_2D(1, rc[1], -pc[0]);
|
|
VTKM_ACCUM_PARAMETRIC_DERIVATIVE_2D(2, pc[1], pc[0]);
|
|
VTKM_ACCUM_PARAMETRIC_DERIVATIVE_2D(3, -pc[1], rc[0]);
|
|
|
|
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)
|
|
{
|
|
typedef typename FieldVecType::ComponentType FieldType;
|
|
typedef vtkm::Vec<FieldType,2> GradientType;
|
|
|
|
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 FieldVecType,
|
|
typename WorldCoordType,
|
|
typename ParametricCoordType,
|
|
typename CellShapeTag>
|
|
VTKM_EXEC
|
|
vtkm::Vec<typename FieldVecType::ComponentType,3>
|
|
CellDerivativeFor3DCell(const FieldVecType &field,
|
|
const WorldCoordType &wCoords,
|
|
const vtkm::Vec<ParametricCoordType,3> &pcoords,
|
|
CellShapeTag)
|
|
{
|
|
typedef typename FieldVecType::ComponentType FieldType;
|
|
typedef vtkm::Vec<FieldType,3> GradientType;
|
|
|
|
// 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, CellShapeTag());
|
|
jacobianTranspose = vtkm::MatrixTranspose(jacobianTranspose);
|
|
|
|
GradientType 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.
|
|
//
|
|
// | | | | | |
|
|
// | 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 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::Vec<vtkm::IdComponent,2> &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::Vec<vtkm::IdComponent,2> &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;
|
|
|
|
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::Vec<vtkm::IdComponent,2> &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 BaseComponentOf<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::Vec<vtkm::IdComponent,2> 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> ¶metricCoords,
|
|
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);
|
|
|
|
typedef vtkm::Vec<typename FieldVecType::ComponentType,3> GradientType;
|
|
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;
|
|
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 BaseComponentOf<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::VecRectilinearPointCoordinates<1> &wCoords,
|
|
const vtkm::Vec<ParametricCoordType,3> &vtkmNotUsed(pcoords),
|
|
vtkm::CellShapeTagLine,
|
|
const vtkm::exec::FunctorBase &vtkmNotUsed(worklet))
|
|
{
|
|
VTKM_ASSERT(field.GetNumberOfComponents() == 2);
|
|
|
|
typedef typename FieldVecType::ComponentType T;
|
|
|
|
return vtkm::Vec<T,3>(
|
|
(field[1]-field[0])/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::Vec<vtkm::IdComponent,3> &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::Vec<vtkm::IdComponent,3> &LUPermutation,
|
|
vtkm::TypeTraitsVectorTag)
|
|
{
|
|
using FieldTraits = vtkm::VecTraits<ValueType>;
|
|
using FieldComponentType = typename FieldTraits::ComponentType;
|
|
|
|
vtkm::Vec<ValueType,3> gradient;
|
|
|
|
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::Vec<vtkm::IdComponent,3> &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 BaseComponentOf<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::Vec<vtkm::IdComponent,3> 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);
|
|
}
|
|
|
|
//-----------------------------------------------------------------------------
|
|
template <typename ParametricCoordType>
|
|
VTKM_EXEC void
|
|
PolygonComputeIndices(const vtkm::Vec<ParametricCoordType, 3>& pcoords,
|
|
vtkm::IdComponent numPoints, vtkm::IdComponent& firstPointIndex,
|
|
vtkm::IdComponent& secondPointIndex)
|
|
{
|
|
ParametricCoordType angle;
|
|
if ((vtkm::Abs(pcoords[0]-0.5f) < 4*vtkm::Epsilon<ParametricCoordType>()) &&
|
|
(vtkm::Abs(pcoords[1]-0.5f) < 4*vtkm::Epsilon<ParametricCoordType>()))
|
|
{
|
|
angle = 0;
|
|
}
|
|
else
|
|
{
|
|
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);
|
|
firstPointIndex =
|
|
static_cast<vtkm::IdComponent>(vtkm::Floor(angle/deltaAngle));
|
|
secondPointIndex = firstPointIndex + 1;
|
|
if (secondPointIndex == numPoints)
|
|
{
|
|
secondPointIndex = 0;
|
|
}
|
|
}
|
|
|
|
//-----------------------------------------------------------------------------
|
|
template<typename FieldVecType,
|
|
typename WorldCoordType>
|
|
VTKM_EXEC
|
|
vtkm::Vec<typename FieldVecType::ComponentType,3>
|
|
PolygonDerivative(const FieldVecType &field,
|
|
const WorldCoordType &wCoords,
|
|
vtkm::IdComponent numPoints,
|
|
vtkm::IdComponent firstPointIndex,
|
|
vtkm::IdComponent secondPointIndex)
|
|
{
|
|
// 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;
|
|
typedef typename WorldCoordType::ComponentType WCoordType;
|
|
|
|
// Find the interpolation for the center point.
|
|
FieldType fieldCenter = field[0];
|
|
WCoordType wcoordCenter = wCoords[0];
|
|
for (vtkm::IdComponent pointIndex = 1; pointIndex < numPoints; pointIndex++)
|
|
{
|
|
fieldCenter = fieldCenter + field[pointIndex];
|
|
wcoordCenter = wcoordCenter + wCoords[pointIndex];
|
|
}
|
|
fieldCenter = fieldCenter*FieldType(1.0f/static_cast<float>(numPoints));
|
|
wcoordCenter = wcoordCenter*WCoordType(1.0f/static_cast<float>(numPoints));
|
|
|
|
// Set up parameters for triangle that pcoords is in
|
|
vtkm::Vec<FieldType,3> triangleField( fieldCenter,
|
|
field[firstPointIndex],
|
|
field[secondPointIndex]);
|
|
|
|
vtkm::Vec<WCoordType,3> triangleWCoords( wcoordCenter,
|
|
wCoords[firstPointIndex],
|
|
wCoords[secondPointIndex] );
|
|
|
|
// Now use the triangle derivative. pcoords is actually invalid for the
|
|
// triangle, but that does not matter as the derivative for a triangle does
|
|
// not depend on it.
|
|
return detail::TriangleDerivative(triangleField, triangleWCoords);
|
|
}
|
|
|
|
//-----------------------------------------------------------------------------
|
|
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::IdComponent firstPointIndex, secondPointIndex;
|
|
PolygonComputeIndices(pcoords,numPoints,firstPointIndex,secondPointIndex);
|
|
return PolygonDerivative(field, wCoords, numPoints, firstPointIndex, secondPointIndex);
|
|
}
|
|
|
|
//-----------------------------------------------------------------------------
|
|
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::VecRectilinearPointCoordinates<2> &wCoords,
|
|
const vtkm::Vec<ParametricCoordType,3> &pcoords,
|
|
vtkm::CellShapeTagQuad,
|
|
const vtkm::exec::FunctorBase &vtkmNotUsed(worklet))
|
|
{
|
|
VTKM_ASSERT(field.GetNumberOfComponents() == 4);
|
|
|
|
typedef typename FieldVecType::ComponentType T;
|
|
typedef vtkm::Vec<T,2> VecT;
|
|
|
|
VecT pc(static_cast<T>(pcoords[0]),
|
|
static_cast<T>(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>(sum[0]/wCoords.GetSpacing()[0],
|
|
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::Vec<vtkm::IdComponent,3> &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::Vec<vtkm::IdComponent,3> &LUPermutation,
|
|
vtkm::TypeTraitsVectorTag)
|
|
{
|
|
using FieldTraits = vtkm::VecTraits<ValueType>;
|
|
using FieldComponentType = typename FieldTraits::ComponentType;
|
|
|
|
vtkm::Vec<ValueType,3> gradient;
|
|
|
|
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::Vec<vtkm::IdComponent,3> &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 BaseComponentOf<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::Vec<vtkm::IdComponent,3> 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::VecRectilinearPointCoordinates<3> &wCoords,
|
|
const vtkm::Vec<ParametricCoordType,3> &pcoords,
|
|
vtkm::CellShapeTagHexahedron,
|
|
const vtkm::exec::FunctorBase &vtkmNotUsed(worklet))
|
|
{
|
|
VTKM_ASSERT(field.GetNumberOfComponents() == 8);
|
|
|
|
typedef typename FieldVecType::ComponentType T;
|
|
typedef vtkm::Vec<T,3> VecT;
|
|
|
|
VecT pc(static_cast<T>(pcoords[0]),
|
|
static_cast<T>(pcoords[1]),
|
|
static_cast<T>(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);
|
|
|
|
return detail::CellDerivativeFor3DCell(
|
|
field, wpoints, pcoords, vtkm::CellShapeTagPyramid());
|
|
}
|
|
|
|
|
|
}
|
|
} // namespace vtkm::exec
|
|
|
|
#endif //vtk_m_exec_Derivative_h
|