Shuenhoy/vtk-m
1
0
Fork 0
mirror of https://gitlab.kitware.com/vtk/vtk-m synced 2024-09-08 13:23:51 +00:00
Code Issues 2 Actions Packages Projects Releases Wiki Activity

Special implementation of cell derivative for rectilinear cells

Browse Source
This commit is contained in:
Kenneth Moreland 2015-09-01 15:30:03 -07:00
parent 284fda03fd
commit b58543297a
2 changed files with 151 additions and 16 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

80
vtkm/exec/CellDerivative.h
View File

@ -23,6 +23,7 @@
#include <vtkm/CellShape.h>
#include <vtkm/Math.h>
#include <vtkm/Matrix.h>
#include <vtkm/VecRectilinearPointCoordinates.h>
#include <vtkm/VectorAnalysis.h>
#include <vtkm/exec/Assert.h>
@ -696,6 +697,24 @@ CellDerivative(const FieldVecType &field,
return (deltaField/vtkm::MagnitudeSquared(vec))*vec;
}
template<typename FieldVecType,
typename ParametricCoordType>
VTKM_EXEC_EXPORT
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 &worklet)
{
VTKM_ASSERT_EXEC(field.GetNumberOfComponents() == 2, worklet);
typedef typename FieldVecType::ComponentType T;
typedef vtkm::Vec<T,2> VecT;
return vtkm::Vec<T,3>((field[1]-field[0])/wCoords.GetSpacing()[0], 0, 0);
}
//-----------------------------------------------------------------------------
template<typename FieldVecType,
typename WorldCoordType,
@ -897,6 +916,35 @@ CellDerivative(const FieldVecType &field,
field, wCoords, pcoords, vtkm::CellShapeTagQuad());
}
template<typename FieldVecType,
typename ParametricCoordType>
VTKM_EXEC_EXPORT
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 &worklet)
{
VTKM_ASSERT_EXEC(field.GetNumberOfComponents() == 4, worklet);
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(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],
0);
}
//-----------------------------------------------------------------------------
template<typename FieldVecType,
typename WorldCoordType,
@ -984,6 +1032,38 @@ CellDerivative(const FieldVecType &field,
field, wCoords, pcoords, vtkm::CellShapeTagHexahedron());
}
template<typename FieldVecType,
typename ParametricCoordType>
VTKM_EXEC_EXPORT
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 &worklet)
{
VTKM_ASSERT_EXEC(field.GetNumberOfComponents() == 8, worklet);
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(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,

87
vtkm/exec/testing/UnitTestCellDerivative.cxx
View File

@ -93,11 +93,11 @@ void GetMinMaxPoints(CellShapeTag,
template<typename FieldType>
struct TestDerivativeFunctor
{
template<typename CellShapeTag>
void DoTest(CellShapeTag shape,
vtkm::IdComponent numPoints,
LinearField field,
vtkm::Vec<FieldType,3> expectedGradient) const
template<typename CellShapeTag, typename WCoordsVecType>
void DoTestWithWCoords(CellShapeTag shape,
const WCoordsVecType worldCoordinates,
LinearField field,
vtkm::Vec<FieldType,3> expectedGradient) const
{
// Stuff to fake running in the execution environment.
char messageBuffer[256];
@ -106,20 +106,12 @@ struct TestDerivativeFunctor
vtkm::exec::FunctorBase workletProxy;
workletProxy.SetErrorMessageBuffer(errorMessage);
vtkm::VecVariable<vtkm::Vec<vtkm::FloatDefault,3>, MAX_POINTS> worldCoordinates;
vtkm::IdComponent numPoints = worldCoordinates.GetNumberOfComponents();
vtkm::VecVariable<FieldType, MAX_POINTS> fieldValues;
for (vtkm::IdComponent pointIndex = 0; pointIndex < numPoints; pointIndex++)
{
vtkm::Vec<vtkm::FloatDefault,3> pcoords =
vtkm::exec::ParametricCoordinatesPoint(numPoints,
pointIndex,
shape,
workletProxy);
VTKM_TEST_ASSERT(!errorMessage.IsErrorRaised(), messageBuffer);
vtkm::Vec<vtkm::FloatDefault,3> wcoords = ParametricToWorld(pcoords);
VTKM_TEST_ASSERT(test_equal(pcoords, WorldToParametric(wcoords)),
"Test world/parametric conversion broken.");
worldCoordinates.Append(wcoords);
vtkm::Vec<vtkm::FloatDefault,3> wcoords = worldCoordinates[pointIndex];
FieldType value = field.GetValue(wcoords);
fieldValues.Append(value);
}
@ -168,6 +160,37 @@ struct TestDerivativeFunctor
}
}
template<typename CellShapeTag>
void DoTest(CellShapeTag shape,
vtkm::IdComponent numPoints,
LinearField field,
vtkm::Vec<FieldType,3> expectedGradient) const
{
// Stuff to fake running in the execution environment.
char messageBuffer[256];
messageBuffer[0] = '\0';
vtkm::exec::internal::ErrorMessageBuffer errorMessage(messageBuffer, 256);
vtkm::exec::FunctorBase workletProxy;
workletProxy.SetErrorMessageBuffer(errorMessage);
vtkm::VecVariable<vtkm::Vec<vtkm::FloatDefault,3>, MAX_POINTS> worldCoordinates;
for (vtkm::IdComponent pointIndex = 0; pointIndex < numPoints; pointIndex++)
{
vtkm::Vec<vtkm::FloatDefault,3> pcoords =
vtkm::exec::ParametricCoordinatesPoint(numPoints,
pointIndex,
shape,
workletProxy);
VTKM_TEST_ASSERT(!errorMessage.IsErrorRaised(), messageBuffer);
vtkm::Vec<vtkm::FloatDefault,3> wcoords = ParametricToWorld(pcoords);
VTKM_TEST_ASSERT(test_equal(pcoords, WorldToParametric(wcoords)),
"Test world/parametric conversion broken.");
worldCoordinates.Append(wcoords);
}
this->DoTestWithWCoords(shape, worldCoordinates, field, expectedGradient);
}
template<typename CellShapeTag>
void DoTest(CellShapeTag shape,
vtkm::IdComponent numPoints,
@ -266,6 +289,38 @@ void TestDerivative()
vtkm::testing::Testing::TryAllCellShapes(TestDerivativeFunctor<vtkm::Float32>());
std::cout << "======== Float64 ==========================" << std::endl;
vtkm::testing::Testing::TryAllCellShapes(TestDerivativeFunctor<vtkm::Float64>());
boost::random::uniform_real_distribution<vtkm::Float64> randomDist(-20.0, 20.0);
LinearField field;
field.OriginValue = randomDist(g_RandomGenerator);
field.Gradient = vtkm::make_Vec(randomDist(g_RandomGenerator),
randomDist(g_RandomGenerator),
randomDist(g_RandomGenerator));
vtkm::Vec<vtkm::Float64,3> expectedGradient = field.Gradient;
TestDerivativeFunctor<vtkm::Float64> testFunctor;
vtkm::Vec<vtkm::FloatDefault,3> origin = vtkm::Vec<vtkm::FloatDefault,3>(0.25f, 0.25f, 0.25f);
vtkm::Vec<vtkm::FloatDefault,3> spacing = vtkm::Vec<vtkm::FloatDefault,3>(2.0f, 2.0f, 2.0f);
std::cout << "======== Uniform Point Coordinates 3D =====" << std::endl;
testFunctor.DoTestWithWCoords(
vtkm::CellShapeTagHexahedron(),
vtkm::VecRectilinearPointCoordinates<3>(origin, spacing),
field,
expectedGradient);
std::cout << "======== Uniform Point Coordinates 2D =====" << std::endl;
expectedGradient[2] = 0.0;
testFunctor.DoTestWithWCoords(
vtkm::CellShapeTagQuad(),
vtkm::VecRectilinearPointCoordinates<2>(origin, spacing),
field,
expectedGradient);
std::cout << "======== Uniform Point Coordinates 1D =====" << std::endl;
expectedGradient[1] = 0.0;
testFunctor.DoTestWithWCoords(
vtkm::CellShapeTagLine(),
vtkm::VecRectilinearPointCoordinates<1>(origin, spacing),
field,
expectedGradient);
}
} // anonymous namespace