External Faces Uniform and Rectilear grids

Added new Run method to External Faces worklet to support
uniform and rectilinear grids, taking advantage of the grids
regular strucuture to speed-up computation.

Changed External Faces filter to call the proper Run method
in the External Faces worklet based on cell set type.

Added tests for uniform and rectilinear grids to External Faces
filter unit test.
This commit is contained in:
Thomas J. Otahal 2017-07-19 13:35:31 -06:00
parent 8c64247707
commit ab25c160ed
3 changed files with 430 additions and 1 deletions

@ -62,7 +62,14 @@ inline VTKM_CONT vtkm::filter::ResultDataSet ExternalFaces::DoExecute(
// external faces worklet
vtkm::cont::CellSetExplicit<> outCellSet(cells.GetName());
vtkm::worklet::ExternalFaces exfaces;
exfaces.Run(vtkm::filter::ApplyPolicyUnstructured(cells, policy), outCellSet, DeviceAdapter());
if (cells.IsSameType(vtkm::cont::CellSetStructured<3>()))
exfaces.Run(cells.Cast<vtkm::cont::CellSetStructured<3>>(),
input.GetCoordinateSystem(this->GetActiveCoordinateSystemIndex()),
outCellSet,
DeviceAdapter());
else
exfaces.Run(vtkm::filter::ApplyPolicyUnstructured(cells, policy), outCellSet, DeviceAdapter());
//3. create the output dataset
vtkm::cont::DataSet output;

@ -49,6 +49,16 @@ vtkm::cont::DataSet MakeDataTestSet2()
return MakeTestDataSet().Make3DExplicitDataSet5();
}
vtkm::cont::DataSet MakeDataTestSet3()
{
return MakeTestDataSet().Make3DUniformDataSet1();
}
vtkm::cont::DataSet MakeDataTestSet4()
{
return MakeTestDataSet().Make3DRectilinearDataSet0();
}
void TestExternalFacesExplicitGrid(const vtkm::cont::DataSet& ds,
bool compactPoints,
vtkm::Id numExpectedExtFaces,
@ -107,10 +117,32 @@ void TestWithHeterogeneousMesh()
TestExternalFacesExplicitGrid(ds, true, 12, 11);
}
void TestWithUniformMesh()
{
std::cout << "Testing with Uniform mesh\n";
vtkm::cont::DataSet ds = MakeDataTestSet3();
std::cout << "Compact Points Off\n";
TestExternalFacesExplicitGrid(ds, false, 16 * 6);
std::cout << "Compact Points On\n";
TestExternalFacesExplicitGrid(ds, true, 16 * 6, 98);
}
void TestWithRectilinearMesh()
{
std::cout << "Testing with Rectilinear mesh\n";
vtkm::cont::DataSet ds = MakeDataTestSet4();
std::cout << "Compact Points Off\n";
TestExternalFacesExplicitGrid(ds, false, 16);
std::cout << "Compact Points On\n";
TestExternalFacesExplicitGrid(ds, true, 16, 18);
}
void TestExternalFacesFilter()
{
TestWithHeterogeneousMesh();
TestWithHexahedraMesh();
TestWithUniformMesh();
TestWithRectilinearMesh();
}
} // anonymous namespace

