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

Fix gradient issue at apex of pyramid cells

Browse Source 
The gradient is malformed at the apex of a pyramid. To get around this,
steal a trick from the VTK source where in this case interpolate some
values a little bit into the interior of the cell.

Also expand the gradient worklet tests to include all cell types.
This commit is contained in:
Kenneth Moreland 2019-08-19 14:21:42 -06:00
parent 484c587221
commit c19b312e58
2 changed files with 65 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

21
vtkm/exec/CellDerivative.h
View File

@ -1165,6 +1165,27 @@ VTKM_EXEC vtkm::Vec<typename FieldVecType::ComponentType, 3> CellDerivative(
vtkm::Vec<WCoordType, 5> wpoints;
wCoords.CopyInto(wpoints);
if (pcoords[2] > static_cast<ParametricCoordType>(0.999))
{
// If we are at the apex of the pyramid we need to do something special.
// As we approach the apex, the derivatives of the parametric shape
// functions in x and y go to 0 while the inverse of the Jacobian
// also goes to 0. This results in 0/0 but using l'Hopital's rule
// we could actually compute the value of the limit, if we had a
// functional expression to compute the gradient. We're on a computer
// so we don't but we can cheat and do a linear extrapolation of the
// derivatives which really ends up as the same thing.
using PCoordType = vtkm::Vec<ParametricCoordType, 3>;
using ReturnType = vtkm::Vec<ValueType, 3>;
PCoordType pcoords1(0.5f, 0.5f, (2 * 0.998f) - pcoords[2]);
ReturnType derivative1 =
detail::CellDerivativeFor3DCell(field, wpoints, pcoords1, vtkm::CellShapeTagPyramid{});
PCoordType pcoords2(0.5f, 0.5f, 0.998f);
ReturnType derivative2 =
detail::CellDerivativeFor3DCell(field, wpoints, pcoords2, vtkm::CellShapeTagPyramid{});
return (ValueType(2) * derivative2) - derivative1;
}
return detail::CellDerivativeFor3DCell(field, wpoints, pcoords, vtkm::CellShapeTagPyramid());
}
}

50
vtkm/worklet/testing/UnitTestPointGradient.cxx
View File

@ -182,12 +182,12 @@ void TestPointGradientUniform3DWithVectorField2()
}
}
void TestPointGradientExplicit()
void TestPointGradientExplicit3D()
{
std::cout << "Testing PointGradient Worklet on Explicit data" << std::endl;
std::cout << "Testing PointGradient Worklet on Explicit 3D data" << std::endl;
vtkm::cont::testing::MakeTestDataSet testDataSet;
vtkm::cont::DataSet dataSet = testDataSet.Make3DExplicitDataSet0();
vtkm::cont::DataSet dataSet = testDataSet.Make3DExplicitDataSet5();
vtkm::cont::ArrayHandle<vtkm::Float32> fieldArray;
dataSet.GetField("pointvar").GetData().CopyTo(fieldArray);
@ -195,21 +195,59 @@ void TestPointGradientExplicit()
vtkm::worklet::PointGradient gradient;
auto result = gradient.Run(dataSet.GetCellSet(), dataSet.GetCoordinateSystem(), fieldArray);
vtkm::Vec3f_32 expected[2] = { { 10.f, 10.1f, 0.0f }, { 10.f, 10.1f, 0.0f } };
for (int i = 0; i < 2; ++i)
//vtkm::cont::printSummary_ArrayHandle(result, std::cout, true);
const int nVerts = 11;
vtkm::Vec3f_32 expected[nVerts] = {
{ 10.0f, 40.2f, 30.1f }, { 27.4f, 40.1f, 10.1f }, { 17.425f, 40.0f, 10.1f },
{ -10.0f, 40.1f, 30.1f }, { 9.9f, -0.0500011f, 30.0f }, { 16.2125f, -4.55f, 10.0f },
{ 6.2f, -4.6f, 10.0f }, { -10.1f, -0.0999985f, 30.0f }, { 22.5125f, -4.575f, 10.025f },
{ 1.0f, -40.3f, 30.0f }, { 0.6f, -49.2f, 10.0f }
};
for (int i = 0; i < nVerts; ++i)
{
VTKM_TEST_ASSERT(test_equal(result.GetPortalConstControl().Get(i), expected[i]),
"Wrong result for PointGradient worklet on 3D explicit data");
}
}
void TestPointGradientExplicit2D()
{
std::cout << "Testing PointGradient Worklet on Explicit 3D data" << std::endl;
vtkm::cont::testing::MakeTestDataSet testDataSet;
vtkm::cont::DataSet dataSet = testDataSet.Make2DExplicitDataSet0();
vtkm::cont::ArrayHandle<vtkm::Float32> fieldArray;
dataSet.GetField("pointvar").GetData().CopyTo(fieldArray);
vtkm::worklet::PointGradient gradient;
auto result = gradient.Run(dataSet.GetCellSet(), dataSet.GetCoordinateSystem(), fieldArray);
//vtkm::cont::printSummary_ArrayHandle(result, std::cout, true);
const int nVerts = 16;
vtkm::Vec3f_32 expected[nVerts] = { { -22.0f, -7.0f, 0.0f }, { -25.5f, -7.0f, 0.0f },
{ -30.5f, 7.0f, 0.0f }, { -32.0f, 16.0f, 0.0f },
{ -23.0f, -42.0f, 0.0f }, { -23.25f, -17.0f, 0.0f },
{ -20.6667f, 1.33333f, 0.0f }, { -23.0f, 14.0f, 0.0f },
{ -8.0f, -42.0f, 0.0f }, { 1.33333f, -19.7222f, 0.0f },
{ -6.2963f, -4.03704f, 0.0f }, { -5.0f, 12.0f, 0.0f },
{ 22.5556f, -13.4444f, 0.0f }, { -19.8889f, 59.1111f, 0.0f },
{ 22.5556f, 15.4444f, 0.0f }, { 45.0f, 26.6667f, 0.0f } };
for (int i = 0; i < nVerts; ++i)
{
VTKM_TEST_ASSERT(test_equal(result.GetPortalConstControl().Get(i), expected[i]),
"Wrong result for PointGradient worklet on 2D explicit data");
}
}
void TestPointGradient()
{
TestPointGradientUniform2D();
TestPointGradientUniform3D();
TestPointGradientUniform3DWithVectorField();
TestPointGradientUniform3DWithVectorField2();
TestPointGradientExplicit();
TestPointGradientExplicit2D();
TestPointGradientExplicit3D();
}
}