Commit e2ebdfea authored by Allison Vacanti's avatar Allison Vacanti

Reduce memory usage in vtkm cell converters.

Also updated the AverageTo... tests to test vtkmCellSetExplicit.

An API typo in the WarpScalar was fixed to match the updated
vtkm checkout.
parent e412ca55
......@@ -12,19 +12,68 @@ import vtk
from vtk.numpy_interface import dataset_adapter as dsa
from vtk.numpy_interface import algorithms as algs
rt = vtk.vtkRTAnalyticSource()
def test_dataset(ds):
p2c = vtk.vtkPointDataToCellData()
p2c.SetInputData(ds)
p2c.Update()
p2c = vtk.vtkPointDataToCellData()
p2c.SetInputConnection(rt.GetOutputPort())
p2c.Update()
d1 = dsa.WrapDataObject(p2c.GetOutput())
d1 = dsa.WrapDataObject(p2c.GetOutput())
vtkm_p2c = vtk.vtkmAverageToCells()
vtkm_p2c.SetInputData(ds)
vtkm_p2c.SetInputArrayToProcess(0, 0, 0, vtk.vtkDataObject.FIELD_ASSOCIATION_POINTS, "RTData")
vtkm_p2c.Update()
vtkm_p2c = vtk.vtkmAverageToCells()
vtkm_p2c.SetInputData(rt.GetOutput())
vtkm_p2c.SetInputArrayToProcess(0, 0, 0, vtk.vtkDataObject.FIELD_ASSOCIATION_POINTS, "RTData")
vtkm_p2c.Update()
d2 = dsa.WrapDataObject(vtkm_p2c.GetOutput())
d2 = dsa.WrapDataObject(vtkm_p2c.GetOutput())
rtD1 = d1.PointData['RTData']
rtD2 = d2.PointData['RTData']
assert (algs.max(algs.abs(d1.CellData['RTData'] - d2.CellData['RTData'])) < 10E-4)
assert (algs.max(algs.abs(rtD1 - rtD2)) < 10E-4)
print("Testing simple debugging grid...")
# This dataset matches the example values in vtkmCellSetExplicit:
dbg = vtk.vtkUnstructuredGrid()
dbg.SetPoints(vtk.vtkPoints())
dbg.GetPoints().SetNumberOfPoints(7)
dbg.InsertNextCell(vtk.VTK_TRIANGLE, 3, [0, 1, 2])
dbg.InsertNextCell(vtk.VTK_QUAD, 4, [0, 1, 3, 4])
dbg.InsertNextCell(vtk.VTK_TRIANGLE, 3, [1, 3, 5])
dbg.InsertNextCell(vtk.VTK_LINE, 2, [5, 6])
dbgRt = vtk.vtkDoubleArray()
dbgRt.SetNumberOfTuples(7)
dbgRt.SetName('RTData')
dbgRt.SetValue(0, 17.40)
dbgRt.SetValue(1, 123.0)
dbgRt.SetValue(2, 28.60)
dbgRt.SetValue(3, 19.47)
dbgRt.SetValue(4, 3.350)
dbgRt.SetValue(5, 0.212)
dbgRt.SetValue(6, 1023.)
dbg.GetPointData().AddArray(dbgRt)
test_dataset(dbg)
print("Success!")
print("Testing homogeneous image data...")
source = vtk.vtkRTAnalyticSource()
source.Update()
imgData = source.GetOutput()
test_dataset(imgData)
print("Success!")
d = dsa.WrapDataObject(imgData)
rtData = d.PointData['RTData']
rtMin = algs.min(rtData)
rtMax = algs.max(rtData)
clipScalar = 0.5 * (rtMin + rtMax)
print("Testing non-homogenous unstructured grid...")
clip = vtk.vtkClipDataSet()
clip.SetInputData(imgData)
clip.SetValue(clipScalar)
clip.Update()
ugrid = clip.GetOutput()
test_dataset(ugrid)
print("Success!")
......@@ -12,23 +12,72 @@ import vtk
from vtk.numpy_interface import dataset_adapter as dsa
from vtk.numpy_interface import algorithms as algs
rt = vtk.vtkRTAnalyticSource()
def test_dataset(ds):
p2c = vtk.vtkPointDataToCellData()
p2c.SetInputData(ds)
p2c.Update()
p2c = vtk.vtkPointDataToCellData()
p2c.SetInputConnection(rt.GetOutputPort())
p2c.Update()
c2p = vtk.vtkCellDataToPointData()
c2p.SetInputConnection(p2c.GetOutputPort())
c2p.Update()
c2p = vtk.vtkCellDataToPointData()
c2p.SetInputConnection(p2c.GetOutputPort())
c2p.Update()
d1 = dsa.WrapDataObject(c2p.GetOutput())
d1 = dsa.WrapDataObject(c2p.GetOutput())
c2p = vtk.vtkmAverageToPoints()
c2p.SetInputData(p2c.GetOutput())
c2p.SetInputArrayToProcess(0, 0, 0, vtk.vtkDataObject.FIELD_ASSOCIATION_CELLS, "RTData")
c2p.Update()
c2p = vtk.vtkmAverageToPoints()
c2p.SetInputData(p2c.GetOutput())
c2p.SetInputArrayToProcess(0, 0, 0, vtk.vtkDataObject.FIELD_ASSOCIATION_CELLS, "RTData")
c2p.Update()
d2 = dsa.WrapDataObject(c2p.GetOutput())
d2 = dsa.WrapDataObject(c2p.GetOutput())
rtD1 = d1.PointData['RTData']
rtD2 = d2.PointData['RTData']
assert (algs.max(algs.abs(d1.PointData['RTData'] - d2.PointData['RTData'])) < 10E-4)
assert (algs.max(algs.abs(rtD1 - rtD2)) < 10E-4)
print("Testing simple debugging grid...")
# This dataset matches the example values in vtkmCellSetExplicit:
dbg = vtk.vtkUnstructuredGrid()
dbg.SetPoints(vtk.vtkPoints())
dbg.GetPoints().SetNumberOfPoints(7)
dbg.InsertNextCell(vtk.VTK_TRIANGLE, 3, [0, 1, 2])
dbg.InsertNextCell(vtk.VTK_QUAD, 4, [0, 1, 3, 4])
dbg.InsertNextCell(vtk.VTK_TRIANGLE, 3, [1, 3, 5])
dbg.InsertNextCell(vtk.VTK_LINE, 2, [5, 6])
dbgRt = vtk.vtkDoubleArray()
dbgRt.SetNumberOfTuples(7)
dbgRt.SetName('RTData')
dbgRt.SetValue(0, 17.40)
dbgRt.SetValue(1, 123.0)
dbgRt.SetValue(2, 28.60)
dbgRt.SetValue(3, 19.47)
dbgRt.SetValue(4, 3.350)
dbgRt.SetValue(5, 0.212)
dbgRt.SetValue(6, 1023.)
dbg.GetPointData().AddArray(dbgRt)
test_dataset(dbg)
print("Success!")
print("Testing homogeneous image data...")
source = vtk.vtkRTAnalyticSource()
source.Update()
imgData = source.GetOutput()
test_dataset(imgData)
print("Success!")
d = dsa.WrapDataObject(imgData)
rtData = d.PointData['RTData']
rtMin = algs.min(rtData)
rtMax = algs.max(rtData)
clipScalar = 0.5 * (rtMin + rtMax)
print("Testing non-homogenous unstructured grid...")
clip = vtk.vtkClipDataSet()
clip.SetInputData(imgData)
clip.SetValue(clipScalar)
clip.Update()
ugrid = clip.GetOutput()
test_dataset(ugrid)
print("Success!")
......@@ -22,48 +22,187 @@
#include "vtkmCellSetExplicit.h"
#include "vtkmConnectivityExec.h"
#include <vtkm/cont/ArrayHandleCast.h>
#include <vtkm/cont/ArrayHandleCounting.h>
#include <vtkm/cont/ArrayHandleConstant.h>
#include <vtkm/cont/ArrayHandleImplicit.h>
#include <vtkm/cont/internal/ReverseConnectivityBuilder.h>
#include <vtkm/worklet/WorkletMapField.h>
#include <vtkm/worklet/DispatcherMapField.h>
namespace {
class ComputeReverseMapping : public vtkm::worklet::WorkletMapField
#include <utility>
namespace
{
// Converts [0, rconnSize) to [0, connSize), skipping cell length entries.
template <typename Device>
struct ExplicitRConnToConn
{
public:
typedef void ControlSignature(FieldIn<IdType> cellId,
FieldIn<IdType> offset,
WholeArrayIn<IdType> connectivity,
WholeArrayOut<IdType> pointIds,
WholeArrayOut<IdType> cellIds);
typedef void ExecutionSignature(_1, _2, _3, _4, _5);
typedef _1 InputDomain;
VTKM_SUPPRESS_EXEC_WARNINGS
template <typename ConnPortalType, typename PortalType>
VTKM_EXEC void operator()(const vtkm::Id& cellId,
const vtkm::Id& read_offset,
const ConnPortalType& connectivity,
const PortalType& pointIdKey,
const PortalType& pointIdValue) const
using OffsetsArray = vtkm::cont::ArrayHandle<vtkm::Id, tovtkm::vtkAOSArrayContainerTag>;
using OffsetsPortal = decltype(std::declval<OffsetsArray>().PrepareForInput(Device()));
// Functor that modifies the offsets array so we can compute point id indices
// Output is:
// modOffset[i] = offsets[i] - i
struct OffsetsModifier
{
OffsetsPortal Offsets;
VTKM_CONT
OffsetsModifier(const OffsetsPortal& offsets = OffsetsPortal{})
: Offsets(offsets)
{
}
VTKM_EXEC
vtkm::Id operator()(vtkm::Id inIdx) const
{
return this->Offsets.Get(inIdx) - inIdx;
}
};
using ModOffsetsArray = vtkm::cont::ArrayHandleImplicit<OffsetsModifier>;
using ModOffsetsPortal = decltype(std::declval<ModOffsetsArray>().PrepareForInput(Device()));
VTKM_CONT
ExplicitRConnToConn(const ModOffsetsPortal& offsets)
: Offsets(offsets)
{
//We can compute the correct write_offset by looking at the current
//read_offset and subtracting the number of cells we already processed
//which is equivalent to our cellid.
//This properly removes from read_offset the vtk cell padding that is
//added to the connectivity array
const vtkm::Id write_offset = read_offset - cellId;
const vtkm::Id numIndices = connectivity.Get(read_offset);
for (vtkm::Id i = 0; i < numIndices; i++)
}
VTKM_EXEC
vtkm::Id operator()(vtkm::Id rConnIdx) const
{
// Compute the connectivity array index (skipping cell length entries)
// The number of cell length entries can be found by taking the index of
// the upper bound of inIdx in Offsets and adding it to inIdx.
//
// Example:
// Conn: | 3 X X X | 4 X X X X | 3 X X X | 2 X X |
// Idx: | 0 1 2 3 | 4 5 6 7 8 | 9 10 11 12 | 13 14 15 |
// InIdx: 0 1 2 3 4 5 6 | 7 8 9 10 11
//
// ModOffset[i] = Offsets[i] - i:
// Offsets: 0 4 9 13 (16)
// ModOffsets: 0 3 7 10 (12)
//
// Define UB(x) to return the index of the upper bound of x in ModOffsets,
// the i'th point index's location in Conn is computed as:
// OutId = UB(InIdx) + InIdx
//
// This gives us the expected out indices:
// InIdx: 0 1 2 3 4 5 6 7 8 9 10 11
// UB(InIdx): 1 1 1 2 2 2 2 3 3 3 4 4
// OutIdx: 1 2 3 5 6 7 8 10 11 12 14 15
const vtkm::Id uBIdx = this->UpperBoundIdx(rConnIdx);
const vtkm::Id connIdx = rConnIdx + uBIdx;
return connIdx;
}
private:
ModOffsetsPortal Offsets;
VTKM_EXEC
vtkm::Id UpperBoundIdx(vtkm::Id inIdx) const
{
vtkm::Id first = 0;
vtkm::Id length = this->Offsets.GetNumberOfValues();
while (length > 0)
{
pointIdKey.Set(write_offset + i, connectivity.Get(read_offset + i + 1) );
pointIdValue.Set(write_offset + i, cellId);
vtkm::Id half = length / 2;
vtkm::Id pos = first + half;
vtkm::Id val = this->Offsets.Get(pos);
if (val <= inIdx)
{
first = pos + 1;
length -= half + 1;
}
else
{
length = half;
}
}
return first;
}
};
} //namespace
// Converts a connectivity index to a cell id:
template <typename Device>
struct ExplicitCellIdCalc
{
using OffsetsArray = vtkm::cont::ArrayHandle<vtkm::Id, tovtkm::vtkAOSArrayContainerTag>;
using OffsetsPortal = decltype(std::declval<OffsetsArray>().PrepareForInput(Device()));
vtkm::Id ConnSize;
OffsetsPortal Offsets;
VTKM_CONT
ExplicitCellIdCalc(vtkm::Id connSize, const OffsetsPortal& offsets)
: ConnSize(connSize)
, Offsets(offsets)
{
}
// Computes the cellid of the connectivity index i.
//
// For a mixed-cell connectivity, the offset table is used to compute
// the cell id.
//
// Example:
// Conn: | 3 X X X | 4 X X X X | 3 X X X | 2 X X |
// Idx: | 1 2 3 | 5 6 7 8 | 10 11 12 | 14 15 |
//
// Offsets: 0 4 9 13
// ModOffsets: 4 9 13 16
//
// Computing the index of the lower bound in ModOffsets for each Idx gives
// the expected cell id values:
// CellId: | 0 0 0 | 1 1 1 1 | 2 2 2 | 3 3 |
VTKM_EXEC
vtkm::Id operator()(vtkm::Id i) const
{
return this->LowerBound(i);
}
private:
/// Returns the i+1 offset, or the full size of the connectivity if at end.
VTKM_EXEC
vtkm::Id GetModifiedOffset(vtkm::Id i) const
{
++i;
return (i >= this->Offsets.GetNumberOfValues()) ? this->ConnSize
: this->Offsets.Get(i);
}
VTKM_EXEC
vtkm::Id LowerBound(vtkm::Id inVal) const
{
vtkm::Id first = 0;
vtkm::Id length = this->Offsets.GetNumberOfValues();
while (length > 0)
{
vtkm::Id half = length / 2;
vtkm::Id pos = first + half;
vtkm::Id val = this->GetModifiedOffset(pos);
if (val < inVal)
{
first = pos + 1;
length -= half + 1;
}
else
{
length = half;
}
}
return first;
}
};
} // end anon namespace
namespace vtkm {
namespace cont {
......@@ -127,38 +266,32 @@ typename vtkm::exec::ReverseConnectivityVTK<Device>
//#2 as it easier to construct
if(!this->ReverseConnectivityBuilt)
{
const vtkm::Id numberOfCells = this->GetNumberOfCells();
const vtkm::Id numberOfPoints = this->GetNumberOfPoints();
const vtkm::Id connectivityLength = this->Connectivity.GetNumberOfValues();
const vtkm::Id rconnSize = connectivityLength - this->IndexOffsets.GetNumberOfValues();
// create a mapping of where each key is the point id and the value
// is the cell id.
using Algorithm = vtkm::cont::DeviceAdapterAlgorithm<Device>;
vtkm::cont::ArrayHandle<vtkm::Id> pointIdKey;
// We need to allocate pointIdKey and RConn to correct length.
// which for this is equal to connectivityLength - len(this->IndexOffsets)
pointIdKey.Allocate(rconnSize);
this->RConn.Allocate(rconnSize);
vtkm::worklet::DispatcherMapField<ComputeReverseMapping, Device> dispatcher;
dispatcher.Invoke(vtkm::cont::make_ArrayHandleCounting<vtkm::Id>(0, 1, numberOfCells),
this->IndexOffsets, this->Connectivity, pointIdKey, this->RConn);
Algorithm::SortByKey(pointIdKey, this->RConn);
// now we can compute the NumIndices
vtkm::cont::ArrayHandle<vtkm::Id> reducedKeys;
Algorithm::ReduceByKey(pointIdKey,
vtkm::cont::make_ArrayHandleConstant(vtkm::IdComponent(1), rconnSize),
reducedKeys, this->RNumIndices, vtkm::Add());
// than a scan exclusive will give us the index offsets
using CastedNumIndices = vtkm::cont::ArrayHandleCast<vtkm::Id,
vtkm::cont::ArrayHandle<vtkm::IdComponent>>;
Algorithm::ScanExclusive(
CastedNumIndices(this->RNumIndices), this->RIndexOffsets);
this->NumberOfPoints = reducedKeys.GetNumberOfValues();
auto offsetsPortal = this->IndexOffsets.PrepareForInput(Device());
typename ExplicitRConnToConn<Device>::OffsetsModifier offsetModifier{offsetsPortal};
auto modOffsets = vtkm::cont::make_ArrayHandleImplicit(offsetModifier,
this->IndexOffsets.GetNumberOfValues());
const ExplicitRConnToConn<Device> rconnToConnCalc{modOffsets.PrepareForInput(Device())};
const ExplicitCellIdCalc<Device> cellIdCalc{this->Connectivity.GetNumberOfValues(),
this->IndexOffsets.PrepareForInput(Device())};
vtkm::cont::internal::ReverseConnectivityBuilder builder;
builder.Run(this->Connectivity,
this->RConn,
this->RNumIndices,
this->RIndexOffsets,
rconnToConnCalc,
cellIdCalc,
numberOfPoints,
rconnSize,
Device{});
this->NumberOfPoints = this->RIndexOffsets.GetNumberOfValues();
this->ReverseConnectivityBuilt = true;
}
......
......@@ -21,49 +21,38 @@
#include "vtkmCellSetSingleType.h"
#include "vtkmConnectivityExec.h"
#include <vtkm/cont/ArrayHandleCast.h>
#include <vtkm/cont/ArrayHandleCounting.h>
#include <vtkm/cont/ArrayHandleConstant.h>
#include <vtkm/cont/internal/ReverseConnectivityBuilder.h>
#include <vtkm/worklet/WorkletMapField.h>
#include <vtkm/worklet/DispatcherMapField.h>
namespace
{
namespace {
class ComputeSingleTypeReverseMapping : public vtkm::worklet::WorkletMapField
// Converts [0, rconnSize) to [0, connSize), skipping cell length entries.
struct SingleTypeRConnToConn
{
vtkm::Id NumberOfPointsPerCell;
public:
ComputeSingleTypeReverseMapping(vtkm::Id numberOfPointsPerCell):
NumberOfPointsPerCell(numberOfPointsPerCell)
{
vtkm::Id PointsPerCell;
VTKM_EXEC
vtkm::Id operator()(vtkm::Id rconnIdx) const
{
return rconnIdx + 1 + (rconnIdx / this->PointsPerCell);
}
typedef void ControlSignature(FieldIn<IdType> cellIndex,
WholeArrayIn<IdType> connectivity,
WholeArrayOut<IdType> pointIds,
WholeArrayOut<IdType> cellIds);
typedef void ExecutionSignature(_1, _2, _3, _4);
typedef _1 InputDomain;
VTKM_SUPPRESS_EXEC_WARNINGS
template <typename ConnPortalType, typename PortalType>
VTKM_EXEC void operator()(const vtkm::Id& cellId,
const ConnPortalType& connectivity,
const PortalType& pointIdKey,
const PortalType& pointIdValue) const
};
// Converts a connectivity index to a cell id:
struct SingleTypeCellIdCalc
{
vtkm::Id EncodedCellSize;
VTKM_EXEC
vtkm::Id operator()(vtkm::Id connIdx) const
{
const vtkm::Id read_offset = (NumberOfPointsPerCell+1) * cellId;
const vtkm::Id write_offset = NumberOfPointsPerCell * cellId;
const vtkm::Id numIndices = connectivity.Get(read_offset);
for (vtkm::Id i = 0; i < numIndices; i++)
{
pointIdKey.Set(write_offset + i, connectivity.Get(read_offset + i + 1) );
pointIdValue.Set(write_offset + i, cellId);
}
return connIdx / this->EncodedCellSize;
}
};
} //namespace
} // end anon namespace
namespace vtkm {
namespace cont {
......@@ -124,43 +113,29 @@ typename vtkm::exec::ReverseConnectivityVTK<Device>
{
if(!this->ReverseConnectivityBuilt)
{
const vtkm::Id numberOfPoints = this->GetNumberOfPoints();
const vtkm::Id numberOfCells = this->GetNumberOfCells();
const vtkm::Id numberOfPointsPerCell = this->DetermineNumberOfPoints();
const vtkm::Id rconnSize = numberOfCells*numberOfPointsPerCell;
// create a mapping of where each key is the point id and the value
// is the cell id.
using Algorithm = vtkm::cont::DeviceAdapterAlgorithm<Device>;
vtkm::cont::ArrayHandle<vtkm::Id> pointIdKey;
// We need to allocate pointIdKey and RConn to correct length.
// which for this is equal to numberOfCells * numberOfPointsPerCell
// as the connectivity has the vtk padding per cell
pointIdKey.Allocate(rconnSize);
this->RConn.Allocate(rconnSize);
vtkm::worklet::DispatcherMapField<ComputeSingleTypeReverseMapping, Device>
dispatcher(
ComputeSingleTypeReverseMapping(this->DetermineNumberOfPoints()));
dispatcher.Invoke(vtkm::cont::make_ArrayHandleCounting<vtkm::Id>(0, 1, numberOfCells),
this->Connectivity, pointIdKey, this->RConn);
Algorithm::SortByKey(pointIdKey, this->RConn);
// now we can compute the NumIndices
vtkm::cont::ArrayHandle<vtkm::Id> reducedKeys;
Algorithm::ReduceByKey(pointIdKey,
vtkm::cont::make_ArrayHandleConstant(vtkm::IdComponent(1), rconnSize),
reducedKeys, this->RNumIndices, vtkm::Add());
// than a scan exclusive will give us the index offsets
using CastedNumIndices = vtkm::cont::ArrayHandleCast<vtkm::Id,
vtkm::cont::ArrayHandle<vtkm::IdComponent>>;
Algorithm::ScanExclusive(
CastedNumIndices(this->RNumIndices), this->RIndexOffsets);
this->NumberOfPoints = reducedKeys.GetNumberOfValues();
const SingleTypeRConnToConn rconnToConnCalc{numberOfPointsPerCell};
const SingleTypeCellIdCalc cellIdCalc{numberOfPointsPerCell + 1}; // +1 for cell length entries
vtkm::cont::internal::ReverseConnectivityBuilder builder;
builder.Run(this->Connectivity,
this->RConn,
this->RNumIndices,
this->RIndexOffsets,
rconnToConnCalc,
cellIdCalc,
numberOfPoints,
rconnSize,
Device{});
this->NumberOfPoints = this->RIndexOffsets.GetNumberOfValues();
this->ReverseConnectivityBuilt = true;
}
} // End if !RConnBuilt
//no need to have a reverse shapes array, as everything has the shape type
//of vertex
......
......@@ -193,11 +193,11 @@ int vtkmWarpScalar::RequestData(vtkInformation* vtkNotUsed(request),
zValues.push_back(input->GetPoints()->GetPoint(i)[2]);
}
vtkm::cont::DataSetFieldAdd::AddPointField(in, "scalarfactor", zValues);
warpScalar.SetScarlarFactorField("scalarfactor");
warpScalar.SetScalarFactorField("scalarfactor");
}
else
{
warpScalar.SetScarlarFactorField(std::string(inScalars->GetName()));
warpScalar.SetScalarFactorField(std::string(inScalars->GetName()));
}
vtkmWarpScalarFilterPolicy policy;
......
Subproject commit 8077b031a8315c21ac0eb97d8e95cc36190d4963
Subproject commit a891e6d90a9f07a08ef4187a659c1f25c353668d
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment