marching tetrahedra generates triangles

This commit is contained in:
Li-Ta Lo 2019-08-15 18:03:53 -06:00
parent 403bca94c7
commit d82d380979
5 changed files with 444 additions and 107 deletions

@ -127,6 +127,16 @@ public:
this->operator()(
vtkm::CellShapeTagHexahedron(), isovalues, fieldIn, numTriangles, numTrianglesTable);
}
else if (shape.Id == CELL_SHAPE_QUAD)
{
this->operator()(
vtkm::CellShapeTagQuad(), isovalues, fieldIn, numTriangles, numTrianglesTable);
}
else if (shape.Id == CELL_SHAPE_TETRA)
{
this->operator()(
vtkm::CellShapeTagTetra(), isovalues, fieldIn, numTriangles, numTrianglesTable);
}
else
{
numTriangles = 0;
@ -135,11 +145,44 @@ public:
template <typename IsoValuesType, typename FieldInType, typename NumTrianglesTablePortalType>
VTKM_EXEC void operator()(vtkm::CellShapeTagQuad vtkmNotUsed(shape),
const IsoValuesType& vtkmNotUsed(isovalues),
const FieldInType& vtkmNotUsed(fieldIn),
vtkm::IdComponent& vtkmNotUsed(numTriangles),
const IsoValuesType& isovalues,
const FieldInType& fieldIn,
vtkm::IdComponent& numTriangles,
const NumTrianglesTablePortalType& vtkmNotUsed(numTrianglesTable)) const
{
static vtkm::IdComponent numLinesTable[16] = { 0, 1, 1, 1, 1, 2, 1, 1, 1, 1, 2, 1, 1, 1, 1, 0 };
vtkm::IdComponent sum = 0;
for (vtkm::Id i = 0; i < isovalues.GetNumberOfValues(); ++i)
{
const vtkm::IdComponent caseNumber =
((fieldIn[0] > isovalues[i]) | (fieldIn[1] > isovalues[i]) << 1 |
(fieldIn[2] > isovalues[i]) << 2 | (fieldIn[3] > isovalues[i]) << 3);
sum += numLinesTable[caseNumber];
}
numTriangles = sum;
std::cout << "num lines: " << sum << std::endl;
}
template <typename IsoValuesType, typename FieldInType, typename NumTrianglesTablePortalType>
VTKM_EXEC void operator()(vtkm::CellShapeTagTetra vtkmNotUsed(shape),
const IsoValuesType& isovalues,
const FieldInType& fieldIn,
vtkm::IdComponent& numTriangles,
const NumTrianglesTablePortalType& vtkmNotUsed(numTrianglesTable)) const
{
static vtkm::IdComponent numTrianglesTable[16] = { 0, 1, 1, 2, 1, 2, 2, 1,
1, 2, 2, 1, 2, 1, 1, 0 };
vtkm::IdComponent caseNumber;
vtkm::IdComponent sum = 0;
for (vtkm::Id i = 0; i < isovalues.GetNumberOfValues(); ++i)
{
caseNumber = ((fieldIn[0] > isovalues[i]) | (fieldIn[1] > isovalues[i]) << 1 |
(fieldIn[2] > isovalues[i]) << 2 | (fieldIn[3] > isovalues[i]) << 3);
sum += numTrianglesTable[caseNumber];
}
numTriangles = sum;
std::cout << "tetra, case: " << caseNumber << ", num triangles: " << numTriangles << std::endl;
}
template <typename IsoValuesType, typename FieldInType, typename NumTrianglesTablePortalType>
@ -321,6 +364,28 @@ public:
visitIndex,
indices);
}
else if (shape.Id == CELL_SHAPE_QUAD)
{
this->operator()(vtkm::CellShapeTagQuad(),
isovalues,
fieldIn,
metaData,
inputCellId,
outputCellId,
visitIndex,
indices);
}
else if (shape.Id == CELL_SHAPE_TETRA)
{
this->operator()(vtkm::CellShapeTagTetra(),
isovalues,
fieldIn,
metaData,
inputCellId,
outputCellId,
visitIndex,
indices);
}
}
template <typename IsoValuesType,
@ -339,6 +404,96 @@ public:
{ //covers when we have quads coming from 2d structured data
}
template <typename IsoValuesType,
typename FieldInType, // Vec-like, one per input point
typename IndicesVecType,
typename DeviceAdapter>
VTKM_EXEC void operator()(CellShapeTagTetra vtkmNotUsed(shape),
const IsoValuesType& isovalues,
const FieldInType& fieldIn, // Input point field defining the contour
EdgeWeightGenerateMetaData::ExecObject<DeviceAdapter>& metaData,
vtkm::Id inputCellId,
vtkm::Id outputCellId,
vtkm::IdComponent visitIndex,
const IndicesVecType& indices) const
{ //covers when we have tetra coming from 3d structured data
const vtkm::Id outputPointId = 3 * outputCellId;
using FieldType = typename vtkm::VecTraits<FieldInType>::ComponentType;
static vtkm::IdComponent numTrianglesTable[16] = { 0, 1, 1, 2, 1, 2, 2, 1,
1, 2, 2, 1, 2, 1, 1, 0 };
// TODO: I don't fully understand the following logic about multiple isovalues
vtkm::IdComponent sum = 0, caseNumber = 0;
vtkm::IdComponent i = 0, size = static_cast<vtkm::IdComponent>(isovalues.GetNumberOfValues());
for (i = 0; i < size; ++i)
{
const FieldType ivalue = isovalues[i];
// Compute the Marching Cubes case number for this cell. We need to iterate
// the isovalues until the sum >= our visit index. But we need to make
// sure the caseNumber is correct before stopping
caseNumber = ((fieldIn[0] > ivalue) | (fieldIn[1] > ivalue) << 1 |
(fieldIn[2] > ivalue) << 2 | (fieldIn[3] > ivalue) << 3);
sum += numTrianglesTable[caseNumber];
if (sum > visitIndex)
{
break;
}
}
//visitIndex = sum - visitIndex - 1;
static const vtkm::IdComponent triTable[16][7] = {
{ -1, -1, -1, -1, -1, -1, -1 }, { 0, 3, 2, -1, -1, -1, -1 },
{ 0, 1, 4, -1, -1, -1, -1 }, { 1, 4, 2, 2, 4, 3, -1 },
{ 1, 2, 5, -1, -1, -1, -1 }, { 0, 3, 5, 0, 5, 1, -1 },
{ 0, 2, 5, 0, 5, 4, -1 }, { 5, 4, 3, -1, -1, -1, -1 },
{ 3, 4, 5, -1, -1, -1, -1 }, { 4, 5, 0, 5, 2, 0, -1 },
{ 1, 5, 0, 5, 3, 0, -1 }, { 5, 2, 1, -1, -1, -1, -1 },
{ 3, 4, 2, 2, 4, 1, -1 }, { 4, 1, 0, -1, -1, -1, -1 },
{ 2, 3, 0, -1, -1, -1, -1 }, { -1, -1, -1, -1, -1, -1, -1 },
};
const vtkm::IdComponent* triTableEntry = &triTable[caseNumber][visitIndex * 3];
std::cout << "caseNumber: " << caseNumber << ", visitIndex: " << visitIndex << std::endl;
const int EdgeTable[] = {
0, 1, // edge 0 : vertex 0 -> vertex 1
1, 2, // edge 1 : vertex 1 -> vertex 2
0, 2, // edge 2 : vertex 0 -> vertex 2
0, 3, // edge 3 : vertex 0 -> vertex 3
1, 3, // edge 4 : vertex 1 -> vertex 3
2, 3 // edge 5 : vertex 2 -> vertex 3
};
#if 1
// Interpolate for vertex positions and associated scalar values
for (vtkm::IdComponent triVertex = 0; triVertex < 3; triVertex++)
{
const vtkm::IdComponent edgeIndex = triTableEntry[triVertex];
const vtkm::IdComponent edgeVertex0 = EdgeTable[2 * edgeIndex + 0];
const vtkm::IdComponent edgeVertex1 = EdgeTable[2 * edgeIndex + 1];
const FieldType fieldValue0 = fieldIn[edgeVertex0];
const FieldType fieldValue1 = fieldIn[edgeVertex1];
// Store the input cell id so that we can properly generate the normals
// in a subsequent call, after we have merged duplicate points
metaData.InterpCellIdPortal.Set(outputPointId + triVertex, inputCellId);
metaData.InterpContourPortal.Set(outputPointId + triVertex, static_cast<vtkm::UInt8>(i));
metaData.InterpIdPortal.Set(outputPointId + triVertex,
vtkm::Id2(indices[edgeVertex0], indices[edgeVertex1]));
vtkm::FloatDefault interpolant = static_cast<vtkm::FloatDefault>(isovalues[i] - fieldValue0) /
static_cast<vtkm::FloatDefault>(fieldValue1 - fieldValue0);
metaData.InterpWeightsPortal.Set(outputPointId + triVertex, interpolant);
}
#endif
}
template <typename IsoValuesType,
typename FieldInType, // Vec-like, one per input point
typename IndicesVecType,
@ -354,6 +509,7 @@ public:
const IndicesVecType& vtkmNotUsed(indices)) const
{ //covers when we have quads coming from 2d structured data
}
template <typename IsoValuesType,
typename FieldInType, // Vec-like, one per input point
typename IndicesVecType,

@ -9,78 +9,80 @@
##============================================================================
set(unit_tests
UnitTestAverageByKey.cxx
UnitTestBoundingIntervalHierarchy.cxx
UnitTestCellAverage.cxx
UnitTestCellDeepCopy.cxx
UnitTestCellGradient.cxx
UnitTestCellSetConnectivity.cxx
UnitTestCellSetDualGraph.cxx
UnitTestCellMeasure.cxx
UnitTestClipping.cxx
UnitTestContourTreeUniform.cxx
UnitTestContourTreeUniformAugmented.cxx
UnitTestCoordinateSystemTransform.cxx
UnitTestCosmoTools.cxx
UnitTestCrossProduct.cxx
UnitTestDotProduct.cxx
UnitTestExternalFaces.cxx
UnitTestExtractGeometry.cxx
UnitTestExtractPoints.cxx
UnitTestExtractStructured.cxx
UnitTestFieldHistogram.cxx
UnitTestFieldStatistics.cxx
UnitTestGraphConnectivity.cxx
UnitTestInnerJoin.cxx
UnitTestImageConnectivity.cxx
UnitTestKdTreeBuildNNS.cxx
UnitTestKeys.cxx
UnitTestMagnitude.cxx
UnitTestMarchingCubes.cxx
UnitTestMask.cxx
UnitTestMaskIndices.cxx
UnitTestMaskPoints.cxx
UnitTestMaskSelect.cxx
UnitTestNormalize.cxx
UnitTestNDimsEntropy.cxx
UnitTestNDimsHistogram.cxx
UnitTestNDimsHistMarginalization.cxx
UnitTestParticleAdvection.cxx
UnitTestPointElevation.cxx
UnitTestPointGradient.cxx
UnitTestPointTransform.cxx
UnitTestProbe.cxx
UnitTestRemoveUnusedPoints.cxx
UnitTestScalarsToColors.cxx
UnitTestScatterAndMask.cxx
UnitTestScatterCounting.cxx
UnitTestScatterPermutation.cxx
UnitTestSplatKernels.cxx
UnitTestSplitSharpEdges.cxx
UnitTestStreamingSine.cxx
UnitTestStreamLineUniformGrid.cxx
UnitTestSurfaceNormals.cxx
UnitTestTemporalAdvection.cxx
UnitTestTetrahedralize.cxx
UnitTestThreshold.cxx
UnitTestThresholdPoints.cxx
UnitTestTriangulate.cxx
UnitTestTube.cxx
UnitTestWholeCellSetIn.cxx
UnitTestWorkletMapField.cxx
UnitTestWorkletMapFieldExecArg.cxx
UnitTestWorkletMapFieldWholeArray.cxx
UnitTestWorkletMapFieldWholeArrayAtomic.cxx
UnitTestWorkletMapPointNeighborhood.cxx
UnitTestWorkletMapTopologyExplicit.cxx
UnitTestWorkletMapTopologyUniform.cxx
UnitTestWorkletReduceByKey.cxx
UnitTestVertexClustering.cxx
UnitTestWarpScalar.cxx
UnitTestWarpVector.cxx
UnitTestWaveletCompressor.cxx
UnitTestWaveletGenerator.cxx
UnitTestZFPCompressor.cxx
# UnitTestAverageByKey.cxx
# UnitTestBoundingIntervalHierarchy.cxx
# UnitTestCellAverage.cxx
# UnitTestCellDeepCopy.cxx
# UnitTestCellGradient.cxx
# UnitTestCellSetConnectivity.cxx
# UnitTestCellSetDualGraph.cxx
# UnitTestCellMeasure.cxx
# UnitTestClipping.cxx
# UnitTestContourTreeUniform.cxx
# UnitTestContourTreeUniformAugmented.cxx
# UnitTestCoordinateSystemTransform.cxx
# UnitTestCosmoTools.cxx
# UnitTestCrossProduct.cxx
# UnitTestDotProduct.cxx
# UnitTestExternalFaces.cxx
# UnitTestExtractGeometry.cxx
# UnitTestExtractPoints.cxx
# UnitTestExtractStructured.cxx
# UnitTestFieldHistogram.cxx
# UnitTestFieldStatistics.cxx
# UnitTestGraphConnectivity.cxx
# UnitTestInnerJoin.cxx
# UnitTestImageConnectivity.cxx
# UnitTestKdTreeBuildNNS.cxx
# UnitTestKeys.cxx
# UnitTestMagnitude.cxx
# UnitTestMarchingCubes.cxx
# UnitTestMarchingSquares.cxx
UnitTestMarchingTetrahedra.cxx
# UnitTestMask.cxx
# UnitTestMaskIndices.cxx
# UnitTestMaskPoints.cxx
# UnitTestMaskSelect.cxx
# UnitTestNormalize.cxx
# UnitTestNDimsEntropy.cxx
# UnitTestNDimsHistogram.cxx
# UnitTestNDimsHistMarginalization.cxx
# UnitTestParticleAdvection.cxx
# UnitTestPointElevation.cxx
# UnitTestPointGradient.cxx
# UnitTestPointTransform.cxx
# UnitTestProbe.cxx
# UnitTestRemoveUnusedPoints.cxx
# UnitTestScalarsToColors.cxx
# UnitTestScatterAndMask.cxx
# UnitTestScatterCounting.cxx
# UnitTestScatterPermutation.cxx
# UnitTestSplatKernels.cxx
# UnitTestSplitSharpEdges.cxx
# UnitTestStreamingSine.cxx
# UnitTestStreamLineUniformGrid.cxx
# UnitTestSurfaceNormals.cxx
# UnitTestTemporalAdvection.cxx
# UnitTestTetrahedralize.cxx
# UnitTestThreshold.cxx
# UnitTestThresholdPoints.cxx
# UnitTestTriangulate.cxx
# UnitTestTube.cxx
# UnitTestWholeCellSetIn.cxx
# UnitTestWorkletMapField.cxx
# UnitTestWorkletMapFieldExecArg.cxx
# UnitTestWorkletMapFieldWholeArray.cxx
# UnitTestWorkletMapFieldWholeArrayAtomic.cxx
# UnitTestWorkletMapPointNeighborhood.cxx
# UnitTestWorkletMapTopologyExplicit.cxx
# UnitTestWorkletMapTopologyUniform.cxx
# UnitTestWorkletReduceByKey.cxx
# UnitTestVertexClustering.cxx
# UnitTestWarpScalar.cxx
# UnitTestWarpVector.cxx
# UnitTestWaveletCompressor.cxx
# UnitTestWaveletGenerator.cxx
# UnitTestZFPCompressor.cxx
)

