mirror of
https://gitlab.kitware.com/vtk/vtk-m
synced 2024-09-16 17:22:55 +00:00
43da67f747
When the implementation for the cell derivative functions was made, the special cases for creating the Jacobian for wedge and pyramid cell shapes was not working. Instead, it just used the hexahedra case for a degenerate cell. This fixes the issues with the special cases. There were multiple issues to be fixed. There were some complaints by the compiler about types. There was a mistake in the pyramid table. But the biggest issue was a problem with macro expansions. It was the classic tale of forgetting that you have to encase parameters in parenthesis to make sure that operator precedence comes out as expected.
1067 lines
43 KiB
C++
1067 lines
43 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 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
|
|
// Copyright 2015 UT-Battelle, LLC.
|
|
// Copyright 2015 Los Alamos National Security.
|
|
//
|
|
// Under the terms of Contract DE-NA0003525 with NTESS,
|
|
// the U.S. Government retains certain rights in this software.
|
|
//
|
|
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
|
|
// Laboratory (LANL), the U.S. Government retains certain rights in
|
|
// this software.
|
|
//============================================================================
|
|
#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 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)
|
|
{
|
|
using FieldType = typename FieldVecType::ComponentType;
|
|
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, 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(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::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 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::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>& 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>((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(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::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 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::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.
|
|
|
|
using FieldType = typename FieldVecType::ComponentType;
|
|
using WCoordType = typename WorldCoordType::ComponentType;
|
|
|
|
// 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::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<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(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::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 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::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::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<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
|