Shuenhoy/vtk-m2
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
737ea7e156
vtk-m2/vtkm/worklet/ExternalFaces.h
Kenneth Moreland edc4c85fd9 Move Scatter from Worklet to Dispatcher 
Previously, when a Worklet needed a scatter, the scatter object was
stored in the Worklet object. That was problematic because that means
the Scatter, which is a control object, was shoved into the execution
environment.

To prevent that, move the Scatter into the Dispatcher object. The
worklet still declares a ScatterType alias, but no longer has a
GetScatter method. Instead, the Dispatcher now takes a Scatter object in
its constructor. If using the default scatter (ScatterIdentity), the
default constructor is used. If using another type of Scatter that
requires data to set up its state, then the caller of the worklet needs
to provide that to the dispatcher. For convenience, worklets are
encouraged to have a MakeScatter method to help construct a proper
scatter object.
2018-04-27 00:43:51 -04:00

978 lines
37 KiB
C++
Raw Blame History

//============================================================================
// 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 2014 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
// Copyright 2014 UT-Battelle, LLC.
// Copyright 2014 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_worklet_ExternalFaces_h
#define vtk_m_worklet_ExternalFaces_h
#include <vtkm/CellShape.h>
#include <vtkm/Hash.h>
#include <vtkm/Math.h>
#include <vtkm/exec/CellFace.h>
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/ArrayHandleConcatenate.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>
#include <vtkm/cont/ArrayHandleTransform.h>
#include <vtkm/cont/CellSetExplicit.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/DeviceAdapterAlgorithm.h>
#include <vtkm/cont/Field.h>
#include <vtkm/cont/Timer.h>
#include <vtkm/worklet/DispatcherMapTopology.h>
#include <vtkm/worklet/DispatcherReduceByKey.h>
#include <vtkm/worklet/Keys.h>
#include <vtkm/worklet/ScatterCounting.h>
#include <vtkm/worklet/WorkletMapTopology.h>
#include <vtkm/worklet/WorkletReduceByKey.h>
namespace vtkm
{
namespace worklet
{
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);
using InputDomain = _1;
VTKM_CONT
NumExternalFacesPerStructuredCell(const vtkm::Vec<vtkm::Float64, 3>& min_point,
const vtkm::Vec<vtkm::Float64, 3>& max_point)
: MinPoint(min_point)
, MaxPoint(max_point)
{
}
VTKM_EXEC
inline vtkm::IdComponent CountExternalFacesOnDimension(vtkm::Float64 grid_min,
vtkm::Float64 grid_max,
vtkm::Float64 cell_min,
vtkm::Float64 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 && !cell_max_at_grid_boundary)
{
count++;
}
else if (!cell_min_at_grid_boundary && cell_max_at_grid_boundary)
{
count++;
}
else if (cell_min_at_grid_boundary && 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
{
(void)shape; // C4100 false positive workaround
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::Float64, 3> MinPoint;
vtkm::Vec<vtkm::Float64, 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);
using InputDomain = _1;
using ScatterType = vtkm::worklet::ScatterCounting;
template <typename CountArrayType, typename Device>
VTKM_CONT static ScatterType MakeScatter(const CountArrayType& countArray, Device)
{
VTKM_IS_ARRAY_HANDLE(CountArrayType);
return ScatterType(countArray, Device());
}
VTKM_CONT
BuildConnectivityStructured(const vtkm::Vec<vtkm::Float64, 3>& min_point,
const vtkm::Vec<vtkm::Float64, 3>& max_point)
: MinPoint(min_point)
, MaxPoint(max_point)
{
}
enum FaceType
{
FACE_GRID_MIN,
FACE_GRID_MAX,
FACE_GRID_MIN_AND_MAX,
FACE_NONE
};
VTKM_EXEC
inline bool FoundFaceOnDimension(vtkm::Float64 grid_min,
vtkm::Float64 grid_max,
vtkm::Float64 cell_min,
vtkm::Float64 cell_max,
vtkm::IdComponent& 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 && !cell_max_at_grid_boundary)
{
Faces = FaceType::FACE_GRID_MIN;
}
else if (!cell_min_at_grid_boundary && cell_max_at_grid_boundary)
{
Faces = FaceType::FACE_GRID_MAX;
}
else if (cell_min_at_grid_boundary && 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::IdComponent 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::IdComponent faceIndex = FindFaceIndexForVisit(visitIndex, pointCoordinates);
const vtkm::IdComponent numFacePoints =
vtkm::exec::CellFaceNumberOfPoints(faceIndex, shape, *this);
VTKM_ASSERT(numFacePoints == faceConnectivity.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[vtkm::exec::CellFaceLocalIndex(facePointIndex, faceIndex, shape, *this)];
}
}
private:
vtkm::Vec<vtkm::Float64, 3> MinPoint;
vtkm::Vec<vtkm::Float64, 3> MaxPoint;
};
//Worklet that returns the number of faces for each cell/shape
class NumFacesPerCell : public vtkm::worklet::WorkletMapPointToCell
{
public:
typedef void ControlSignature(CellSetIn inCellSet, FieldOut<> numFacesInCell);
typedef _2 ExecutionSignature(CellShape);
using InputDomain = _1;
template <typename CellShapeTag>
VTKM_EXEC vtkm::IdComponent operator()(CellShapeTag shape) const
{
return vtkm::exec::CellFaceNumberOfFaces(shape, *this);
}
};
//Worklet that identifies a cell face by a hash value. Not necessarily completely unique.
class FaceHash : public vtkm::worklet::WorkletMapPointToCell
{
public:
typedef void ControlSignature(CellSetIn cellset,
FieldOut<> faceHashes,
FieldOut<> originCells,
FieldOut<> originFaces);
typedef void ExecutionSignature(_2, _3, _4, CellShape, FromIndices, InputIndex, VisitIndex);
using InputDomain = _1;
using ScatterType = vtkm::worklet::ScatterCounting;
template <typename CellShapeTag, typename CellNodeVecType>
VTKM_EXEC void operator()(vtkm::HashType& faceHash,
vtkm::Id& cellIndex,
vtkm::IdComponent& faceIndex,
CellShapeTag shape,
const CellNodeVecType& cellNodeIds,
vtkm::Id inputIndex,
vtkm::IdComponent visitIndex) const
{
faceHash = vtkm::Hash(vtkm::exec::CellFaceCanonicalId(visitIndex, shape, cellNodeIds, *this));
cellIndex = inputIndex;
faceIndex = visitIndex;
}
};
// Worklet that identifies the number of cells written out per face.
// Because there can be collisions in the face ids, this instance might
// represent multiple faces, which have to be checked. The resulting
// number is the total number of external faces.
class FaceCounts : public vtkm::worklet::WorkletReduceByKey
{
public:
typedef void ControlSignature(KeysIn keys,
WholeCellSetIn<> inputCells,
ValuesIn<> originCells,
ValuesIn<> originFaces,
ReducedValuesOut<> numOutputCells);
typedef _5 ExecutionSignature(_2, _3, _4);
using InputDomain = _1;
template <typename CellSetType, typename OriginCellsType, typename OriginFacesType>
VTKM_EXEC vtkm::IdComponent operator()(const CellSetType& cellSet,
const OriginCellsType& originCells,
const OriginFacesType& originFaces) const
{
vtkm::IdComponent numCellsOnHash = originCells.GetNumberOfComponents();
VTKM_ASSERT(originFaces.GetNumberOfComponents() == numCellsOnHash);
// Start by assuming all faces are unique, then remove one for each
// face we find a duplicate for.
vtkm::IdComponent numExternalFaces = numCellsOnHash;
for (vtkm::IdComponent myIndex = 0;
myIndex < numCellsOnHash - 1; // Don't need to check last face
myIndex++)
{
vtkm::Id3 myFace =
vtkm::exec::CellFaceCanonicalId(originFaces[myIndex],
cellSet.GetCellShape(originCells[myIndex]),
cellSet.GetIndices(originCells[myIndex]),
*this);
for (vtkm::IdComponent otherIndex = myIndex + 1; otherIndex < numCellsOnHash; otherIndex++)
{
vtkm::Id3 otherFace =
vtkm::exec::CellFaceCanonicalId(originFaces[otherIndex],
cellSet.GetCellShape(originCells[otherIndex]),
cellSet.GetIndices(originCells[otherIndex]),
*this);
if (myFace == otherFace)
{
// Faces are the same. Must be internal. Remove 2, one for each face. We don't have to
// worry about otherFace matching anything else because a proper topology will have at
// most 2 cells sharing a face, so there should be no more matches.
numExternalFaces -= 2;
break;
}
}
}
return numExternalFaces;
}
};
private:
// Resolves duplicate hashes by finding a specified unique face for a given hash.
// Given a cell set (from a WholeCellSetIn) and the cell/face id pairs for each face
// associated with a given hash, returns the index of the cell/face provided of the
// visitIndex-th unique face. Basically, this method searches through all the cell/face
// pairs looking for unique sets and returns the one associated with visitIndex.
template <typename CellSetType, typename OriginCellsType, typename OriginFacesType>
VTKM_EXEC static vtkm::IdComponent FindUniqueFace(const CellSetType& cellSet,
const OriginCellsType& originCells,
const OriginFacesType& originFaces,
vtkm::IdComponent visitIndex,
const vtkm::exec::FunctorBase* self)
{
vtkm::IdComponent numCellsOnHash = originCells.GetNumberOfComponents();
VTKM_ASSERT(originFaces.GetNumberOfComponents() == numCellsOnHash);
// Find the visitIndex-th unique face.
vtkm::IdComponent numFound = 0;
vtkm::IdComponent myIndex = 0;
while (true)
{
VTKM_ASSERT(myIndex < numCellsOnHash);
vtkm::Id3 myFace = vtkm::exec::CellFaceCanonicalId(originFaces[myIndex],
cellSet.GetCellShape(originCells[myIndex]),
cellSet.GetIndices(originCells[myIndex]),
*self);
bool foundPair = false;
for (vtkm::IdComponent otherIndex = myIndex + 1; otherIndex < numCellsOnHash; otherIndex++)
{
vtkm::Id3 otherFace =
vtkm::exec::CellFaceCanonicalId(originFaces[otherIndex],
cellSet.GetCellShape(originCells[otherIndex]),
cellSet.GetIndices(originCells[otherIndex]),
*self);
if (myFace == otherFace)
{
// Faces are the same. Must be internal.
foundPair = true;
break;
}
}
if (!foundPair)
{
if (numFound == visitIndex)
{
break;
}
else
{
numFound++;
}
}
myIndex++;
}
return myIndex;
}
public:
// Worklet that returns the number of points for each outputted face.
// Have to manage the case where multiple faces have the same hash.
class NumPointsPerFace : public vtkm::worklet::WorkletReduceByKey
{
public:
typedef void ControlSignature(KeysIn keys,
WholeCellSetIn<> inputCells,
ValuesIn<> originCells,
ValuesIn<> originFaces,
ReducedValuesOut<> numPointsInFace);
typedef _5 ExecutionSignature(_2, _3, _4, VisitIndex);
using InputDomain = _1;
using ScatterType = vtkm::worklet::ScatterCounting;
template <typename CountArrayType, typename Device>
VTKM_CONT static ScatterType MakeScatter(const CountArrayType& countArray, Device)
{
VTKM_IS_ARRAY_HANDLE(CountArrayType);
return ScatterType(countArray, Device());
}
template <typename CellSetType, typename OriginCellsType, typename OriginFacesType>
VTKM_EXEC vtkm::IdComponent operator()(const CellSetType& cellSet,
const OriginCellsType& originCells,
const OriginFacesType& originFaces,
vtkm::IdComponent visitIndex) const
{
vtkm::IdComponent myIndex =
ExternalFaces::FindUniqueFace(cellSet, originCells, originFaces, visitIndex, this);
return vtkm::exec::CellFaceNumberOfPoints(
originFaces[myIndex], cellSet.GetCellShape(originCells[myIndex]), *this);
}
};
// Worklet that returns the shape and connectivity for each external face
class BuildConnectivity : public vtkm::worklet::WorkletReduceByKey
{
public:
typedef void ControlSignature(KeysIn keys,
WholeCellSetIn<> inputCells,
ValuesIn<> originCells,
ValuesIn<> originFaces,
ReducedValuesOut<> shapesOut,
ReducedValuesOut<> connectivityOut,
ReducedValuesOut<> cellIdMapOut);
typedef void ExecutionSignature(_2, _3, _4, VisitIndex, _5, _6, _7);
using InputDomain = _1;
using ScatterType = vtkm::worklet::ScatterCounting;
template <typename CellSetType,
typename OriginCellsType,
typename OriginFacesType,
typename ConnectivityType>
VTKM_EXEC void operator()(const CellSetType& cellSet,
const OriginCellsType& originCells,
const OriginFacesType& originFaces,
vtkm::IdComponent visitIndex,
vtkm::UInt8& shapeOut,
ConnectivityType& connectivityOut,
vtkm::Id& cellIdMapOut) const
{
const vtkm::IdComponent myIndex =
ExternalFaces::FindUniqueFace(cellSet, originCells, originFaces, visitIndex, this);
const vtkm::IdComponent myFace = originFaces[myIndex];
typename CellSetType::CellShapeTag shapeIn = cellSet.GetCellShape(originCells[myIndex]);
shapeOut = vtkm::exec::CellFaceShape(myFace, shapeIn, *this);
cellIdMapOut = originCells[myIndex];
const vtkm::IdComponent numFacePoints =
vtkm::exec::CellFaceNumberOfPoints(myFace, shapeIn, *this);
VTKM_ASSERT(numFacePoints == connectivityOut.GetNumberOfComponents());
typename CellSetType::IndicesType inCellIndices = cellSet.GetIndices(originCells[myIndex]);
for (vtkm::IdComponent facePointIndex = 0; facePointIndex < numFacePoints; facePointIndex++)
{
connectivityOut[facePointIndex] =
inCellIndices[vtkm::exec::CellFaceLocalIndex(facePointIndex, myFace, shapeIn, *this)];
}
}
};
class IsPolyDataCell : public vtkm::worklet::WorkletMapPointToCell
{
public:
typedef void ControlSignature(CellSetIn inCellSet, FieldOut<> isPolyDataCell);
typedef _2 ExecutionSignature(CellShape);
using InputDomain = _1;
template <typename CellShapeTag>
VTKM_EXEC vtkm::IdComponent operator()(CellShapeTag shape) const
{
return !vtkm::exec::CellFaceNumberOfFaces(shape, *this);
}
};
class CountPolyDataCellPoints : public vtkm::worklet::WorkletMapPointToCell
{
public:
using ScatterType = vtkm::worklet::ScatterCounting;
typedef void ControlSignature(CellSetIn inCellSet, FieldOut<> numPoints);
typedef _2 ExecutionSignature(PointCount);
using InputDomain = _1;
VTKM_EXEC vtkm::Id operator()(vtkm::Id count) const { return count; }
};
class PassPolyDataCells : public vtkm::worklet::WorkletMapPointToCell
{
public:
using ScatterType = vtkm::worklet::ScatterCounting;
typedef void ControlSignature(CellSetIn inputTopology,
FieldOut<> shapes,
FieldOut<> pointIndices,
FieldOut<> cellIdMapOut);
typedef void ExecutionSignature(CellShape, PointIndices, InputIndex, _2, _3, _4);
template <typename CellShape, typename InPointIndexType, typename OutPointIndexType>
VTKM_EXEC void operator()(const CellShape& inShape,
const InPointIndexType& inPoints,
vtkm::Id inputIndex,
vtkm::UInt8& outShape,
OutPointIndexType& outPoints,
vtkm::Id& cellIdMapOut) const
{
cellIdMapOut = inputIndex;
outShape = inShape.Id;
vtkm::IdComponent numPoints = inPoints.GetNumberOfComponents();
VTKM_ASSERT(numPoints == outPoints.GetNumberOfComponents());
for (vtkm::IdComponent pointIndex = 0; pointIndex < numPoints; pointIndex++)
{
outPoints[pointIndex] = inPoints[pointIndex];
}
}
};
template <typename T>
struct BiasFunctor
{
VTKM_EXEC_CONT
BiasFunctor(T bias = T(0))
: Bias(bias)
{
}
VTKM_EXEC_CONT
T operator()(T x) const { return x + this->Bias; }
T Bias;
};
public:
VTKM_CONT
ExternalFaces()
: PassPolyData(true)
{
}
VTKM_CONT
void SetPassPolyData(bool flag) { this->PassPolyData = flag; }
VTKM_CONT
bool GetPassPolyData() const { return this->PassPolyData; }
//----------------------------------------------------------------------------
template <typename ValueType, typename StorageType, typename DeviceAdapter>
vtkm::cont::ArrayHandle<ValueType> ProcessCellField(
const vtkm::cont::ArrayHandle<ValueType, StorageType>& in,
const DeviceAdapter&) const
{
using Algo = vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapter>;
// Use a temporary permutation array to simplify the mapping:
auto tmp = vtkm::cont::make_ArrayHandlePermutation(this->CellIdMap, in);
// Copy into an array with default storage:
vtkm::cont::ArrayHandle<ValueType> result;
Algo::Copy(tmp, result);
return result;
}
void ReleaseCellMapArrays() { this->CellIdMap.ReleaseResources(); }
///////////////////////////////////////////////////
/// \brief ExternalFaces: Extract Faces on outside of geometry for regular grids.
///
/// Faster Run() method for uniform and rectilinear grid types.
/// Uses grid extents to find cells on the boundaries of the grid.
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::Float64, 3> MinPoint;
vtkm::Vec<vtkm::Float64, 3> MaxPoint;
vtkm::Id3 PointDimensions = inCellSet.GetPointDimensions();
using DefaultHandle = vtkm::cont::ArrayHandle<vtkm::FloatDefault>;
using CartesianArrayHandle =
vtkm::cont::ArrayHandleCartesianProduct<DefaultHandle, DefaultHandle, DefaultHandle>;
auto coordData = coord.GetData();
if (coordData.IsType<CartesianArrayHandle>())
{
auto vertices = coordData.Cast<CartesianArrayHandle>();
MinPoint[0] =
static_cast<vtkm::Float64>(vertices.GetPortalConstControl().GetFirstPortal().Get(0));
MinPoint[1] =
static_cast<vtkm::Float64>(vertices.GetPortalConstControl().GetSecondPortal().Get(0));
MinPoint[2] =
static_cast<vtkm::Float64>(vertices.GetPortalConstControl().GetThirdPortal().Get(0));
MaxPoint[0] = static_cast<vtkm::Float64>(
vertices.GetPortalConstControl().GetFirstPortal().Get(PointDimensions[0] - 1));
MaxPoint[1] = static_cast<vtkm::Float64>(
vertices.GetPortalConstControl().GetSecondPortal().Get(PointDimensions[1] - 1));
MaxPoint[2] = static_cast<vtkm::Float64>(
vertices.GetPortalConstControl().GetThirdPortal().Get(PointDimensions[2] - 1));
}
else
{
auto vertices = coordData.Cast<vtkm::cont::ArrayHandleUniformPointCoordinates>();
auto Coordinates = vertices.GetPortalConstControl();
MinPoint = Coordinates.GetOrigin();
vtkm::Vec<vtkm::Float64, 3> spacing = Coordinates.GetSpacing();
vtkm::Vec<vtkm::Float64, 3> unitLength;
unitLength[0] = static_cast<vtkm::Float64>(PointDimensions[0] - 1);
unitLength[1] = static_cast<vtkm::Float64>(PointDimensions[1] - 1);
unitLength[2] = static_cast<vtkm::Float64>(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, DeviceAdapter>
numExternalFacesDispatcher((NumExternalFacesPerStructuredCell(MinPoint, MaxPoint)));
numExternalFacesDispatcher.Invoke(inCellSet, numExternalFaces, coordData);
using DeviceAlgorithms = typename vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapter>;
vtkm::Id numberOfExternalFaces = DeviceAlgorithms::Reduce(numExternalFaces, 0, vtkm::Sum());
auto scatterCellToExternalFace =
BuildConnectivityStructured::MakeScatter(numExternalFaces, DeviceAdapter());
// Maps output cells to input cells. Store this for cell field mapping.
this->CellIdMap = scatterCellToExternalFace.GetOutputToInputMap();
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, DeviceAdapter>
buildConnectivityStructuredDispatcher(BuildConnectivityStructured(MinPoint, MaxPoint),
scatterCellToExternalFace);
buildConnectivityStructuredDispatcher.Invoke(
inCellSet,
inCellSet,
faceShapes,
facePointCount,
vtkm::cont::make_ArrayHandleGroupVec<4>(faceConnectivity),
coordData);
outCellSet.Fill(inCellSet.GetNumberOfPoints(), faceShapes, facePointCount, faceConnectivity);
}
///////////////////////////////////////////////////
/// \brief ExternalFaces: Extract Faces on outside of geometry
template <typename InCellSetType,
typename ShapeStorage,
typename NumIndicesStorage,
typename ConnectivityStorage,
typename OffsetsStorage,
typename DeviceAdapter>
VTKM_CONT void Run(const InCellSetType& inCellSet,
vtkm::cont::CellSetExplicit<ShapeStorage,
NumIndicesStorage,
ConnectivityStorage,
OffsetsStorage>& outCellSet,
DeviceAdapter)
{
using PointCountArrayType = vtkm::cont::ArrayHandle<vtkm::IdComponent, NumIndicesStorage>;
using ShapeArrayType = vtkm::cont::ArrayHandle<vtkm::UInt8, ShapeStorage>;
using OffsetsArrayType = vtkm::cont::ArrayHandle<vtkm::Id, OffsetsStorage>;
using ConnectivityArrayType = vtkm::cont::ArrayHandle<vtkm::Id, ConnectivityStorage>;
//Create a worklet to map the number of faces to each cell
vtkm::cont::ArrayHandle<vtkm::IdComponent> facesPerCell;
vtkm::worklet::DispatcherMapTopology<NumFacesPerCell, DeviceAdapter> numFacesDispatcher;
numFacesDispatcher.Invoke(inCellSet, facesPerCell);
vtkm::worklet::ScatterCounting scatterCellToFace(facesPerCell, DeviceAdapter());
facesPerCell.ReleaseResources();
PointCountArrayType polyDataPointCount;
ShapeArrayType polyDataShapes;
OffsetsArrayType polyDataOffsets;
ConnectivityArrayType polyDataConnectivity;
vtkm::cont::ArrayHandle<vtkm::Id> polyDataCellIdMap;
vtkm::Id polyDataConnectivitySize = 0;
if (this->PassPolyData)
{
vtkm::cont::ArrayHandle<vtkm::IdComponent> isPolyDataCell;
vtkm::worklet::DispatcherMapTopology<IsPolyDataCell, DeviceAdapter> isPolyDataCellDispatcher;
isPolyDataCellDispatcher.Invoke(inCellSet, isPolyDataCell);
vtkm::worklet::ScatterCounting scatterPolyDataCells(isPolyDataCell, DeviceAdapter());
isPolyDataCell.ReleaseResources();
if (scatterPolyDataCells.GetOutputRange(inCellSet.GetNumberOfCells()) != 0)
{
vtkm::worklet::DispatcherMapTopology<CountPolyDataCellPoints, DeviceAdapter>
countPolyDataCellPointsDispatcher(scatterPolyDataCells);
countPolyDataCellPointsDispatcher.Invoke(inCellSet, polyDataPointCount);
vtkm::cont::ConvertNumComponentsToOffsets(
polyDataPointCount, polyDataOffsets, polyDataConnectivitySize);
vtkm::worklet::DispatcherMapTopology<PassPolyDataCells, DeviceAdapter>
passPolyDataCellsDispatcher(scatterPolyDataCells);
polyDataConnectivity.Allocate(polyDataConnectivitySize);
passPolyDataCellsDispatcher.Invoke(
inCellSet,
polyDataShapes,
vtkm::cont::make_ArrayHandleGroupVecVariable(polyDataConnectivity, polyDataOffsets),
polyDataCellIdMap);
}
}
if (scatterCellToFace.GetOutputRange(inCellSet.GetNumberOfCells()) == 0)
{
if (!polyDataConnectivitySize)
{
// Data has no faces. Output is empty.
outCellSet.PrepareToAddCells(0, 0);
outCellSet.CompleteAddingCells(inCellSet.GetNumberOfPoints());
return;
}
else
{
// Pass only input poly data to output
outCellSet.Fill(inCellSet.GetNumberOfPoints(),
polyDataShapes,
polyDataPointCount,
polyDataConnectivity,
polyDataOffsets);
this->CellIdMap = polyDataCellIdMap;
return;
}
}
vtkm::cont::ArrayHandle<vtkm::HashType> faceHashes;
vtkm::cont::ArrayHandle<vtkm::Id> originCells;
vtkm::cont::ArrayHandle<vtkm::IdComponent> originFaces;
vtkm::worklet::DispatcherMapTopology<FaceHash, DeviceAdapter> faceHashDispatcher(
scatterCellToFace);
faceHashDispatcher.Invoke(inCellSet, faceHashes, originCells, originFaces);
vtkm::worklet::Keys<vtkm::HashType> faceKeys(faceHashes, DeviceAdapter());
vtkm::cont::ArrayHandle<vtkm::IdComponent> faceOutputCount;
vtkm::worklet::DispatcherReduceByKey<FaceCounts, DeviceAdapter> faceCountDispatcher;
faceCountDispatcher.Invoke(faceKeys, inCellSet, originCells, originFaces, faceOutputCount);
auto scatterCullInternalFaces = NumPointsPerFace::MakeScatter(faceOutputCount, DeviceAdapter());
PointCountArrayType facePointCount;
vtkm::worklet::DispatcherReduceByKey<NumPointsPerFace, DeviceAdapter> pointsPerFaceDispatcher(
scatterCullInternalFaces);
pointsPerFaceDispatcher.Invoke(faceKeys, inCellSet, originCells, originFaces, facePointCount);
ShapeArrayType faceShapes;
OffsetsArrayType faceOffsets;
vtkm::Id connectivitySize;
vtkm::cont::ConvertNumComponentsToOffsets(facePointCount, faceOffsets, connectivitySize);
ConnectivityArrayType faceConnectivity;
// Must pre allocate because worklet invocation will not have enough
// information to.
faceConnectivity.Allocate(connectivitySize);
vtkm::worklet::DispatcherReduceByKey<BuildConnectivity, DeviceAdapter>
buildConnectivityDispatcher(scatterCullInternalFaces);
vtkm::cont::ArrayHandle<vtkm::Id> faceToCellIdMap;
buildConnectivityDispatcher.Invoke(
faceKeys,
inCellSet,
originCells,
originFaces,
faceShapes,
vtkm::cont::make_ArrayHandleGroupVecVariable(faceConnectivity, faceOffsets),
faceToCellIdMap);
if (!polyDataConnectivitySize)
{
outCellSet.Fill(
inCellSet.GetNumberOfPoints(), faceShapes, facePointCount, faceConnectivity, faceOffsets);
this->CellIdMap = faceToCellIdMap;
}
else
{
// Join poly data to face data output
using DeviceAlgorithm = typename vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapter>;
vtkm::cont::ArrayHandleConcatenate<ShapeArrayType, ShapeArrayType> faceShapesArray(
faceShapes, polyDataShapes);
ShapeArrayType joinedShapesArray;
DeviceAlgorithm::Copy(faceShapesArray, joinedShapesArray);
vtkm::cont::ArrayHandleConcatenate<PointCountArrayType, PointCountArrayType> pointCountArray(
facePointCount, polyDataPointCount);
PointCountArrayType joinedPointCountArray;
DeviceAlgorithm::Copy(pointCountArray, joinedPointCountArray);
vtkm::cont::ArrayHandleConcatenate<ConnectivityArrayType, ConnectivityArrayType>
connectivityArray(faceConnectivity, polyDataConnectivity);
ConnectivityArrayType joinedConnectivity;
DeviceAlgorithm::Copy(connectivityArray, joinedConnectivity);
// Adjust poly data offsets array with face connectivity size before join
using TransformBiasArrayType =
vtkm::cont::ArrayHandleTransform<OffsetsArrayType, BiasFunctor<vtkm::Id>>;
TransformBiasArrayType adjustedPolyDataOffsets =
vtkm::cont::make_ArrayHandleTransform<OffsetsArrayType>(
polyDataOffsets, BiasFunctor<vtkm::Id>(faceConnectivity.GetNumberOfValues()));
vtkm::cont::ArrayHandleConcatenate<OffsetsArrayType, TransformBiasArrayType> offsetsArray(
faceOffsets, adjustedPolyDataOffsets);
OffsetsArrayType joinedOffsets;
DeviceAlgorithm::Copy(offsetsArray, joinedOffsets);
vtkm::cont::ArrayHandleConcatenate<vtkm::cont::ArrayHandle<vtkm::Id>,
vtkm::cont::ArrayHandle<vtkm::Id>>
cellIdMapArray(faceToCellIdMap, polyDataCellIdMap);
vtkm::cont::ArrayHandle<vtkm::Id> joinedCellIdMap;
DeviceAlgorithm::Copy(cellIdMapArray, joinedCellIdMap);
outCellSet.Fill(inCellSet.GetNumberOfPoints(),
joinedShapesArray,
joinedPointCountArray,
joinedConnectivity,
joinedOffsets);
this->CellIdMap = joinedCellIdMap;
}
}
private:
vtkm::cont::ArrayHandle<vtkm::Id> CellIdMap;
bool PassPolyData;
}; //struct ExternalFaces
}
} //namespace vtkm::worklet
#endif //vtk_m_worklet_ExternalFaces_h
Reference in New Issue View Git Blame Copy Permalink