@ -0,0 +1,86 @@
//============================================================================
// 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.
//============================================================================
#include <vtkm/VectorAnalysis.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/filter/MarchingCubes.h>
#include <vtkm/worklet/WorkletMapField.h>
class SineWave : public vtkm::worklet::WorkletMapField
{
public:
using ControlSignature = void(FieldIn coord, FieldOut v);
using ExecutionSignature = void(_1, _2);
using InputDomain = _1;
VTKM_EXEC
void operator()(const vtkm::Vec<vtkm::FloatDefault, 3> coord, vtkm::Float32& v) const
{
v = vtkm::MagnitudeSquared(coord);
}
};
vtkm::cont::DataSet MakeSineWaveDataSet(vtkm::Id3 dims)
{
vtkm::cont::DataSet dataset;
// create a point coordinate system
const vtkm::Id3 vdims{ dims[0] + 1, dims[1] + 1, dims[2] + 1 };
vtkm::cont::ArrayHandleUniformPointCoordinates coordinates{
vdims,
vtkm::Vec<vtkm::FloatDefault, 3>(-1.0, -1.0, 0),
vtkm::Vec<vtkm::FloatDefault, 3>(2.0f / vdims[0], 2.0f / vdims[1], 0.0)
};
dataset.AddCoordinateSystem(vtkm::cont::CoordinateSystem("coordinates", coordinates));
// generate point field from point coordinates.
vtkm::cont::ArrayHandle<vtkm::Float32> pointFieldArray;
vtkm::worklet::DispatcherMapField<SineWave> dispatcherMapField(SineWave{});
dispatcherMapField.Invoke(coordinates, pointFieldArray);
dataset.AddField(
vtkm::cont::Field("nodevar", vtkm::cont::Field::Association::POINTS, pointFieldArray));
// add cell set
vtkm::cont::CellSetStructured<2> cellSet("cells");
cellSet.SetPointDimensions({ vdims[0], vdims[1] });
dataset.AddCellSet(cellSet);
return dataset;
}
void TestMarchingSquares()
{
vtkm::cont::DataSet input = MakeSineWaveDataSet({ 2, 2, 0 });
vtkm::cont::ArrayHandle<vtkm::Float32> pointFieldArray;
input.GetField("nodevar").GetData().CopyTo(pointFieldArray);
vtkm::Float32 isoValue = 0.5;
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float32, 3>> verticesArray;
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float32, 3>> normalsArray;
vtkm::cont::ArrayHandle<vtkm::Float32> scalarsArray;
vtkm::worklet::MarchingCubes isosurfaceFilter;
isosurfaceFilter.SetMergeDuplicatePoints(false);
auto result = isosurfaceFilter.Run(&isoValue,
1,
input.GetCellSet(0),
input.GetCoordinateSystem(),
pointFieldArray,
verticesArray,
normalsArray);
}
int UnitTestMarchingSquares(int argc, char* argv[])
{
return vtkm::cont::testing::Testing::Run(TestMarchingSquares, argc, argv);
}

