5a473b0012
C has a feature where if you perform arithmetic on small integers (like char and short), it will automatically promote the result to a 32 bit integer. If you then store that back in the same type you started with GCC will warn you that you are loosing the precision (that you didn't ask for in the first place). This is particularly annoying in templated code. Anyway, fixed yet another instance of that happening.
718 lines
23 KiB
C++
718 lines
23 KiB
C++
//============================================================================
|
|
// 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 Sandia Corporation.
|
|
// Copyright 2014 UT-Battelle, LLC.
|
|
// Copyright 2014 Los Alamos National Security.
|
|
//
|
|
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
|
|
// 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 vtkm_m_worklet_Clip_h
|
|
#define vtkm_m_worklet_Clip_h
|
|
|
|
#include <vtkm/worklet/DispatcherMapTopology.h>
|
|
#include <vtkm/worklet/WorkletMapTopology.h>
|
|
#include <vtkm/worklet/internal/ClipTables.h>
|
|
|
|
#include <vtkm/cont/CellSetExplicit.h>
|
|
#include <vtkm/cont/CoordinateSystem.h>
|
|
#include <vtkm/cont/DeviceAdapterAlgorithm.h>
|
|
#include <vtkm/cont/DynamicArrayHandle.h>
|
|
#include <vtkm/cont/DynamicCellSet.h>
|
|
#include <vtkm/cont/Timer.h>
|
|
#include <vtkm/ImplicitFunctions.h>
|
|
|
|
#include <vtkm/exec/ExecutionWholeArray.h>
|
|
#include <vtkm/exec/FunctorBase.h>
|
|
|
|
#if defined(THRUST_MAJOR_VERSION) && THRUST_MAJOR_VERSION == 1 && \
|
|
THRUST_MINOR_VERSION == 8 && THRUST_SUBMINOR_VERSION < 3
|
|
// Workaround a bug in thrust 1.8.0 - 1.8.2 scan implementations which produces
|
|
// wrong results
|
|
#include <thrust/detail/type_traits.h>
|
|
#define THRUST_SCAN_WORKAROUND
|
|
#endif
|
|
|
|
|
|
namespace vtkm {
|
|
namespace worklet {
|
|
|
|
namespace internal {
|
|
|
|
// TODO: The following code is meant to be temporary until similar functionality is
|
|
// implemented elsewhere.
|
|
// The current DeviceAdapterAlgorithm::Copy resizes the destination array, and
|
|
// therefore, is not usable for our purpose here.
|
|
template <typename T, typename StorageIn, typename StorageOut, typename DeviceAdapter>
|
|
class CopyArray : public vtkm::exec::FunctorBase
|
|
{
|
|
public:
|
|
typedef typename vtkm::cont::ArrayHandle<T, StorageIn>
|
|
::template ExecutionTypes<DeviceAdapter>::PortalConst SourcePortalType;
|
|
typedef typename vtkm::cont::ArrayHandle<T, StorageOut>
|
|
::template ExecutionTypes<DeviceAdapter>::Portal DestPortalType;
|
|
|
|
VTKM_CONT_EXPORT
|
|
CopyArray(SourcePortalType input, DestPortalType output)
|
|
: Input(input), Output(output)
|
|
{
|
|
}
|
|
|
|
VTKM_EXEC_EXPORT
|
|
void operator()(vtkm::Id idx) const
|
|
{
|
|
this->Output.Set(idx, this->Input.Get(idx));
|
|
}
|
|
|
|
private:
|
|
SourcePortalType Input;
|
|
DestPortalType Output;
|
|
};
|
|
|
|
template <typename T, typename StorageIn, typename StorageOut, typename DeviceAdapter>
|
|
VTKM_CONT_EXPORT
|
|
void ResizeArrayHandle(const vtkm::cont::ArrayHandle<T, StorageIn> &input,
|
|
vtkm::Id size,
|
|
vtkm::cont::ArrayHandle<T, StorageOut> &output,
|
|
DeviceAdapter)
|
|
{
|
|
typedef vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapter> Algorithm;
|
|
output.Allocate(size);
|
|
CopyArray<T, StorageIn, StorageOut, DeviceAdapter>
|
|
copyArray(input.PrepareForInput(DeviceAdapter()),
|
|
output.PrepareForInPlace(DeviceAdapter()));
|
|
Algorithm::Schedule(copyArray, input.GetNumberOfValues());
|
|
}
|
|
|
|
template<typename T>
|
|
VTKM_EXEC_EXPORT
|
|
void swap(T &v1, T &v2)
|
|
{
|
|
T temp = v1;
|
|
v1 = v2;
|
|
v2 = temp;
|
|
}
|
|
|
|
template <typename T>
|
|
VTKM_EXEC_CONT_EXPORT
|
|
T Scale(const T &val, vtkm::Float64 scale)
|
|
{
|
|
return static_cast<T>(scale * static_cast<vtkm::Float64>(val));
|
|
}
|
|
|
|
template <typename T, vtkm::IdComponent NumComponents>
|
|
VTKM_EXEC_CONT_EXPORT
|
|
vtkm::Vec<T, NumComponents> Scale(const vtkm::Vec<T, NumComponents> &val,
|
|
vtkm::Float64 scale)
|
|
{
|
|
return val * scale;
|
|
}
|
|
|
|
template <typename DeviceAdapter>
|
|
class ExecutionConnectivityExplicit : vtkm::exec::ExecutionObjectBase
|
|
{
|
|
private:
|
|
typedef typename vtkm::cont::ArrayHandle<vtkm::UInt8>
|
|
::template ExecutionTypes<DeviceAdapter>::Portal UInt8Portal;
|
|
|
|
typedef typename vtkm::cont::ArrayHandle<vtkm::IdComponent>
|
|
::template ExecutionTypes<DeviceAdapter>::Portal IdComponentPortal;
|
|
|
|
typedef typename vtkm::cont::ArrayHandle<vtkm::Id>
|
|
::template ExecutionTypes<DeviceAdapter>::Portal IdPortal;
|
|
|
|
public:
|
|
VTKM_CONT_EXPORT
|
|
ExecutionConnectivityExplicit()
|
|
: Shapes(), NumIndices(), Connectivity(), IndexOffsets()
|
|
{
|
|
}
|
|
|
|
VTKM_CONT_EXPORT
|
|
ExecutionConnectivityExplicit(const UInt8Portal &shapes,
|
|
const IdComponentPortal &numIndices,
|
|
const IdPortal &connectivity,
|
|
const IdPortal &indexOffsets)
|
|
: Shapes(shapes),
|
|
NumIndices(numIndices),
|
|
Connectivity(connectivity),
|
|
IndexOffsets(indexOffsets)
|
|
{
|
|
}
|
|
|
|
VTKM_EXEC_EXPORT
|
|
void SetCellShape(vtkm::Id cellIndex, vtkm::UInt8 shape)
|
|
{
|
|
this->Shapes.Set(cellIndex, shape);
|
|
}
|
|
|
|
VTKM_EXEC_EXPORT
|
|
void SetNumberOfIndices(vtkm::Id cellIndex, vtkm::IdComponent numIndices)
|
|
{
|
|
this->NumIndices.Set(cellIndex, numIndices);
|
|
}
|
|
|
|
VTKM_EXEC_EXPORT
|
|
void SetIndexOffset(vtkm::Id cellIndex, vtkm::Id indexOffset)
|
|
{
|
|
this->IndexOffsets.Set(cellIndex, indexOffset);
|
|
}
|
|
|
|
VTKM_EXEC_EXPORT
|
|
void SetConnectivity(vtkm::Id connectivityIndex, vtkm::Id pointIndex)
|
|
{
|
|
this->Connectivity.Set(connectivityIndex, pointIndex);
|
|
}
|
|
|
|
private:
|
|
UInt8Portal Shapes;
|
|
IdComponentPortal NumIndices;
|
|
IdPortal Connectivity;
|
|
IdPortal IndexOffsets;
|
|
};
|
|
|
|
} // namespace internal
|
|
|
|
|
|
|
|
struct ClipStats
|
|
{
|
|
vtkm::Id NumberOfCells;
|
|
vtkm::Id NumberOfIndices;
|
|
vtkm::Id NumberOfNewPoints;
|
|
|
|
struct SumOp
|
|
{
|
|
VTKM_EXEC_CONT_EXPORT
|
|
ClipStats operator()(const ClipStats &cs1, const ClipStats &cs2) const
|
|
{
|
|
ClipStats sum = cs1;
|
|
sum.NumberOfCells += cs2.NumberOfCells;
|
|
sum.NumberOfIndices += cs2.NumberOfIndices;
|
|
sum.NumberOfNewPoints += cs2.NumberOfNewPoints;
|
|
return sum;
|
|
}
|
|
};
|
|
};
|
|
|
|
struct EdgeInterpolation
|
|
{
|
|
vtkm::Id Vertex1, Vertex2;
|
|
vtkm::Float64 Weight;
|
|
|
|
struct LessThanOp
|
|
{
|
|
VTKM_EXEC_EXPORT
|
|
bool operator()(const EdgeInterpolation &v1, const EdgeInterpolation &v2) const
|
|
{
|
|
return (v1.Vertex1 < v2.Vertex1) ||
|
|
(v1.Vertex1 == v2.Vertex1 && v1.Vertex2 < v2.Vertex2);
|
|
}
|
|
};
|
|
|
|
struct EqualToOp
|
|
{
|
|
VTKM_EXEC_EXPORT
|
|
bool operator()(const EdgeInterpolation &v1, const EdgeInterpolation &v2) const
|
|
{
|
|
return v1.Vertex1 == v2.Vertex1 && v1.Vertex2 == v2.Vertex2;
|
|
}
|
|
};
|
|
};
|
|
|
|
|
|
VTKM_EXEC_CONT_EXPORT
|
|
std::ostream& operator<<(std::ostream &out, const ClipStats &val)
|
|
{
|
|
return out << std::endl << "{ Cells: " << val.NumberOfCells
|
|
<< ", Indices: " << val.NumberOfIndices
|
|
<< ", NewPoints: " << val.NumberOfNewPoints << " }";
|
|
}
|
|
|
|
VTKM_EXEC_CONT_EXPORT
|
|
std::ostream& operator<<(std::ostream &out, const EdgeInterpolation &val)
|
|
{
|
|
return out << std::endl << "{ Vertex1: " << val.Vertex1
|
|
<< ", Vertex2: " << val.Vertex2
|
|
<< ", Weight: " << val.Weight << " }";
|
|
}
|
|
|
|
|
|
template<typename DeviceAdapter>
|
|
class Clip
|
|
{
|
|
private:
|
|
typedef internal::ClipTables::DevicePortal<DeviceAdapter> ClipTablesPortal;
|
|
|
|
typedef typename vtkm::cont::ArrayHandle<vtkm::Id>
|
|
::template ExecutionTypes<DeviceAdapter>::Portal IdPortal;
|
|
typedef typename vtkm::cont::ArrayHandle<vtkm::Id>
|
|
::template ExecutionTypes<DeviceAdapter>::PortalConst IdPortalConst;
|
|
typedef typename vtkm::cont::ArrayHandle<EdgeInterpolation>
|
|
::template ExecutionTypes<DeviceAdapter>::Portal EdgeInterpolationPortal;
|
|
typedef typename vtkm::cont::ArrayHandle<EdgeInterpolation>
|
|
::template ExecutionTypes<DeviceAdapter>::PortalConst
|
|
EdgeInterpolationPortalConst;
|
|
|
|
typedef vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapter> Algorithm;
|
|
|
|
|
|
public:
|
|
struct TypeClipStats : vtkm::ListTagBase<ClipStats> { };
|
|
|
|
|
|
class ComputeStats : public vtkm::worklet::WorkletMapPointToCell
|
|
{
|
|
public:
|
|
typedef void ControlSignature(CellSetIn cellset,
|
|
FieldInPoint<ScalarAll> scalars,
|
|
FieldOutCell<IdType> clipTableIdxs,
|
|
FieldOutCell<TypeClipStats> stats);
|
|
typedef void ExecutionSignature(_2, CellShape, PointCount, _3, _4);
|
|
|
|
VTKM_CONT_EXPORT
|
|
ComputeStats(vtkm::Float64 value, const ClipTablesPortal &clipTables)
|
|
: Value(value), ClipTables(clipTables)
|
|
{
|
|
}
|
|
|
|
template<typename ScalarsVecType, typename CellShapeTag>
|
|
VTKM_EXEC_EXPORT
|
|
void operator()(const ScalarsVecType &scalars,
|
|
CellShapeTag shape,
|
|
vtkm::Id count,
|
|
vtkm::Id& clipTableIdx,
|
|
ClipStats& stats) const
|
|
{
|
|
const vtkm::Id mask[] = { 1, 2, 4, 8, 16, 32, 64, 128 };
|
|
|
|
vtkm::Id caseId = 0;
|
|
for (vtkm::IdComponent i = 0; i < count; ++i)
|
|
{
|
|
caseId |= (static_cast<vtkm::Float64>(scalars[i]) > this->Value) ?
|
|
mask[i] : 0;
|
|
}
|
|
|
|
vtkm::Id idx = this->ClipTables.GetCaseIndex(shape.Id, caseId);
|
|
clipTableIdx = idx;
|
|
|
|
vtkm::Id numberOfCells = this->ClipTables.ValueAt(idx++);
|
|
vtkm::Id numberOfIndices = 0;
|
|
vtkm::Id numberOfNewPoints = 0;
|
|
for (vtkm::Id cell = 0; cell < numberOfCells; ++cell)
|
|
{
|
|
++idx; // skip shape-id
|
|
vtkm::Id npts = this->ClipTables.ValueAt(idx++);
|
|
numberOfIndices += npts;
|
|
while (npts--)
|
|
{
|
|
// value < 100 means a new point needs to be generated by clipping an edge
|
|
numberOfNewPoints += (this->ClipTables.ValueAt(idx++) < 100) ? 1 : 0;
|
|
}
|
|
}
|
|
|
|
stats.NumberOfCells = numberOfCells;
|
|
stats.NumberOfIndices = numberOfIndices;
|
|
stats.NumberOfNewPoints = numberOfNewPoints;
|
|
}
|
|
|
|
private:
|
|
vtkm::Float64 Value;
|
|
ClipTablesPortal ClipTables;
|
|
};
|
|
|
|
|
|
class GenerateCellSet : public vtkm::worklet::WorkletMapPointToCell
|
|
{
|
|
public:
|
|
typedef void ControlSignature(CellSetIn cellset,
|
|
FieldInPoint<ScalarAll> scalars,
|
|
FieldInCell<IdType> clipTableIdxs,
|
|
FieldInCell<TypeClipStats> cellSetIdxs,
|
|
ExecObject connectivityExplicit,
|
|
ExecObject interpolation,
|
|
ExecObject newPointsConnectivityReverseMap);
|
|
typedef void ExecutionSignature(CellShape, _2, FromIndices, _3, _4, _5, _6, _7);
|
|
|
|
VTKM_CONT_EXPORT
|
|
GenerateCellSet(vtkm::Float64 value, const ClipTablesPortal clipTables)
|
|
: Value(value), ClipTables(clipTables)
|
|
{
|
|
}
|
|
|
|
template<typename CellShapeTag,
|
|
typename ScalarsVecType,
|
|
typename IndicesVecType,
|
|
typename Storage>
|
|
VTKM_EXEC_EXPORT
|
|
void operator()(
|
|
CellShapeTag shape,
|
|
const ScalarsVecType &scalars,
|
|
const IndicesVecType &indices,
|
|
vtkm::Id clipTableIdx,
|
|
ClipStats cellSetIndices,
|
|
internal::ExecutionConnectivityExplicit<DeviceAdapter> connectivityExplicit,
|
|
vtkm::exec::ExecutionWholeArray<EdgeInterpolation, Storage, DeviceAdapter>
|
|
interpolation,
|
|
vtkm::exec::ExecutionWholeArray<vtkm::Id, Storage, DeviceAdapter>
|
|
newPointsConnectivityReverseMap) const
|
|
{
|
|
vtkm::Id idx = clipTableIdx;
|
|
|
|
// index of first cell
|
|
vtkm::Id cellIdx = cellSetIndices.NumberOfCells;
|
|
// index of first cell in connectivity array
|
|
vtkm::Id connectivityIdx = cellSetIndices.NumberOfIndices;
|
|
// index of new points generated by first cell
|
|
vtkm::Id newPtsIdx = cellSetIndices.NumberOfNewPoints;
|
|
|
|
vtkm::Id numberOfCells = this->ClipTables.ValueAt(idx++);
|
|
for (vtkm::Id cell = 0; cell < numberOfCells; ++cell, ++cellIdx)
|
|
{
|
|
connectivityExplicit.SetCellShape(cellIdx, this->ClipTables.ValueAt(idx++));
|
|
vtkm::IdComponent numPoints = this->ClipTables.ValueAt(idx++);
|
|
connectivityExplicit.SetNumberOfIndices(cellIdx, numPoints);
|
|
connectivityExplicit.SetIndexOffset(cellIdx, connectivityIdx);
|
|
|
|
for (vtkm::Id pt = 0; pt < numPoints; ++pt, ++idx)
|
|
{
|
|
vtkm::IdComponent entry =
|
|
static_cast<vtkm::IdComponent>(this->ClipTables.ValueAt(idx));
|
|
if (entry >= 100) // existing point
|
|
{
|
|
connectivityExplicit.SetConnectivity(connectivityIdx++,
|
|
indices[entry - 100]);
|
|
}
|
|
else // edge, new point to be generated by cutting the egde
|
|
{
|
|
internal::ClipTables::EdgeVec edge =
|
|
this->ClipTables.GetEdge(shape.Id, entry);
|
|
|
|
EdgeInterpolation ei;
|
|
ei.Vertex1 = indices[edge[0]];
|
|
ei.Vertex2 = indices[edge[1]];
|
|
if (ei.Vertex1 > ei.Vertex2)
|
|
{
|
|
internal::swap(ei.Vertex1, ei.Vertex2);
|
|
internal::swap(edge[0], edge[1]);
|
|
}
|
|
ei.Weight = (static_cast<vtkm::Float64>(scalars[edge[0]]) - this->Value) /
|
|
static_cast<vtkm::Float64>(scalars[edge[0]] - scalars[edge[1]]);
|
|
|
|
interpolation.Set(newPtsIdx , ei);
|
|
newPointsConnectivityReverseMap.Set(newPtsIdx, connectivityIdx++);
|
|
++newPtsIdx;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
private:
|
|
vtkm::Float64 Value;
|
|
ClipTablesPortal ClipTables;
|
|
};
|
|
|
|
|
|
// The following can be done using DeviceAdapterAlgorithm::LowerBounds followed by
|
|
// a worklet for updating connectivity. We are going with a custom worklet, that
|
|
// combines lower-bounds computation and connectivity update, because this is
|
|
// currently faster and uses less memory.
|
|
class AmendConnectivity : public vtkm::exec::FunctorBase
|
|
{
|
|
public:
|
|
VTKM_CONT_EXPORT
|
|
AmendConnectivity(EdgeInterpolationPortalConst newPoints,
|
|
EdgeInterpolationPortalConst uniqueNewPoints,
|
|
IdPortalConst newPointsConnectivityReverseMap,
|
|
vtkm::Id newPointsOffset,
|
|
IdPortal connectivity)
|
|
: NewPoints(newPoints),
|
|
UniqueNewPoints(uniqueNewPoints),
|
|
NewPointsConnectivityReverseMap(newPointsConnectivityReverseMap),
|
|
NewPointsOffset(newPointsOffset),
|
|
Connectivity(connectivity)
|
|
{
|
|
}
|
|
|
|
VTKM_EXEC_EXPORT
|
|
void operator()(vtkm::Id idx) const
|
|
{
|
|
EdgeInterpolation current = this->NewPoints.Get(idx);
|
|
typename EdgeInterpolation::LessThanOp lt;
|
|
|
|
// find point index by looking up in the unique points array (binary search)
|
|
vtkm::Id count = UniqueNewPoints.GetNumberOfValues();
|
|
vtkm::Id first = 0;
|
|
while (count > 0)
|
|
{
|
|
vtkm::Id step = count/2;
|
|
vtkm::Id mid = first + step;
|
|
if (lt(this->UniqueNewPoints.Get(mid), current))
|
|
{
|
|
first = ++mid;
|
|
count -= step + 1;
|
|
}
|
|
else
|
|
{
|
|
count = step;
|
|
}
|
|
}
|
|
|
|
this->Connectivity.Set(this->NewPointsConnectivityReverseMap.Get(idx),
|
|
this->NewPointsOffset + first);
|
|
}
|
|
|
|
private:
|
|
EdgeInterpolationPortalConst NewPoints;
|
|
EdgeInterpolationPortalConst UniqueNewPoints;
|
|
IdPortalConst NewPointsConnectivityReverseMap;
|
|
vtkm::Id NewPointsOffset;
|
|
IdPortal Connectivity;
|
|
};
|
|
|
|
|
|
|
|
Clip()
|
|
: ClipTablesInstance(), NewPointsInterpolation(), NewPointsOffset()
|
|
{
|
|
}
|
|
|
|
template <typename ScalarsArrayHandle>
|
|
vtkm::cont::CellSetExplicit<> Run(const vtkm::cont::DynamicCellSet &cellSet,
|
|
const ScalarsArrayHandle &scalars,
|
|
vtkm::Float64 value)
|
|
{
|
|
DeviceAdapter device;
|
|
ClipTablesPortal clipTablesDevicePortal =
|
|
this->ClipTablesInstance.GetDevicePortal(device);
|
|
|
|
|
|
// Step 1. compute counts for the elements of the cell set data structure
|
|
vtkm::cont::ArrayHandle<vtkm::Id> clipTableIdxs;
|
|
vtkm::cont::ArrayHandle<ClipStats> stats;
|
|
|
|
ComputeStats computeStats(value, clipTablesDevicePortal);
|
|
DispatcherMapTopology<ComputeStats>(computeStats).Invoke(cellSet, scalars,
|
|
clipTableIdxs, stats);
|
|
|
|
// compute offsets for each invocation
|
|
ClipStats zero = { 0, 0, 0 };
|
|
vtkm::cont::ArrayHandle<ClipStats> cellSetIndices;
|
|
ClipStats total = Algorithm::ScanExclusive(stats, cellSetIndices,
|
|
ClipStats::SumOp(), zero);
|
|
stats.ReleaseResources();
|
|
|
|
|
|
// Step 2. generate the output cell set
|
|
vtkm::cont::ArrayHandle<vtkm::UInt8> shapes;
|
|
vtkm::cont::ArrayHandle<vtkm::IdComponent> numIndices;
|
|
vtkm::cont::ArrayHandle<vtkm::Id> connectivity;
|
|
vtkm::cont::ArrayHandle<vtkm::Id> cellToConnectivityMap;
|
|
internal::ExecutionConnectivityExplicit<DeviceAdapter> outConnectivity(
|
|
shapes.PrepareForOutput(total.NumberOfCells, device),
|
|
numIndices.PrepareForOutput(total.NumberOfCells, device),
|
|
connectivity.PrepareForOutput(total.NumberOfIndices, device),
|
|
cellToConnectivityMap.PrepareForOutput(total.NumberOfCells, device));
|
|
|
|
vtkm::cont::ArrayHandle<EdgeInterpolation> newPoints;
|
|
// reverse map from the new points to connectivity array
|
|
vtkm::cont::ArrayHandle<vtkm::Id> newPointsConnectivityReverseMap;
|
|
|
|
GenerateCellSet generateCellSet(value, clipTablesDevicePortal);
|
|
|
|
DispatcherMapTopology<GenerateCellSet>(generateCellSet).Invoke(cellSet, scalars,
|
|
clipTableIdxs, cellSetIndices, outConnectivity,
|
|
vtkm::exec::ExecutionWholeArray<
|
|
EdgeInterpolation,
|
|
VTKM_DEFAULT_STORAGE_TAG,
|
|
DeviceAdapter>(newPoints, total.NumberOfNewPoints),
|
|
vtkm::exec::ExecutionWholeArray<
|
|
vtkm::Id,
|
|
VTKM_DEFAULT_STORAGE_TAG,
|
|
DeviceAdapter>(newPointsConnectivityReverseMap, total.NumberOfNewPoints));
|
|
cellSetIndices.ReleaseResources();
|
|
|
|
// Step 3. remove duplicates from the list of new points
|
|
vtkm::cont::ArrayHandle<vtkm::worklet::EdgeInterpolation> uniqueNewPoints;
|
|
|
|
Algorithm::SortByKey(newPoints, newPointsConnectivityReverseMap,
|
|
EdgeInterpolation::LessThanOp());
|
|
Algorithm::Copy(newPoints, uniqueNewPoints);
|
|
Algorithm::Unique(uniqueNewPoints, EdgeInterpolation::EqualToOp());
|
|
|
|
this->NewPointsInterpolation = uniqueNewPoints;
|
|
this->NewPointsOffset = scalars.GetNumberOfValues();
|
|
|
|
|
|
// Step 4. update the connectivity array with indexes to the new, unique points
|
|
AmendConnectivity computeNewPointsConnectivity(
|
|
newPoints.PrepareForInput(device),
|
|
uniqueNewPoints.PrepareForInput(device),
|
|
newPointsConnectivityReverseMap.PrepareForInput(device),
|
|
this->NewPointsOffset,
|
|
connectivity.PrepareForInPlace(device));
|
|
Algorithm::Schedule(computeNewPointsConnectivity, total.NumberOfNewPoints);
|
|
|
|
|
|
vtkm::cont::CellSetExplicit<> output;
|
|
output.Fill(shapes, numIndices, connectivity);
|
|
|
|
return output;
|
|
}
|
|
|
|
template<typename ImplicitFunction>
|
|
class ClipWithImplicitFunction
|
|
{
|
|
public:
|
|
VTKM_CONT_EXPORT
|
|
ClipWithImplicitFunction(Clip *clipper, const vtkm::cont::DynamicCellSet &cellSet,
|
|
ImplicitFunction function, vtkm::cont::CellSetExplicit<> *result)
|
|
: Clipper(clipper), CellSet(&cellSet), Function(function), Result(result)
|
|
{ }
|
|
|
|
template <typename ArrayHandleType>
|
|
VTKM_CONT_EXPORT
|
|
void operator()(const ArrayHandleType &handle) const
|
|
{
|
|
vtkm::cont::ArrayHandleTransform<
|
|
vtkm::FloatDefault,
|
|
ArrayHandleType,
|
|
vtkm::ImplicitFunctionValue<ImplicitFunction> > clipScalars(handle,
|
|
this->Function);
|
|
|
|
*this->Result = this->Clipper->Run(*this->CellSet, clipScalars, 0.0);
|
|
}
|
|
|
|
private:
|
|
Clip *Clipper;
|
|
const vtkm::cont::DynamicCellSet *CellSet;
|
|
const vtkm::ImplicitFunctionValue<ImplicitFunction> Function;
|
|
vtkm::cont::CellSetExplicit<> *Result;
|
|
};
|
|
|
|
template <typename ImplicitFunction>
|
|
vtkm::cont::CellSetExplicit<> Run(const vtkm::cont::DynamicCellSet &cellSet,
|
|
const ImplicitFunction &clipFunction,
|
|
const vtkm::cont::CoordinateSystem &coords)
|
|
{
|
|
vtkm::cont::CellSetExplicit<> output;
|
|
ClipWithImplicitFunction<ImplicitFunction> clip(this, cellSet, clipFunction,
|
|
&output);
|
|
coords.GetData().CastAndCall(clip);
|
|
|
|
return output;
|
|
}
|
|
|
|
class InterpolateField
|
|
{
|
|
public:
|
|
template <typename T>
|
|
class Kernel : public vtkm::exec::FunctorBase
|
|
{
|
|
public:
|
|
typedef typename vtkm::cont::ArrayHandle<T>
|
|
::template ExecutionTypes<DeviceAdapter>::Portal FieldPortal;
|
|
|
|
VTKM_CONT_EXPORT
|
|
Kernel(EdgeInterpolationPortalConst interpolation,
|
|
vtkm::Id newPointsOffset,
|
|
FieldPortal field)
|
|
: Interpolation(interpolation), NewPointsOffset(newPointsOffset), Field(field)
|
|
{
|
|
}
|
|
|
|
VTKM_EXEC_EXPORT
|
|
void operator()(vtkm::Id idx) const
|
|
{
|
|
EdgeInterpolation ei = this->Interpolation.Get(idx);
|
|
T v1 = Field.Get(ei.Vertex1);
|
|
T v2 = Field.Get(ei.Vertex2);
|
|
Field.Set(this->NewPointsOffset + idx,
|
|
static_cast<T>(internal::Scale(T(v2 - v1), ei.Weight) + v1));
|
|
}
|
|
|
|
private:
|
|
EdgeInterpolationPortalConst Interpolation;
|
|
vtkm::Id NewPointsOffset;
|
|
FieldPortal Field;
|
|
};
|
|
|
|
template <typename T, typename Storage>
|
|
VTKM_CONT_EXPORT
|
|
void operator()(const vtkm::cont::ArrayHandle<T, Storage> &field) const
|
|
{
|
|
vtkm::Id count = this->InterpolationArray.GetNumberOfValues();
|
|
|
|
vtkm::cont::ArrayHandle<T> result;
|
|
internal::ResizeArrayHandle(field, field.GetNumberOfValues() + count,
|
|
result, DeviceAdapter());
|
|
|
|
Kernel<T> kernel(
|
|
this->InterpolationArray.PrepareForInput(DeviceAdapter()),
|
|
this->NewPointsOffset,
|
|
result.PrepareForInPlace(DeviceAdapter()));
|
|
|
|
Algorithm::Schedule(kernel, count);
|
|
*(this->Output) = vtkm::cont::DynamicArrayHandle(result);
|
|
}
|
|
|
|
VTKM_CONT_EXPORT
|
|
InterpolateField(
|
|
vtkm::cont::ArrayHandle<EdgeInterpolation> interpolationArray,
|
|
vtkm::Id newPointsOffset,
|
|
vtkm::cont::DynamicArrayHandle *output)
|
|
: InterpolationArray(interpolationArray),
|
|
NewPointsOffset(newPointsOffset),
|
|
Output(output)
|
|
{
|
|
}
|
|
|
|
private:
|
|
vtkm::cont::ArrayHandle<EdgeInterpolation> InterpolationArray;
|
|
vtkm::Id NewPointsOffset;
|
|
vtkm::cont::DynamicArrayHandle *Output;
|
|
};
|
|
|
|
template <typename DynamicArrayHandleType>
|
|
vtkm::cont::DynamicArrayHandle ProcessField(
|
|
const DynamicArrayHandleType &fieldData) const
|
|
{
|
|
vtkm::cont::DynamicArrayHandle output;
|
|
fieldData.CastAndCall(InterpolateField(this->NewPointsInterpolation,
|
|
this->NewPointsOffset,
|
|
&output));
|
|
return output;
|
|
}
|
|
|
|
private:
|
|
internal::ClipTables ClipTablesInstance;
|
|
vtkm::cont::ArrayHandle<EdgeInterpolation> NewPointsInterpolation;
|
|
vtkm::Id NewPointsOffset;
|
|
};
|
|
|
|
}
|
|
} // namespace vtkm::worklet
|
|
|
|
#if defined(THRUST_SCAN_WORKAROUND)
|
|
namespace thrust {
|
|
namespace detail {
|
|
|
|
// causes a different code path which does not have the bug
|
|
template<> struct is_integral<vtkm::worklet::ClipStats> : public true_type {};
|
|
|
|
}
|
|
} // namespace thrust::detail
|
|
#endif
|
|
|
|
#endif // vtkm_m_worklet_Clip_h
|