@ -27,6 +27,7 @@
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/ArrayHandleConstant.h>
#include <vtkm/cont/ArrayHandleGroupVec.h>
#include <vtkm/cont/ArrayHandleGroupVecVariable.h>
#include <vtkm/cont/ArrayHandleIndex.h>
#include <vtkm/cont/ArrayHandlePermutation.h>
@ -93,6 +94,266 @@ struct ExternalFaces
}
};
//Worklet that returns the number of external faces for each structured cell
class NumExternalFacesPerStructuredCell : public vtkm::worklet::WorkletMapPointToCell
{
public:
typedef void ControlSignature(CellSetIn inCellSet,
FieldOut<> numFacesInCell,
FieldInPoint<Vec3> pointCoordinates);
typedef _2 ExecutionSignature(CellShape, _3);
typedef _1 InputDomain;
VTKM_CONT
NumExternalFacesPerStructuredCell(const vtkm::Vec<vtkm::Float32, 3>& min_point,
const vtkm::Vec<vtkm::Float32, 3>& max_point)
: MinPoint(min_point)
, MaxPoint(max_point)
{
}
VTKM_EXEC
inline vtkm::IdComponent CountExternalFacesOnDimension(vtkm::Float32 grid_min,
vtkm::Float32 grid_max,
vtkm::Float32 cell_min,
vtkm::Float32 cell_max) const
{
vtkm::IdComponent count = 0;
bool cell_min_at_grid_boundary = cell_min <= grid_min;
bool cell_max_at_grid_boundary = cell_max >= grid_max;
if (cell_min_at_grid_boundary and !cell_max_at_grid_boundary)
count++;
else if (!cell_min_at_grid_boundary and cell_max_at_grid_boundary)
count++;
else if (cell_min_at_grid_boundary and cell_max_at_grid_boundary)
count += 2;
return count;
}
template <typename CellShapeTag, typename PointCoordVecType>
VTKM_EXEC vtkm::IdComponent operator()(CellShapeTag shape,
const PointCoordVecType& pointCoordinates) const
{
VTKM_ASSERT(shape.Id == CELL_SHAPE_HEXAHEDRON);
vtkm::IdComponent count = 0;
count += CountExternalFacesOnDimension(
MinPoint[0], MaxPoint[0], pointCoordinates[0][0], pointCoordinates[1][0]);
count += CountExternalFacesOnDimension(
MinPoint[1], MaxPoint[1], pointCoordinates[0][1], pointCoordinates[3][1]);
count += CountExternalFacesOnDimension(
MinPoint[2], MaxPoint[2], pointCoordinates[0][2], pointCoordinates[4][2]);
return count;
}
private:
vtkm::Vec<vtkm::Float32, 3> MinPoint;
vtkm::Vec<vtkm::Float32, 3> MaxPoint;
};
//Worklet that finds face connectivity for each structured cell
class BuildConnectivityStructured : public vtkm::worklet::WorkletMapPointToCell
{
public:
typedef void ControlSignature(CellSetIn inCellSet,
WholeCellSetIn<> inputCell,
FieldOut<> faceShapes,
FieldOut<> facePointCount,
FieldOut<> faceConnectivity,
FieldInPoint<Vec3> pointCoordinates);
typedef void ExecutionSignature(CellShape, VisitIndex, InputIndex, _2, _3, _4, _5, _6);
typedef _1 InputDomain;
using ScatterType = vtkm::worklet::ScatterCounting;
VTKM_CONT
ScatterType GetScatter() const { return this->Scatter; }
template <typename CountArrayType, typename Device>
VTKM_CONT BuildConnectivityStructured(const vtkm::Vec<vtkm::Float32, 3>& min_point,
const vtkm::Vec<vtkm::Float32, 3>& max_point,
const CountArrayType& countArray,
Device)
: MinPoint(min_point)
, MaxPoint(max_point)
, Scatter(countArray, Device())
{
VTKM_IS_ARRAY_HANDLE(CountArrayType);
}
VTKM_CONT
BuildConnectivityStructured(const vtkm::Vec<vtkm::Float32, 3>& min_point,
const vtkm::Vec<vtkm::Float32, 3>& max_point,
const ScatterType& scatter)
: MinPoint(min_point)
, MaxPoint(max_point)
, Scatter(scatter)
{
}
enum FaceType
{
FACE_GRID_MIN,
FACE_GRID_MAX,
FACE_GRID_MIN_AND_MAX,
FACE_NONE
};
VTKM_EXEC
inline bool FoundFaceOnDimension(vtkm::Float32 grid_min,
vtkm::Float32 grid_max,
vtkm::Float32 cell_min,
vtkm::Float32 cell_max,
vtkm::Id& faceIndex,
vtkm::IdComponent& count,
vtkm::IdComponent dimensionFaceOffset,
vtkm::IdComponent visitIndex) const
{
bool cell_min_at_grid_boundary = cell_min <= grid_min;
bool cell_max_at_grid_boundary = cell_max >= grid_max;
FaceType Faces = FaceType::FACE_NONE;
if (cell_min_at_grid_boundary and !cell_max_at_grid_boundary)
Faces = FaceType::FACE_GRID_MIN;
else if (!cell_min_at_grid_boundary and cell_max_at_grid_boundary)
Faces = FaceType::FACE_GRID_MAX;
else if (cell_min_at_grid_boundary and cell_max_at_grid_boundary)
Faces = FaceType::FACE_GRID_MIN_AND_MAX;
if (Faces == FaceType::FACE_NONE)
return false;
if (Faces == FaceType::FACE_GRID_MIN)
{
if (visitIndex == count)
{
faceIndex = dimensionFaceOffset;
return true;
}
else
count++;
}
else if (Faces == FaceType::FACE_GRID_MAX)
{
if (visitIndex == count)
{
faceIndex = dimensionFaceOffset + 1;
return true;
}
else
count++;
}
else if (Faces == FaceType::FACE_GRID_MIN_AND_MAX)
{
if (visitIndex == count)
{
faceIndex = dimensionFaceOffset;
return true;
}
count++;
if (visitIndex == count)
{
faceIndex = dimensionFaceOffset + 1;
return true;
}
count++;
}
return false;
}
template <typename PointCoordVecType>
VTKM_EXEC inline vtkm::IdComponent FindFaceIndexForVisit(
vtkm::IdComponent visitIndex,
const PointCoordVecType& pointCoordinates) const
{
vtkm::IdComponent count = 0;
vtkm::Id faceIndex = 0;
// Search X dimension
if (!FoundFaceOnDimension(MinPoint[0],
MaxPoint[0],
pointCoordinates[0][0],
pointCoordinates[1][0],
faceIndex,
count,
0,
visitIndex))
{
// Search Y dimension
if (!FoundFaceOnDimension(MinPoint[1],
MaxPoint[1],
pointCoordinates[0][1],
pointCoordinates[3][1],
faceIndex,
count,
2,
visitIndex))
{
// Search Z dimension
FoundFaceOnDimension(MinPoint[2],
MaxPoint[2],
pointCoordinates[0][2],
pointCoordinates[4][2],
faceIndex,
count,
4,
visitIndex);
}
}
return faceIndex;
}
template <typename CellShapeTag,
typename CellSetType,
typename PointCoordVecType,
typename ConnectivityType>
VTKM_EXEC void operator()(CellShapeTag shape,
vtkm::IdComponent visitIndex,
vtkm::Id inputIndex,
const CellSetType& cellSet,
vtkm::UInt8& shapeOut,
vtkm::IdComponent& numFacePointsOut,
ConnectivityType& faceConnectivity,
const PointCoordVecType& pointCoordinates) const
{
VTKM_ASSERT(shape.Id == CELL_SHAPE_HEXAHEDRON);
vtkm::Id faceIndex = FindFaceIndexForVisit(visitIndex, pointCoordinates);
vtkm::VecCConst<vtkm::IdComponent> localFaceIndices =
vtkm::exec::CellFaceLocalIndices(faceIndex, shape, *this);
vtkm::IdComponent numFacePoints = localFaceIndices.GetNumberOfComponents();
VTKM_ASSERT(numFacePoints == connectivityOut.GetNumberOfComponents());
typename CellSetType::IndicesType inCellIndices = cellSet.GetIndices(inputIndex);
shapeOut = vtkm::CELL_SHAPE_QUAD;
numFacePointsOut = 4;
for (vtkm::IdComponent facePointIndex = 0; facePointIndex < numFacePoints; facePointIndex++)
{
faceConnectivity[facePointIndex] = inCellIndices[localFaceIndices[facePointIndex]];
}
}
private:
vtkm::Vec<vtkm::Float32, 3> MinPoint;
vtkm::Vec<vtkm::Float32, 3> MaxPoint;
ScatterType Scatter;
};
//Worklet that returns the number of faces for each cell/shape
class NumFacesPerCell : public vtkm::worklet::WorkletMapPointToCell
{
@ -355,6 +616,135 @@ struct ExternalFaces
};
public:
template <typename ShapeStorage,
typename NumIndicesStorage,
typename ConnectivityStorage,
typename OffsetsStorage,
typename DeviceAdapter>
VTKM_CONT void Run(const vtkm::cont::CellSetStructured<3>& inCellSet,
const vtkm::cont::CoordinateSystem& coord,
vtkm::cont::CellSetExplicit<ShapeStorage,
NumIndicesStorage,
ConnectivityStorage,
OffsetsStorage>& outCellSet,
DeviceAdapter)
{
vtkm::Vec<vtkm::Float32, 3> MinPoint;
vtkm::Vec<vtkm::Float32, 3> MaxPoint;
vtkm::exec::ConnectivityStructured<vtkm::TopologyElementTagPoint,
vtkm::TopologyElementTagCell,
3>
Conn;
Conn = inCellSet.PrepareForInput(
DeviceAdapter(), vtkm::TopologyElementTagPoint(), vtkm::TopologyElementTagCell());
vtkm::Id3 PointDimensions = Conn.GetPointDimensions();
typedef vtkm::cont::ArrayHandle<vtkm::FloatDefault> DefaultHandle;
typedef vtkm::cont::ArrayHandleCartesianProduct<DefaultHandle, DefaultHandle, DefaultHandle>
CartesianArrayHandle;
if (coord.GetData().IsSameType(CartesianArrayHandle()))
{
CartesianArrayHandle vertices;
vertices = coord.GetData().Cast<CartesianArrayHandle>();
MinPoint[0] =
static_cast<vtkm::Float32>(vertices.GetPortalConstControl().GetFirstPortal().Get(0));
MinPoint[1] =
static_cast<vtkm::Float32>(vertices.GetPortalConstControl().GetSecondPortal().Get(0));
MinPoint[2] =
static_cast<vtkm::Float32>(vertices.GetPortalConstControl().GetThirdPortal().Get(0));
MaxPoint[0] = static_cast<vtkm::Float32>(
vertices.GetPortalConstControl().GetFirstPortal().Get(PointDimensions[0] - 1));
MaxPoint[1] = static_cast<vtkm::Float32>(
vertices.GetPortalConstControl().GetSecondPortal().Get(PointDimensions[1] - 1));
MaxPoint[2] = static_cast<vtkm::Float32>(
vertices.GetPortalConstControl().GetThirdPortal().Get(PointDimensions[2] - 1));
}
else
{
vtkm::cont::ArrayHandleUniformPointCoordinates vertices;
vertices = coord.GetData().Cast<vtkm::cont::ArrayHandleUniformPointCoordinates>();
typedef typename vtkm::cont::ArrayHandleUniformPointCoordinates UniformArrayHandle;
typedef
typename UniformArrayHandle::ExecutionTypes<DeviceAdapter>::PortalConst UniformConstPortal;
UniformConstPortal Coordinates = vertices.PrepareForInput(DeviceAdapter());
MinPoint = Coordinates.GetOrigin();
vtkm::Vec<vtkm::Float32, 3> spacing = Coordinates.GetSpacing();
vtkm::Vec<vtkm::Float32, 3> unitLength;
unitLength[0] = static_cast<vtkm::Float32>(PointDimensions[0] - 1);
unitLength[1] = static_cast<vtkm::Float32>(PointDimensions[1] - 1);
unitLength[2] = static_cast<vtkm::Float32>(PointDimensions[2] - 1);
MaxPoint = MinPoint + spacing * unitLength;
}
// Create a worklet to count the number of external faces on each cell
vtkm::cont::ArrayHandle<vtkm::IdComponent> numExternalFaces;
vtkm::worklet::DispatcherMapTopology<NumExternalFacesPerStructuredCell>
numExternalFacesDispatcher((NumExternalFacesPerStructuredCell(MinPoint, MaxPoint)));
#ifdef __VTKM_EXTERNAL_FACES_BENCHMARK
vtkm::cont::Timer<DeviceAdapter> timer;
#endif
numExternalFacesDispatcher.Invoke(inCellSet, numExternalFaces, coord.GetData());
#ifdef __VTKM_EXTERNAL_FACES_BENCHMARK
std::cout << "NumExternalFacesPerStructuredCell_Worklet," << timer.GetElapsedTime() << "\n";
#endif
#ifdef __VTKM_EXTERNAL_FACES_BENCHMARK
timer.Reset();
#endif
typedef typename vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapter> DeviceAlgorithms;
vtkm::Id numberOfExternalFaces = DeviceAlgorithms::Reduce(numExternalFaces, 0, vtkm::Sum());
#ifdef __VTKM_EXTERNAL_FACES_BENCHMARK
std::cout << "numberOfExternalFaces_Reduce," << timer.GetElapsedTime() << "\n";
#endif
#ifdef __VTKM_EXTERNAL_FACES_BENCHMARK
timer.Reset();
#endif
vtkm::worklet::ScatterCounting scatterCellToExternalFace(numExternalFaces, DeviceAdapter());
#ifdef __VTKM_EXTERNAL_FACES_BENCHMARK
std::cout << "numExternalFaces_ScatterCounting," << timer.GetElapsedTime() << "\n";
#endif
numExternalFaces.ReleaseResources();
vtkm::Id connectivitySize = 4 * numberOfExternalFaces;
vtkm::cont::ArrayHandle<vtkm::Id, ConnectivityStorage> faceConnectivity;
vtkm::cont::ArrayHandle<vtkm::UInt8, ShapeStorage> faceShapes;
vtkm::cont::ArrayHandle<vtkm::IdComponent, NumIndicesStorage> facePointCount;
// Must pre allocate because worklet invocation will not have enough
// information to.
faceConnectivity.Allocate(connectivitySize);
vtkm::worklet::DispatcherMapTopology<BuildConnectivityStructured>
buildConnectivityStructuredDispatcher(
(BuildConnectivityStructured(MinPoint, MaxPoint, scatterCellToExternalFace)));
#ifdef __VTKM_EXTERNAL_FACES_BENCHMARK
timer.Reset();
#endif
buildConnectivityStructuredDispatcher.Invoke(
inCellSet,
inCellSet,
faceShapes,
facePointCount,
vtkm::cont::make_ArrayHandleGroupVec<4>(faceConnectivity),
coord.GetData());
#ifdef __VTKM_EXTERNAL_FACES_BENCHMARK
std::cout << "BuildConnectivityStructured_Worklet," << timer.GetElapsedTime() << "\n";
#endif
outCellSet.Fill(inCellSet.GetNumberOfPoints(), faceShapes, facePointCount, faceConnectivity);
#ifdef __VTKM_EXTERNAL_FACES_BENCHMARK
std::cout << "Total External Faces = " << outCellSet.GetNumberOfCells() << std::endl;
#endif
}
///////////////////////////////////////////////////
/// \brief ExternalFaces: Extract Faces on outside of geometry
template <typename InCellSetType,