@ -0,0 +1,93 @@
//============================================================================
// 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.
//============================================================================
#include <vtkm/VectorAnalysis.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/filter/MarchingCubes.h>
#include <vtkm/filter/Tetrahedralize.h>
#include <vtkm/worklet/WorkletMapField.h>
class Tangle : public vtkm::worklet::WorkletMapField
{
public:
using ControlSignature = void(FieldIn coord, FieldOut v);
using ExecutionSignature = void(_1, _2);
using InputDomain = _1;
VTKM_EXEC
void operator()(const vtkm::Vec<vtkm::FloatDefault, 3> coord, vtkm::Float32& v) const
{
v = vtkm::MagnitudeSquared(coord);
}
};
vtkm::cont::DataSet MakeTangleDataSet(vtkm::Id3 dims)
{
vtkm::cont::DataSet dataset;
// create a point coordinate system
const vtkm::Id3 vdims{ dims[0] + 1, dims[1] + 1, dims[2] + 1 };
vtkm::cont::ArrayHandleUniformPointCoordinates coordinates{
vdims,
vtkm::Vec<vtkm::FloatDefault, 3>(-1.0, -1.0, -1.0),
vtkm::Vec<vtkm::FloatDefault, 3>(2.0f / vdims[0], 2.0f / vdims[1], 2.0f / vdims[2])
};
dataset.AddCoordinateSystem(vtkm::cont::CoordinateSystem("coordinates", coordinates));
// generate point field from point coordinates.
vtkm::cont::ArrayHandle<vtkm::Float32> pointFieldArray;
vtkm::worklet::DispatcherMapField<Tangle> dispatcherMapField(Tangle{});
dispatcherMapField.Invoke(coordinates, pointFieldArray);
dataset.AddField(
vtkm::cont::Field("nodevar", vtkm::cont::Field::Association::POINTS, pointFieldArray));
// add cell set
vtkm::cont::CellSetStructured<3> cellSet("cells");
cellSet.SetPointDimensions(vdims);
dataset.AddCellSet(cellSet);
return dataset;
}
void TestMarchingTetrahedra()
{
vtkm::cont::DataSet input = MakeTangleDataSet({ 2, 2, 2 });
vtkm::filter::Tetrahedralize tetrahedralize;
vtkm::cont::DataSet tetra = tetrahedralize.Execute(input);
std::cout << "number of tetras: " << tetra.GetCellSet().GetNumberOfCells() << std::endl;
vtkm::cont::ArrayHandle<vtkm::Float32> pointFieldArray;
tetra.GetField("nodevar").GetData().CopyTo(pointFieldArray);
vtkm::Float32 isoValue = 0.5;
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float32, 3>> verticesArray;
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float32, 3>> normalsArray;
vtkm::cont::ArrayHandle<vtkm::Float32> scalarsArray;
vtkm::worklet::MarchingCubes isosurfaceFilter;
isosurfaceFilter.SetMergeDuplicatePoints(false);
auto result = isosurfaceFilter.Run(&isoValue,
1,
tetra.GetCellSet(0),
tetra.GetCoordinateSystem(),
pointFieldArray,
verticesArray,
normalsArray);
result.PrintSummary(std::cout);
}
int UnitTestMarchingTetrahedra(int argc, char* argv[])
{
return vtkm::cont::testing::Testing::Run(TestMarchingTetrahedra, argc, argv);
}

