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

Optimize StructuredPointGradient for non boundary points.

Browse Source 
We can remove lots of clamp calls by checking if we are on
a boundary and using a non-clamped API.
This commit is contained in:
Robert Maynard 2020-04-21 14:05:30 -04:00
parent 3074b6a149
commit 93d87e06fd
3 changed files with 128 additions and 23 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

36
vtkm/exec/BoundaryState.h
View File

@ -176,6 +176,24 @@ struct BoundaryState
}
//@}
//@{
/// Takes a local neighborhood index (in the ranges of -neighborhood size to neighborhood size)
/// and returns the ijk of the equivalent point in the full data set. If the given value is out
/// of range, the returned value is undefined.
///
VTKM_EXEC vtkm::Id3 NeighborIndexToFullIndex(const vtkm::IdComponent3& neighbor) const
{
return this->IJK + neighbor;
}
VTKM_EXEC vtkm::Id3 NeighborIndexToFullIndex(vtkm::IdComponent neighborI,
vtkm::IdComponent neighborJ,
vtkm::IdComponent neighborK) const
{
return this->NeighborIndexToFullIndex(vtkm::make_Vec(neighborI, neighborJ, neighborK));
}
//@}
//@{
/// Takes a local neighborhood index (in the ranges of -neighborhood size to
/// neighborhood size), clamps it to the dataset bounds, and returns a new
@ -221,6 +239,24 @@ struct BoundaryState
}
//@}
//@{
/// Takes a local neighborhood index (in the ranges of -neighborhood size to neighborhood size)
/// and returns the flat index of the equivalent point in the full data set. If the given value
/// is out of range, the result is undefined.
///
VTKM_EXEC vtkm::Id NeighborIndexToFlatIndex(const vtkm::IdComponent3& neighbor) const
{
vtkm::Id3 full = this->IJK + neighbor;
return (full[2] * this->PointDimensions[1] + full[1]) * this->PointDimensions[0] + full[0];
}
VTKM_EXEC vtkm::Id NeighborIndexToFlatIndex(vtkm::IdComponent neighborI,
vtkm::IdComponent neighborJ,
vtkm::IdComponent neighborK) const
{
return this->NeighborIndexToFlatIndex(vtkm::make_Vec(neighborI, neighborJ, neighborK));
}
//@}
vtkm::Id3 IJK;
vtkm::Id3 PointDimensions;
};

24
vtkm/exec/FieldNeighborhood.h
View File

@ -50,12 +50,24 @@ struct FieldNeighborhood
return Portal.Get(this->Boundary->NeighborIndexToFlatIndexClamp(i, j, k));
}
VTKM_EXEC
ValueType GetUnchecked(vtkm::IdComponent i, vtkm::IdComponent j, vtkm::IdComponent k) const
{
return Portal.Get(this->Boundary->NeighborIndexToFlatIndex(i, j, k));
}
VTKM_EXEC
ValueType Get(const vtkm::Id3& ijk) const
{
return Portal.Get(this->Boundary->NeighborIndexToFlatIndexClamp(ijk));
}
VTKM_EXEC
ValueType GetUnchecked(const vtkm::Id3& ijk) const
{
return Portal.Get(this->Boundary->NeighborIndexToFlatIndex(ijk));
}
vtkm::exec::BoundaryState const* const Boundary;
FieldPortalType Portal;
};
@ -82,12 +94,24 @@ struct FieldNeighborhood<vtkm::internal::ArrayPortalUniformPointCoordinates>
return Portal.Get(this->Boundary->NeighborIndexToFullIndexClamp(i, j, k));
}
VTKM_EXEC
ValueType GetUnchecked(vtkm::IdComponent i, vtkm::IdComponent j, vtkm::IdComponent k) const
{
return Portal.Get(this->Boundary->NeighborIndexToFullIndex(i, j, k));
}
VTKM_EXEC
ValueType Get(const vtkm::IdComponent3& ijk) const
{
return Portal.Get(this->Boundary->NeighborIndexToFullIndexClamp(ijk));
}
VTKM_EXEC
ValueType GetUnchecked(const vtkm::IdComponent3& ijk) const
{
return Portal.Get(this->Boundary->NeighborIndexToFullIndex(ijk));
}
vtkm::exec::BoundaryState const* const Boundary;
vtkm::internal::ArrayPortalUniformPointCoordinates Portal;
};

91
vtkm/worklet/gradient/StructuredPointGradient.h
View File