@ -110,39 +110,39 @@ public:
return outCellSet;
}
vtkm::cont::CellSetSingleType<> Run(const vtkm::cont::CellSetExplicit<>& cellSet,
vtkm::cont::ArrayHandle<vtkm::IdComponent>& outCellsPerCell)
{
vtkm::cont::CellSetSingleType<> outCellSet(cellSet.GetName());
// Input topology
auto inShapes =
cellSet.GetShapesArray(vtkm::TopologyElementTagPoint(), vtkm::TopologyElementTagCell());
auto inNumIndices =
cellSet.GetNumIndicesArray(vtkm::TopologyElementTagPoint(), vtkm::TopologyElementTagCell());
// Output topology
vtkm::cont::ArrayHandle<vtkm::Id> outConnectivity;
vtkm::worklet::internal::TetrahedralizeTables tables;
// Determine the number of output cells each input cell will generate
vtkm::worklet::DispatcherMapField<TetrahedraPerCell> tetPerCellDispatcher;
tetPerCellDispatcher.Invoke(inShapes, tables.PrepareForInput(), outCellsPerCell);
// Build new cells
vtkm::worklet::DispatcherMapTopology<TetrahedralizeCell> tetrahedralizeDispatcher(
TetrahedralizeCell::MakeScatter(outCellsPerCell));
tetrahedralizeDispatcher.Invoke(
cellSet, tables.PrepareForInput(), vtkm::cont::make_ArrayHandleGroupVec<4>(outConnectivity));
// Add cells to output cellset
outCellSet.Fill(cellSet.GetNumberOfPoints(), vtkm::CellShapeTagTetra::Id, 4, outConnectivity);
return outCellSet;
}
};
template <>
vtkm::cont::CellSetSingleType<> TetrahedralizeExplicit::Run(
const vtkm::cont::CellSetExplicit<>& cellSet,
vtkm::cont::ArrayHandle<vtkm::IdComponent>& outCellsPerCell)
{
vtkm::cont::CellSetSingleType<> outCellSet(cellSet.GetName());
// Input topology
auto inShapes =
cellSet.GetShapesArray(vtkm::TopologyElementTagPoint(), vtkm::TopologyElementTagCell());
auto inNumIndices =
cellSet.GetNumIndicesArray(vtkm::TopologyElementTagPoint(), vtkm::TopologyElementTagCell());
// Output topology
vtkm::cont::ArrayHandle<vtkm::Id> outConnectivity;
vtkm::worklet::internal::TetrahedralizeTables tables;
// Determine the number of output cells each input cell will generate
vtkm::worklet::DispatcherMapField<TetrahedraPerCell> tetPerCellDispatcher;
tetPerCellDispatcher.Invoke(inShapes, tables.PrepareForInput(), outCellsPerCell);
// Build new cells
vtkm::worklet::DispatcherMapTopology<TetrahedralizeCell> tetrahedralizeDispatcher(
TetrahedralizeCell::MakeScatter(outCellsPerCell));
tetrahedralizeDispatcher.Invoke(
cellSet, tables.PrepareForInput(), vtkm::cont::make_ArrayHandleGroupVec<4>(outConnectivity));
// Add cells to output cellset
outCellSet.Fill(cellSet.GetNumberOfPoints(), vtkm::CellShapeTagTetra::Id, 4, outConnectivity);
return outCellSet;
}
}
} // namespace vtkm::worklet