@ -45,15 +45,19 @@ struct StructuredPointGradient : public vtkm::worklet::WorkletPointNeighborhood
using OT = typename GradientOutType::ComponentType;
vtkm::Vec<CT, 3> xi, eta, zeta;
this->Jacobian(inputPoints, boundary, xi, eta, zeta); //store the metrics in xi,eta,zeta
vtkm::Vec<bool, 3> onBoundary{ !boundary.IsRadiusInXBoundary(1),
!boundary.IsRadiusInYBoundary(1),
!boundary.IsRadiusInZBoundary(1) };
this->Jacobian(inputPoints, onBoundary, xi, eta, zeta); //store the metrics in xi,eta,zeta
auto dxi = inputField.Get(1, 0, 0) - inputField.Get(-1, 0, 0);
auto deta = inputField.Get(0, 1, 0) - inputField.Get(0, -1, 0);
auto dzeta = inputField.Get(0, 0, 1) - inputField.Get(0, 0, -1);
dxi = (boundary.IsRadiusInXBoundary(1) ? dxi * 0.5f : dxi);
deta = (boundary.IsRadiusInYBoundary(1) ? deta * 0.5f : deta);
dzeta = (boundary.IsRadiusInZBoundary(1) ? dzeta * 0.5f : dzeta);
dxi = (onBoundary[0] ? dxi : dxi * 0.5f);
deta = (onBoundary[1] ? deta : deta * 0.5f);
dzeta = (onBoundary[2] ? dzeta : dzeta * 0.5f);
outputGradient[0] = static_cast<OT>(xi[0] * dxi + eta[0] * deta + zeta[0] * dzeta);
outputGradient[1] = static_cast<OT>(xi[1] * dxi + eta[1] * deta + zeta[1] * dzeta);
@ -75,24 +79,44 @@ struct StructuredPointGradient : public vtkm::worklet::WorkletPointNeighborhood
using CoordType = typename PointsIn::ValueType;
using OT = typename GradientOutType::ComponentType;
CoordType r = inputPoints.Portal.GetSpacing();
r[0] = (boundary.IsRadiusInXBoundary(1) ? r[0] * 0.5f : r[0]);
r[1] = (boundary.IsRadiusInYBoundary(1) ? r[1] * 0.5f : r[1]);
r[2] = (boundary.IsRadiusInZBoundary(1) ? r[2] * 0.5f : r[2]);
auto dx = inputField.Get(1, 0, 0) - inputField.Get(-1, 0, 0);
auto dy = inputField.Get(0, 1, 0) - inputField.Get(0, -1, 0);
auto dz = inputField.Get(0, 0, 1) - inputField.Get(0, 0, -1);
#if (defined(VTKM_CUDA) && defined(VTKM_GCC))
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wconversion"
#endif
outputGradient[0] = static_cast<OT>(dx * r[0]);
outputGradient[1] = static_cast<OT>(dy * r[1]);
outputGradient[2] = static_cast<OT>(dz * r[2]);
if (boundary.IsRadiusInXBoundary(1))
{
auto dx = inputField.GetUnchecked(1, 0, 0) - inputField.GetUnchecked(-1, 0, 0);
outputGradient[0] = static_cast<OT>(dx * (r[0] * 0.5f));
}
else
{
auto dx = inputField.Get(1, 0, 0) - inputField.Get(-1, 0, 0);
outputGradient[0] = static_cast<OT>(dx * r[0]);
}
if (boundary.IsRadiusInYBoundary(1))
{
auto dy = inputField.GetUnchecked(0, 1, 0) - inputField.GetUnchecked(0, -1, 0);
outputGradient[1] = static_cast<OT>(dy * r[1] * 0.5f);
}
else
{
auto dy = inputField.Get(0, 1, 0) - inputField.Get(0, -1, 0);
outputGradient[1] = static_cast<OT>(dy * (r[1]));
}
if (boundary.IsRadiusInZBoundary(1))
{
auto dz = inputField.GetUnchecked(0, 0, 1) - inputField.GetUnchecked(0, 0, -1);
outputGradient[2] = static_cast<OT>(dz * r[2] * 0.5f);
}
else
{
auto dz = inputField.Get(0, 0, 1) - inputField.Get(0, 0, -1);
outputGradient[2] = static_cast<OT>(dz * (r[2]));
}
#if (defined(VTKM_CUDA) && defined(VTKM_GCC))
#pragma GCC diagnostic pop
#endif
@ -103,20 +127,41 @@ struct StructuredPointGradient : public vtkm::worklet::WorkletPointNeighborhood
//will be float,3 even when T is a 3 component field
template <typename PointsIn, typename CT>
VTKM_EXEC void Jacobian(const PointsIn& inputPoints,
const vtkm::exec::BoundaryState& boundary,
const vtkm::Vec<bool, 3>& onBoundary,
vtkm::Vec<CT, 3>& m_xi,
vtkm::Vec<CT, 3>& m_eta,
vtkm::Vec<CT, 3>& m_zeta) const
{
using CoordType = typename PointsIn::ValueType;
CoordType xi, eta, zeta;
CoordType xi = inputPoints.Get(1, 0, 0) - inputPoints.Get(-1, 0, 0);
CoordType eta = inputPoints.Get(0, 1, 0) - inputPoints.Get(0, -1, 0);
CoordType zeta = inputPoints.Get(0, 0, 1) - inputPoints.Get(0, 0, -1);
xi = (boundary.IsRadiusInXBoundary(1) ? xi * 0.5f : xi);
eta = (boundary.IsRadiusInYBoundary(1) ? eta * 0.5f : eta);
zeta = (boundary.IsRadiusInZBoundary(1) ? zeta * 0.5f : zeta);
if (onBoundary[0])
{
xi = (inputPoints.Get(1, 0, 0) - inputPoints.Get(-1, 0, 0)) * 0.5f;
}
else
{
xi = inputPoints.GetUnchecked(1, 0, 0) - inputPoints.GetUnchecked(-1, 0, 0);
}
if (onBoundary[1])
{
eta = (inputPoints.Get(0, 1, 0) - inputPoints.Get(0, -1, 0)) * 0.5f;
}
else
{
eta = inputPoints.GetUnchecked(0, 1, 0) - inputPoints.GetUnchecked(0, -1, 0);
}
if (onBoundary[2])
{
zeta = (inputPoints.Get(0, 0, 1) - inputPoints.Get(0, 0, -1)) * 0.5f;
}
else
{
zeta = inputPoints.GetUnchecked(0, 0, 1) - inputPoints.GetUnchecked(0, 0, -1);
}
CT aj = xi[0] * eta[1] * zeta[2] + xi[1] * eta[2] * zeta[0] + xi[2] * eta[0] * zeta[1] -
xi[2] * eta[1] * zeta[0] - xi[1] * eta[0] * zeta[2] - xi[0] * eta[2] * zeta[1];