b1343474c1
For no particularly good reason, there were two functions that converted and array of counts to an array of offsets: `ConvertNumComponentsToOffsets` and `ConvertNumIndicesToOffsets`. These functions were identical, except one was defined in `ArrayHandleGroupVecVariable.h` and the other was defined in `CellSetExplicit.h`. These two functions have been consolidated into one (which is now called `ConvertNumComponentsToOffsets`). The consolidated function has also been put in its own header file: `ConvertNumComponentsToOffsets.h`. Normally, backward compatibility would be established using deprecated features. However, one of the things being worked on is the removal of device-specific code (e.g. `vtkm::cont::Algorithm`) from core classes like `CellSetExplicit` so that less code needs to use the device compiler (especially downstream code). Part of this change removed unnecessary includes of `Algorithm.h` in `ArrayHandleGroupVecVariable.h` and `CellSetExplicit.h`. This header had to be added to some classes that were not including it themselves.
2427 lines
90 KiB
C++
2427 lines
90 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.
|
|
//
|
|
//=============================================================================
|
|
|
|
#ifndef vtk_m_worklet_MIR_h
|
|
#define vtk_m_worklet_MIR_h
|
|
|
|
#include <vtkm/CellShape.h>
|
|
#include <vtkm/Types.h>
|
|
#include <vtkm/VectorAnalysis.h>
|
|
|
|
#include <vtkm/cont/Algorithm.h>
|
|
#include <vtkm/cont/ArrayCopy.h>
|
|
#include <vtkm/cont/ArrayHandle.h>
|
|
#include <vtkm/cont/ArrayHandlePermutation.h>
|
|
#include <vtkm/cont/ArrayHandleView.h>
|
|
#include <vtkm/cont/CellSetExplicit.h>
|
|
#include <vtkm/cont/CellSetPermutation.h>
|
|
#include <vtkm/cont/ConvertNumComponentsToOffsets.h>
|
|
#include <vtkm/cont/CoordinateSystem.h>
|
|
#include <vtkm/cont/DynamicCellSet.h>
|
|
#include <vtkm/cont/ErrorFilterExecution.h>
|
|
#include <vtkm/cont/ExecutionObjectBase.h>
|
|
#include <vtkm/cont/Timer.h>
|
|
#include <vtkm/cont/VariantArrayHandle.h>
|
|
|
|
#include <vtkm/exec/FunctorBase.h>
|
|
|
|
#include <vtkm/worklet/DispatcherMapField.h>
|
|
#include <vtkm/worklet/DispatcherMapTopology.h>
|
|
#include <vtkm/worklet/DispatcherReduceByKey.h>
|
|
#include <vtkm/worklet/Keys.h>
|
|
#include <vtkm/worklet/ScatterCounting.h>
|
|
#include <vtkm/worklet/WorkletMapField.h>
|
|
#include <vtkm/worklet/WorkletMapTopology.h>
|
|
#include <vtkm/worklet/WorkletReduceByKey.h>
|
|
|
|
#include <vtkm/worklet/clip/ClipTables.h>
|
|
#include <vtkm/worklet/mir/MIRTables.h>
|
|
|
|
namespace vtkm
|
|
{
|
|
namespace worklet
|
|
{
|
|
|
|
struct MIRStats
|
|
{
|
|
vtkm::Id NumberOfCells = 0;
|
|
vtkm::Id NumberOfIndices = 0;
|
|
vtkm::Id NumberOfEdgeIndices = 0;
|
|
|
|
//VTKM New point stats
|
|
vtkm::Id NumberOfInCellPoints = 0;
|
|
vtkm::Id NumberOfInCellIndices = 0;
|
|
vtkm::Id NumberOfInCellInterpPoints = 0;
|
|
vtkm::Id NumberOfInCellEdgeIndices = 0;
|
|
|
|
struct SumOp
|
|
{
|
|
VTKM_EXEC_CONT
|
|
MIRStats operator()(const MIRStats& stat1, const MIRStats& stat2) const
|
|
{
|
|
MIRStats sum = stat1;
|
|
sum.NumberOfCells += stat2.NumberOfCells;
|
|
sum.NumberOfIndices += stat2.NumberOfIndices;
|
|
sum.NumberOfEdgeIndices += stat2.NumberOfEdgeIndices;
|
|
sum.NumberOfInCellPoints += stat2.NumberOfInCellPoints;
|
|
sum.NumberOfInCellIndices += stat2.NumberOfInCellIndices;
|
|
sum.NumberOfInCellInterpPoints += stat2.NumberOfInCellInterpPoints;
|
|
sum.NumberOfInCellEdgeIndices += stat2.NumberOfInCellEdgeIndices;
|
|
return sum;
|
|
}
|
|
};
|
|
};
|
|
struct EdgeInterpolation
|
|
{
|
|
vtkm::Id Vertex1 = -1;
|
|
vtkm::Id Vertex2 = -1;
|
|
vtkm::Float64 Weight = 0;
|
|
|
|
struct LessThanOp
|
|
{
|
|
VTKM_EXEC
|
|
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
|
|
bool operator()(const EdgeInterpolation& v1, const EdgeInterpolation& v2) const
|
|
{
|
|
return v1.Vertex1 == v2.Vertex1 && v1.Vertex2 == v2.Vertex2;
|
|
}
|
|
};
|
|
};
|
|
namespace MIRinternal
|
|
{
|
|
template <typename T>
|
|
VTKM_EXEC_CONT 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 vtkm::Vec<T, NumComponents> Scale(const vtkm::Vec<T, NumComponents>& val,
|
|
vtkm::Float64 scale)
|
|
{
|
|
return val * scale;
|
|
}
|
|
}
|
|
|
|
class ExecutionConnectivityExplicit
|
|
{
|
|
private:
|
|
using UInt8Portal = typename vtkm::cont::ArrayHandle<vtkm::UInt8>::WritePortalType;
|
|
using IdComponentPortal = typename vtkm::cont::ArrayHandle<vtkm::IdComponent>::WritePortalType;
|
|
using IdPortal = typename vtkm::cont::ArrayHandle<vtkm::Id>::WritePortalType;
|
|
|
|
public:
|
|
VTKM_CONT
|
|
ExecutionConnectivityExplicit() = default;
|
|
|
|
VTKM_CONT
|
|
ExecutionConnectivityExplicit(vtkm::cont::ArrayHandle<vtkm::UInt8> shapes,
|
|
vtkm::cont::ArrayHandle<vtkm::IdComponent> numberOfIndices,
|
|
vtkm::cont::ArrayHandle<vtkm::Id> connectivity,
|
|
vtkm::cont::ArrayHandle<vtkm::Id> offsets,
|
|
MIRStats stats,
|
|
vtkm::cont::DeviceAdapterId device,
|
|
vtkm::cont::Token& token)
|
|
: Shapes(shapes.PrepareForOutput(stats.NumberOfCells, device, token))
|
|
, NumberOfIndices(numberOfIndices.PrepareForOutput(stats.NumberOfCells, device, token))
|
|
, Connectivity(connectivity.PrepareForOutput(stats.NumberOfIndices, device, token))
|
|
, Offsets(offsets.PrepareForOutput(stats.NumberOfCells, device, token))
|
|
{
|
|
}
|
|
|
|
VTKM_EXEC
|
|
void SetCellShape(vtkm::Id cellIndex, vtkm::UInt8 shape) { this->Shapes.Set(cellIndex, shape); }
|
|
|
|
VTKM_EXEC
|
|
void SetNumberOfIndices(vtkm::Id cellIndex, vtkm::IdComponent numIndices)
|
|
{
|
|
this->NumberOfIndices.Set(cellIndex, numIndices);
|
|
}
|
|
|
|
VTKM_EXEC
|
|
void SetIndexOffset(vtkm::Id cellIndex, vtkm::Id indexOffset)
|
|
{
|
|
this->Offsets.Set(cellIndex, indexOffset);
|
|
}
|
|
|
|
VTKM_EXEC
|
|
void SetConnectivity(vtkm::Id connectivityIndex, vtkm::Id pointIndex)
|
|
{
|
|
this->Connectivity.Set(connectivityIndex, pointIndex);
|
|
}
|
|
|
|
private:
|
|
UInt8Portal Shapes;
|
|
IdComponentPortal NumberOfIndices;
|
|
IdPortal Connectivity;
|
|
IdPortal Offsets;
|
|
};
|
|
|
|
class ConnectivityExplicit : vtkm::cont::ExecutionObjectBase
|
|
{
|
|
public:
|
|
VTKM_CONT
|
|
ConnectivityExplicit() = default;
|
|
|
|
VTKM_CONT
|
|
ConnectivityExplicit(const vtkm::cont::ArrayHandle<vtkm::UInt8>& shapes,
|
|
const vtkm::cont::ArrayHandle<vtkm::IdComponent>& numberOfIndices,
|
|
const vtkm::cont::ArrayHandle<vtkm::Id>& connectivity,
|
|
const vtkm::cont::ArrayHandle<vtkm::Id>& offsets,
|
|
const MIRStats& stats)
|
|
: Shapes(shapes)
|
|
, NumberOfIndices(numberOfIndices)
|
|
, Connectivity(connectivity)
|
|
, Offsets(offsets)
|
|
, Stats(stats)
|
|
{
|
|
}
|
|
|
|
VTKM_CONT ExecutionConnectivityExplicit PrepareForExecution(vtkm::cont::DeviceAdapterId device,
|
|
vtkm::cont::Token& token) const
|
|
{
|
|
ExecutionConnectivityExplicit execConnectivity(this->Shapes,
|
|
this->NumberOfIndices,
|
|
this->Connectivity,
|
|
this->Offsets,
|
|
this->Stats,
|
|
device,
|
|
token);
|
|
return execConnectivity;
|
|
}
|
|
|
|
private:
|
|
vtkm::cont::ArrayHandle<vtkm::UInt8> Shapes;
|
|
vtkm::cont::ArrayHandle<vtkm::IdComponent> NumberOfIndices;
|
|
vtkm::cont::ArrayHandle<vtkm::Id> Connectivity;
|
|
vtkm::cont::ArrayHandle<vtkm::Id> Offsets;
|
|
vtkm::worklet::MIRStats Stats;
|
|
};
|
|
|
|
class ComputeStats : public vtkm::worklet::WorkletVisitCellsWithPoints
|
|
{
|
|
public:
|
|
vtkm::Id targ;
|
|
VTKM_CONT ComputeStats(vtkm::Id target)
|
|
: targ(target)
|
|
{
|
|
}
|
|
VTKM_CONT ComputeStats() = default;
|
|
using ControlSignature = void(CellSetIn,
|
|
WholeArrayIn curVals,
|
|
WholeArrayIn prevVals,
|
|
FieldInCell offsets,
|
|
ExecObject mirTables,
|
|
FieldInCell parentObj,
|
|
FieldInCell prevCol,
|
|
FieldOutCell stats,
|
|
FieldOutCell caseID);
|
|
using ExecutionSignature = void(CellShape, PointCount, _3, _2, _4, _5, _6, _7, _8, _9);
|
|
using InputDomain = _1;
|
|
|
|
template <typename CellShapeTag,
|
|
typename ScalarFieldVec,
|
|
typename ScalarFieldVec1,
|
|
typename DeviceAdapter,
|
|
typename ScalarPos,
|
|
typename ParentObj,
|
|
typename PreCol>
|
|
VTKM_EXEC void operator()(
|
|
const CellShapeTag shape,
|
|
const vtkm::IdComponent pointCount,
|
|
const ScalarFieldVec& prevVals,
|
|
const ScalarFieldVec1& newVals,
|
|
const ScalarPos& valPositionStart,
|
|
const vtkm::worklet::MIRCases::MIRTables::MIRDevicePortal<DeviceAdapter>& MIRData,
|
|
const ParentObj&,
|
|
const PreCol& prevCol,
|
|
MIRStats& MIRStat,
|
|
vtkm::Id& MIRDataIndex) const
|
|
{
|
|
(void)shape;
|
|
vtkm::Id caseId = 0;
|
|
if (prevCol == vtkm::Id(-1))
|
|
{
|
|
// In case of this being the first material for the cell, automatically set it to the furthest case (that is, same shape, color 1)
|
|
for (vtkm::IdComponent iter = pointCount - 1; iter >= 0; iter--)
|
|
{
|
|
caseId++;
|
|
if (iter > 0)
|
|
{
|
|
caseId *= 2;
|
|
}
|
|
}
|
|
}
|
|
else
|
|
{
|
|
for (vtkm::IdComponent iter = pointCount - 1; iter >= 0; iter--)
|
|
{
|
|
if (static_cast<vtkm::Float64>(prevVals[valPositionStart + iter]) <=
|
|
static_cast<vtkm::Float64>(newVals[valPositionStart + iter]))
|
|
{
|
|
caseId++;
|
|
}
|
|
if (iter > 0)
|
|
{
|
|
caseId *= 2;
|
|
}
|
|
}
|
|
}
|
|
// Reinitialize all struct values to 0, experienced weird memory bug otherwise, might be an issue with development environment
|
|
MIRStat.NumberOfCells = 0;
|
|
MIRStat.NumberOfEdgeIndices = 0;
|
|
MIRStat.NumberOfInCellEdgeIndices = 0;
|
|
MIRStat.NumberOfInCellIndices = 0;
|
|
MIRStat.NumberOfInCellInterpPoints = 0;
|
|
MIRStat.NumberOfInCellPoints = 0;
|
|
MIRStat.NumberOfIndices = 0;
|
|
vtkm::Id index = MIRData.GetCaseIndex(shape.Id, caseId, pointCount);
|
|
MIRDataIndex = vtkm::Id(caseId);
|
|
vtkm::Id numberOfCells = MIRData.GetNumberOfShapes(shape.Id, caseId, pointCount);
|
|
if (numberOfCells == -1)
|
|
{
|
|
this->RaiseError("Getting a size index of a polygon with more points than 8 or less points "
|
|
"than 3. Bad case.");
|
|
return;
|
|
}
|
|
MIRStat.NumberOfCells = numberOfCells;
|
|
|
|
for (vtkm::IdComponent shapes = 0; shapes < numberOfCells; shapes++)
|
|
{
|
|
vtkm::UInt8 cellType = MIRData.ValueAt(index++);
|
|
// SH_PNT is a specification that a center point is to be used
|
|
// Note: It is only possible to support 1 midpoint with the current code format
|
|
if (cellType == MIRCases::SH_PNT)
|
|
{
|
|
MIRStat.NumberOfCells = numberOfCells - 1;
|
|
vtkm::UInt8 numberOfIndices = MIRData.ValueAt(index + 2);
|
|
index += 3;
|
|
MIRStat.NumberOfInCellPoints = 1;
|
|
MIRStat.NumberOfInCellInterpPoints = numberOfIndices;
|
|
for (vtkm::IdComponent points = 0; points < numberOfIndices; points++)
|
|
{
|
|
vtkm::Id elem = MIRData.ValueAt(index);
|
|
// If the midpoint needs to reference an edge point, record it.
|
|
MIRStat.NumberOfInCellEdgeIndices += (elem >= MIRCases::EA) ? 1 : 0;
|
|
index++;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
vtkm::Id numberOfIndices = MIRData.GetNumberOfIndices(cellType);
|
|
index++;
|
|
MIRStat.NumberOfIndices += numberOfIndices;
|
|
for (vtkm::IdComponent points = 0; points < numberOfIndices; points++, index++)
|
|
{
|
|
vtkm::IdComponent element = MIRData.ValueAt(index);
|
|
if (element >= MIRCases::EA && element <= MIRCases::EL)
|
|
{
|
|
MIRStat.NumberOfEdgeIndices++;
|
|
}
|
|
else if (element == MIRCases::N0)
|
|
{
|
|
// N0 stands for the midpoint. Technically it could be N0->N3, but with the current
|
|
// setup, only N0 is supported/present in the MIRCases tables.
|
|
MIRStat.NumberOfInCellIndices++;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
};
|
|
class MIRParentObject : public vtkm::cont::ExecutionAndControlObjectBase
|
|
{
|
|
public:
|
|
VTKM_CONT MIRParentObject() = default;
|
|
VTKM_CONT MIRParentObject(vtkm::Id numCells,
|
|
vtkm::cont::ArrayHandle<vtkm::Id> celllook,
|
|
vtkm::cont::ArrayHandle<vtkm::Id> cellCol,
|
|
vtkm::cont::ArrayHandle<vtkm::Id> newCellCol,
|
|
vtkm::cont::ArrayHandle<vtkm::Id> newcellLook)
|
|
: newCellColors(newCellCol)
|
|
, newCellLookback(newcellLook)
|
|
, numberOfInd(numCells)
|
|
, cellLookback(celllook)
|
|
, cellColors(cellCol){};
|
|
|
|
class MIRParentPortal
|
|
{
|
|
public:
|
|
VTKM_EXEC void SetNewCellLookback(vtkm::Id index, vtkm::Id originalIndex)
|
|
{
|
|
this->NewCellLookback.Set(index, originalIndex);
|
|
}
|
|
VTKM_EXEC void SetNewCellColor(vtkm::Id index, vtkm::Id col)
|
|
{
|
|
this->NewCellColors.Set(index, col);
|
|
}
|
|
VTKM_EXEC vtkm::Id GetParentCellIndex(vtkm::Id index) { return this->CellLookback.Get(index); }
|
|
VTKM_EXEC vtkm::Id GetParentCellColor(vtkm::Id index) { return this->CellColors.Get(index); }
|
|
|
|
private:
|
|
typename vtkm::cont::ArrayHandle<vtkm::Id>::ReadPortalType CellLookback;
|
|
typename vtkm::cont::ArrayHandle<vtkm::Id>::ReadPortalType CellColors;
|
|
typename vtkm::cont::ArrayHandle<vtkm::Id>::WritePortalType NewCellColors;
|
|
typename vtkm::cont::ArrayHandle<vtkm::Id>::WritePortalType NewCellLookback;
|
|
friend class MIRParentObject;
|
|
};
|
|
|
|
VTKM_CONT MIRParentPortal PrepareForExecution(vtkm::cont::DeviceAdapterId device,
|
|
vtkm::cont::Token& token)
|
|
{
|
|
MIRParentPortal dev;
|
|
dev.CellLookback = this->cellLookback.PrepareForInput(device, token);
|
|
dev.CellColors = this->cellColors.PrepareForInput(device, token);
|
|
dev.NewCellColors = this->newCellColors.PrepareForOutput(this->numberOfInd, device, token);
|
|
dev.NewCellLookback = this->newCellLookback.PrepareForOutput(this->numberOfInd, device, token);
|
|
return dev;
|
|
}
|
|
vtkm::cont::ArrayHandle<vtkm::Id> newCellColors;
|
|
vtkm::cont::ArrayHandle<vtkm::Id> newCellLookback;
|
|
|
|
private:
|
|
vtkm::Id numberOfInd;
|
|
vtkm::cont::ArrayHandle<vtkm::Id> cellLookback;
|
|
vtkm::cont::ArrayHandle<vtkm::Id> cellColors;
|
|
};
|
|
class GenerateCellSet : public vtkm::worklet::WorkletVisitCellsWithPoints
|
|
{
|
|
public:
|
|
VTKM_EXEC_CONT
|
|
GenerateCellSet(vtkm::Id tar)
|
|
: target(tar)
|
|
{
|
|
}
|
|
|
|
using ControlSignature = void(CellSetIn,
|
|
WholeArrayIn prevVals,
|
|
WholeArrayIn newVals,
|
|
FieldInCell vf_pos,
|
|
FieldInCell mirTableIndices,
|
|
FieldInCell mirStats,
|
|
ExecObject mirTables,
|
|
ExecObject connectivityObject,
|
|
WholeArrayOut edgePointReverseConnectivity,
|
|
WholeArrayOut edgePointInterpolation,
|
|
WholeArrayOut inCellReverseConnectivity,
|
|
WholeArrayOut inCellEdgeReverseConnectivity,
|
|
WholeArrayOut inCellEdgeInterpolation,
|
|
WholeArrayOut inCellInterpolationKeys,
|
|
WholeArrayOut inCellInterpolationInfo,
|
|
ExecObject cellLookbackObj,
|
|
WholeArrayOut simpleLookback);
|
|
|
|
using ExecutionSignature = void(CellShape,
|
|
InputIndex,
|
|
PointCount,
|
|
PointIndices,
|
|
_2,
|
|
_3,
|
|
_4,
|
|
_5,
|
|
_6,
|
|
_7,
|
|
_8,
|
|
_9,
|
|
_10,
|
|
_11,
|
|
_12,
|
|
_13,
|
|
_14,
|
|
_15,
|
|
_16,
|
|
_17); // 20! NO MORE ROOM!
|
|
|
|
template <typename CellShapeTag,
|
|
typename PointVecType,
|
|
typename ScalarVecType1,
|
|
typename ScalarVecType2,
|
|
typename ConnectivityObject,
|
|
typename IdArrayType,
|
|
typename EdgeInterpolationPortalType,
|
|
typename DeviceAdapter,
|
|
typename ScalarPos,
|
|
typename CellLookbackArr>
|
|
VTKM_EXEC void operator()(
|
|
const CellShapeTag shape,
|
|
const vtkm::Id workIndex,
|
|
const vtkm::IdComponent pointcount,
|
|
const PointVecType points,
|
|
const ScalarVecType1& curScalars, // Previous VF
|
|
const ScalarVecType2& newScalars, // New VF
|
|
const ScalarPos& valPositionStart, // Offsets into the ^ arrays for indexing
|
|
const vtkm::Id& clipDataIndex,
|
|
const MIRStats mirStats,
|
|
const worklet::MIRCases::MIRTables::MIRDevicePortal<DeviceAdapter>& MIRData,
|
|
ConnectivityObject& connectivityObject,
|
|
IdArrayType& edgePointReverseConnectivity,
|
|
EdgeInterpolationPortalType& edgePointInterpolation,
|
|
IdArrayType& inCellReverseConnectivity,
|
|
IdArrayType& inCellEdgeReverseConnectivity,
|
|
EdgeInterpolationPortalType& inCellEdgeInterpolation,
|
|
IdArrayType& inCellInterpolationKeys,
|
|
IdArrayType& inCellInterpolationInfo,
|
|
worklet::MIRParentObject::MIRParentPortal& parentObj,
|
|
CellLookbackArr& cellLookbackArray) const
|
|
{
|
|
|
|
(void)shape;
|
|
vtkm::Id clipIndex = MIRData.GetCaseIndex(shape.Id, clipDataIndex, pointcount);
|
|
|
|
// Start index for the cells of this case.
|
|
vtkm::Id cellIndex = mirStats.NumberOfCells;
|
|
// Start index to store connevtivity of this case.
|
|
vtkm::Id connectivityIndex = mirStats.NumberOfIndices;
|
|
// Start indices for reverse mapping into connectivity for this case.
|
|
vtkm::Id edgeIndex = mirStats.NumberOfEdgeIndices;
|
|
vtkm::Id inCellIndex = mirStats.NumberOfInCellIndices;
|
|
vtkm::Id inCellPoints = mirStats.NumberOfInCellPoints;
|
|
// Start Indices to keep track of interpolation points for new cell.
|
|
vtkm::Id inCellInterpPointIndex = mirStats.NumberOfInCellInterpPoints;
|
|
vtkm::Id inCellEdgeInterpIndex = mirStats.NumberOfInCellEdgeIndices;
|
|
|
|
// Iterate over the shapes for the current cell and begin to fill connectivity.
|
|
vtkm::Id numberOfCells = MIRData.GetNumberOfShapes(shape.Id, clipDataIndex, pointcount);
|
|
|
|
for (vtkm::Id cell = 0; cell < numberOfCells; ++cell)
|
|
{
|
|
vtkm::UInt8 cellShape = MIRData.ValueAt(clipIndex++);
|
|
if (cellShape == MIRCases::SH_PNT)
|
|
{
|
|
clipIndex += 2;
|
|
vtkm::IdComponent numberOfPoints = MIRData.ValueAt(clipIndex);
|
|
clipIndex++;
|
|
// Case for a new cell point
|
|
|
|
// 1. Output the input cell id for which we need to generate new point.
|
|
// 2. Output number of points used for interpolation.
|
|
// 3. If vertex
|
|
// - Add vertex to connectivity interpolation information.
|
|
// 4. If edge
|
|
// - Add edge interpolation information for new points.
|
|
// - Reverse connectivity map for new points.
|
|
// Make an array which has all the elements that need to be used
|
|
// for interpolation.
|
|
for (vtkm::IdComponent point = 0; point < numberOfPoints;
|
|
point++, inCellInterpPointIndex++, clipIndex++)
|
|
{
|
|
vtkm::IdComponent entry = static_cast<vtkm::IdComponent>(MIRData.ValueAt(clipIndex));
|
|
inCellInterpolationKeys.Set(inCellInterpPointIndex, workIndex);
|
|
if (entry <= MIRCases::P7)
|
|
{
|
|
inCellInterpolationInfo.Set(inCellInterpPointIndex, points[entry]);
|
|
}
|
|
else
|
|
{
|
|
internal::ClipTables::EdgeVec edge =
|
|
MIRData.GetEdge(shape.Id, entry - MIRCases::EA, pointcount);
|
|
if (edge[0] == 255 || edge[1] == 255)
|
|
{
|
|
this->RaiseError("Edge vertices are assigned incorrect values.");
|
|
return;
|
|
}
|
|
|
|
EdgeInterpolation ei;
|
|
ei.Vertex1 = points[edge[0]];
|
|
ei.Vertex2 = points[edge[1]];
|
|
// For consistency purposes keep the points ordered.
|
|
if (ei.Vertex1 > ei.Vertex2)
|
|
{
|
|
this->swap(ei.Vertex1, ei.Vertex2);
|
|
this->swap(edge[0], edge[1]);
|
|
}
|
|
// need to swap the weight of the point to be A-C / ((D-C) - (B-A)),
|
|
// where A and C are edge0 mats 1 and 2, and B and D are edge1 mats 1 and 2.
|
|
ei.Weight = vtkm::Float64(1) +
|
|
((static_cast<vtkm::Float64>(curScalars[valPositionStart + edge[0]] -
|
|
newScalars[valPositionStart + edge[0]])) /
|
|
static_cast<vtkm::Float64>(
|
|
curScalars[valPositionStart + edge[1]] - curScalars[valPositionStart + edge[0]] +
|
|
newScalars[valPositionStart + edge[0]] - newScalars[valPositionStart + edge[1]]));
|
|
|
|
inCellEdgeReverseConnectivity.Set(inCellEdgeInterpIndex, inCellInterpPointIndex);
|
|
inCellEdgeInterpolation.Set(inCellEdgeInterpIndex, ei);
|
|
inCellEdgeInterpIndex++;
|
|
}
|
|
}
|
|
}
|
|
else
|
|
{
|
|
vtkm::IdComponent numberOfPoints =
|
|
static_cast<vtkm::IdComponent>(MIRData.GetNumberOfIndices(cellShape));
|
|
vtkm::IdComponent colorQ = static_cast<vtkm::IdComponent>(MIRData.ValueAt(clipIndex++));
|
|
vtkm::Id color = colorQ == vtkm::IdComponent(MIRCases::COLOR0)
|
|
? parentObj.GetParentCellColor(workIndex)
|
|
: target;
|
|
parentObj.SetNewCellColor(cellIndex, color);
|
|
parentObj.SetNewCellLookback(cellIndex, parentObj.GetParentCellIndex(workIndex));
|
|
connectivityObject.SetCellShape(cellIndex, cellShape);
|
|
connectivityObject.SetNumberOfIndices(cellIndex, numberOfPoints);
|
|
connectivityObject.SetIndexOffset(cellIndex, connectivityIndex);
|
|
|
|
for (vtkm::IdComponent point = 0; point < numberOfPoints; point++, clipIndex++)
|
|
{
|
|
vtkm::IdComponent entry = static_cast<vtkm::IdComponent>(MIRData.ValueAt(clipIndex));
|
|
if (entry == MIRCases::N0) // case of cell point interpolation
|
|
{
|
|
// Add index of the corresponding cell point.
|
|
inCellReverseConnectivity.Set(inCellIndex++, connectivityIndex);
|
|
connectivityObject.SetConnectivity(connectivityIndex, inCellPoints);
|
|
connectivityIndex++;
|
|
}
|
|
else if (entry <= MIRCases::P7) // existing vertex
|
|
{
|
|
connectivityObject.SetConnectivity(connectivityIndex, points[entry]);
|
|
connectivityIndex++;
|
|
}
|
|
else // case of a new edge point
|
|
{
|
|
internal::ClipTables::EdgeVec edge =
|
|
MIRData.GetEdge(shape.Id, entry - MIRCases::EA, pointcount);
|
|
if (edge[0] == 255 || edge[1] == 255)
|
|
{
|
|
this->RaiseError("Edge vertices are assigned incorrect values.");
|
|
return;
|
|
}
|
|
EdgeInterpolation ei;
|
|
ei.Vertex1 = points[edge[0]];
|
|
ei.Vertex2 = points[edge[1]];
|
|
// For consistency purposes keep the points ordered.
|
|
if (ei.Vertex1 > ei.Vertex2)
|
|
{
|
|
this->swap(ei.Vertex1, ei.Vertex2);
|
|
this->swap(edge[0], edge[1]);
|
|
}
|
|
|
|
ei.Weight = vtkm::Float64(1) +
|
|
((static_cast<vtkm::Float64>(curScalars[valPositionStart + edge[0]] -
|
|
newScalars[valPositionStart + edge[0]])) /
|
|
static_cast<vtkm::Float64>(
|
|
curScalars[valPositionStart + edge[1]] - curScalars[valPositionStart + edge[0]] +
|
|
newScalars[valPositionStart + edge[0]] - newScalars[valPositionStart + edge[1]]));
|
|
//Add to set of new edge points
|
|
//Add reverse connectivity;
|
|
edgePointReverseConnectivity.Set(edgeIndex, connectivityIndex++);
|
|
edgePointInterpolation.Set(edgeIndex, ei);
|
|
edgeIndex++;
|
|
}
|
|
}
|
|
// Set cell matID...
|
|
cellLookbackArray.Set(cellIndex, workIndex);
|
|
++cellIndex;
|
|
}
|
|
}
|
|
}
|
|
|
|
template <typename T>
|
|
VTKM_EXEC void swap(T& v1, T& v2) const
|
|
{
|
|
T temp = v1;
|
|
v1 = v2;
|
|
v2 = temp;
|
|
}
|
|
|
|
private:
|
|
vtkm::Id target;
|
|
};
|
|
class ScatterEdgeConnectivity : public vtkm::worklet::WorkletMapField
|
|
{
|
|
public:
|
|
VTKM_CONT
|
|
ScatterEdgeConnectivity(vtkm::Id edgePointOffset)
|
|
: EdgePointOffset(edgePointOffset)
|
|
{
|
|
}
|
|
|
|
using ControlSignature = void(FieldIn sourceValue,
|
|
FieldIn destinationIndices,
|
|
WholeArrayOut destinationData);
|
|
|
|
using ExecutionSignature = void(_1, _2, _3);
|
|
|
|
using InputDomain = _1;
|
|
|
|
template <typename ConnectivityDataType>
|
|
VTKM_EXEC void operator()(const vtkm::Id sourceValue,
|
|
const vtkm::Id destinationIndex,
|
|
ConnectivityDataType& destinationData) const
|
|
{
|
|
destinationData.Set(destinationIndex, (sourceValue + EdgePointOffset));
|
|
}
|
|
|
|
private:
|
|
vtkm::Id EdgePointOffset;
|
|
};
|
|
class ScatterInCellConnectivity : public vtkm::worklet::WorkletMapField
|
|
{
|
|
public:
|
|
VTKM_CONT
|
|
ScatterInCellConnectivity(vtkm::Id inCellPointOffset)
|
|
: InCellPointOffset(inCellPointOffset)
|
|
{
|
|
}
|
|
|
|
using ControlSignature = void(FieldIn destinationIndices, WholeArrayOut destinationData);
|
|
|
|
using ExecutionSignature = void(_1, _2);
|
|
|
|
using InputDomain = _1;
|
|
|
|
template <typename ConnectivityDataType>
|
|
VTKM_EXEC void operator()(const vtkm::Id destinationIndex,
|
|
ConnectivityDataType& destinationData) const
|
|
{
|
|
auto sourceValue = destinationData.Get(destinationIndex);
|
|
destinationData.Set(destinationIndex, (sourceValue + InCellPointOffset));
|
|
}
|
|
|
|
private:
|
|
vtkm::Id InCellPointOffset;
|
|
};
|
|
class MIR
|
|
{
|
|
public:
|
|
MIR()
|
|
: MIRTablesInstance()
|
|
, EdgePointsInterpolation()
|
|
, InCellInterpolationKeys()
|
|
, InCellInterpolationInfo()
|
|
, CellMapOutputToInput()
|
|
, EdgePointsOffset()
|
|
, InCellPointsOffset()
|
|
{
|
|
}
|
|
template <typename VFList1, typename VFList2, typename CellSet, typename VFLocs, typename IDList>
|
|
vtkm::cont::CellSetExplicit<> Run(const CellSet& cellSet,
|
|
const VFList1& prevValues,
|
|
const VFList2& curValues,
|
|
const VFLocs& offsets,
|
|
const IDList& prevIDs,
|
|
const vtkm::Id& newID,
|
|
const IDList& prevLookback,
|
|
IDList& newIDs,
|
|
IDList& newLookback)
|
|
{
|
|
// First compute the stats for the MIR algorithm & build the offsets
|
|
//{
|
|
ComputeStats statWorklet(newID);
|
|
vtkm::worklet::DispatcherMapTopology<ComputeStats> statsDispatch(statWorklet);
|
|
|
|
// Output variables
|
|
vtkm::cont::ArrayHandle<MIRStats> mirStats;
|
|
vtkm::cont::ArrayHandle<vtkm::Id> mirInd;
|
|
|
|
statsDispatch.Invoke(cellSet,
|
|
curValues,
|
|
prevValues,
|
|
offsets,
|
|
this->MIRTablesInstance,
|
|
prevLookback,
|
|
prevIDs,
|
|
mirStats,
|
|
mirInd);
|
|
// Sum all stats to form an offset array (for indexing in the MIR algorithm)
|
|
MIRStats zero;
|
|
vtkm::cont::ArrayHandle<MIRStats> cellSetStats;
|
|
MIRStats total =
|
|
vtkm::cont::Algorithm::ScanExclusive(mirStats, cellSetStats, MIRStats::SumOp(), zero);
|
|
mirStats.ReleaseResources();
|
|
//}
|
|
// Secondly, build the sets.
|
|
//{
|
|
// CellSetExplicit sets
|
|
vtkm::cont::ArrayHandle<vtkm::UInt8> shapes;
|
|
vtkm::cont::ArrayHandle<vtkm::IdComponent> numberOfIndices;
|
|
vtkm::cont::ArrayHandle<vtkm::Id> connectivity;
|
|
vtkm::cont::ArrayHandle<vtkm::Id> offset;
|
|
ConnectivityExplicit connectivityObject(shapes, numberOfIndices, connectivity, offset, total);
|
|
// Connectivity related sets
|
|
vtkm::cont::ArrayHandle<vtkm::Id> edgePointReverseConnectivity;
|
|
edgePointReverseConnectivity.Allocate(total.NumberOfEdgeIndices);
|
|
vtkm::cont::ArrayHandle<EdgeInterpolation> edgeInterpolation;
|
|
edgeInterpolation.Allocate(total.NumberOfEdgeIndices);
|
|
vtkm::cont::ArrayHandle<vtkm::Id> cellPointReverseConnectivity;
|
|
cellPointReverseConnectivity.Allocate(total.NumberOfInCellIndices);
|
|
vtkm::cont::ArrayHandle<vtkm::Id> cellPointEdgeReverseConnectivity;
|
|
cellPointEdgeReverseConnectivity.Allocate(total.NumberOfInCellEdgeIndices);
|
|
vtkm::cont::ArrayHandle<EdgeInterpolation> cellPointEdgeInterpolation;
|
|
cellPointEdgeInterpolation.Allocate(total.NumberOfInCellEdgeIndices);
|
|
this->InCellInterpolationKeys.Allocate(total.NumberOfInCellInterpPoints);
|
|
this->InCellInterpolationInfo.Allocate(total.NumberOfInCellInterpPoints);
|
|
this->CellMapOutputToInput.Allocate(total.NumberOfCells);
|
|
|
|
|
|
//}
|
|
// Thirdly, call the MIR generator
|
|
//{
|
|
GenerateCellSet cellSetWorklet(newID);
|
|
vtkm::worklet::DispatcherMapTopology<GenerateCellSet> cellSetDispatcher(cellSetWorklet);
|
|
// Output arrays storing information about cell lookbacks and cell material IDs
|
|
vtkm::cont::ArrayHandle<vtkm::Id> nextID, nextLookback;
|
|
nextID.Allocate(total.NumberOfCells);
|
|
nextLookback.Allocate(total.NumberOfCells);
|
|
MIRParentObject po(total.NumberOfCells, prevLookback, prevIDs, nextID, nextLookback);
|
|
|
|
|
|
// Perform the MIR step
|
|
cellSetDispatcher.Invoke(cellSet,
|
|
prevValues,
|
|
curValues,
|
|
offsets,
|
|
mirInd,
|
|
cellSetStats,
|
|
this->MIRTablesInstance,
|
|
connectivityObject,
|
|
edgePointReverseConnectivity,
|
|
edgeInterpolation,
|
|
cellPointReverseConnectivity,
|
|
cellPointEdgeReverseConnectivity,
|
|
cellPointEdgeInterpolation,
|
|
this->InCellInterpolationKeys,
|
|
this->InCellInterpolationInfo,
|
|
po,
|
|
this->CellMapOutputToInput);
|
|
|
|
//}
|
|
// Forthly, create the output set and clean up connectivity objects.
|
|
//{
|
|
// Get unique keys for all shared edges
|
|
vtkm::cont::Algorithm::SortByKey(
|
|
edgeInterpolation, edgePointReverseConnectivity, EdgeInterpolation::LessThanOp());
|
|
vtkm::cont::Algorithm::Copy(edgeInterpolation, this->EdgePointsInterpolation);
|
|
vtkm::cont::Algorithm::Unique(this->EdgePointsInterpolation, EdgeInterpolation::EqualToOp());
|
|
vtkm::cont::ArrayHandle<vtkm::Id> edgeInterpolationIndexToUnique;
|
|
vtkm::cont::Algorithm::LowerBounds(this->EdgePointsInterpolation,
|
|
edgeInterpolation,
|
|
edgeInterpolationIndexToUnique,
|
|
EdgeInterpolation::LessThanOp());
|
|
|
|
vtkm::cont::ArrayHandle<vtkm::Id> cellInterpolationIndexToUnique;
|
|
vtkm::cont::Algorithm::LowerBounds(this->EdgePointsInterpolation,
|
|
cellPointEdgeInterpolation,
|
|
cellInterpolationIndexToUnique,
|
|
EdgeInterpolation::LessThanOp());
|
|
this->EdgePointsOffset = cellSet.GetNumberOfPoints();
|
|
this->InCellPointsOffset =
|
|
this->EdgePointsOffset + this->EdgePointsInterpolation.GetNumberOfValues();
|
|
|
|
ScatterEdgeConnectivity scatterEdgePointConnectivity(this->EdgePointsOffset);
|
|
vtkm::worklet::DispatcherMapField<ScatterEdgeConnectivity> scatterEdgeDispatcher(
|
|
scatterEdgePointConnectivity);
|
|
scatterEdgeDispatcher.Invoke(
|
|
edgeInterpolationIndexToUnique, edgePointReverseConnectivity, connectivity);
|
|
scatterEdgeDispatcher.Invoke(cellInterpolationIndexToUnique,
|
|
cellPointEdgeReverseConnectivity,
|
|
this->InCellInterpolationInfo);
|
|
// Add offset in connectivity of all new in-cell points.
|
|
ScatterInCellConnectivity scatterInCellPointConnectivity(this->InCellPointsOffset);
|
|
vtkm::worklet::DispatcherMapField<ScatterInCellConnectivity> scatterInCellDispatcher(
|
|
scatterInCellPointConnectivity);
|
|
scatterInCellDispatcher.Invoke(cellPointReverseConnectivity, connectivity);
|
|
|
|
vtkm::cont::CellSetExplicit<> output;
|
|
vtkm::Id numberOfPoints = cellSet.GetNumberOfPoints() +
|
|
this->EdgePointsInterpolation.GetNumberOfValues() + total.NumberOfInCellPoints;
|
|
|
|
vtkm::cont::ConvertNumComponentsToOffsets(numberOfIndices, offset);
|
|
// Create explicit cell set output
|
|
output.Fill(numberOfPoints, shapes, connectivity, offset);
|
|
//}
|
|
vtkm::cont::ArrayCopy(po.newCellColors, newIDs);
|
|
vtkm::cont::ArrayCopy(po.newCellLookback, newLookback);
|
|
|
|
return output;
|
|
}
|
|
|
|
template <typename ArrayHandleType>
|
|
class InterpolateField
|
|
{
|
|
public:
|
|
using ValueType = typename ArrayHandleType::ValueType;
|
|
using TypeMappedValue = vtkm::List<ValueType>;
|
|
|
|
InterpolateField(vtkm::cont::ArrayHandle<EdgeInterpolation> edgeInterpolationArray,
|
|
vtkm::cont::ArrayHandle<vtkm::Id> inCellInterpolationKeys,
|
|
vtkm::cont::ArrayHandle<vtkm::Id> inCellInterpolationInfo,
|
|
vtkm::Id edgePointsOffset,
|
|
vtkm::Id inCellPointsOffset,
|
|
ArrayHandleType* output)
|
|
: EdgeInterpolationArray(edgeInterpolationArray)
|
|
, InCellInterpolationKeys(inCellInterpolationKeys)
|
|
, InCellInterpolationInfo(inCellInterpolationInfo)
|
|
, EdgePointsOffset(edgePointsOffset)
|
|
, InCellPointsOffset(inCellPointsOffset)
|
|
, Output(output)
|
|
{
|
|
}
|
|
|
|
class PerformEdgeInterpolations : public vtkm::worklet::WorkletMapField
|
|
{
|
|
public:
|
|
PerformEdgeInterpolations(vtkm::Id edgePointsOffset)
|
|
: EdgePointsOffset(edgePointsOffset)
|
|
{
|
|
}
|
|
|
|
using ControlSignature = void(FieldIn edgeInterpolations, WholeArrayInOut outputField);
|
|
|
|
using ExecutionSignature = void(_1, _2, WorkIndex);
|
|
|
|
template <typename EdgeInterp, typename OutputFieldPortal>
|
|
VTKM_EXEC void operator()(const EdgeInterp& ei,
|
|
OutputFieldPortal& field,
|
|
const vtkm::Id workIndex) const
|
|
{
|
|
using T = typename OutputFieldPortal::ValueType;
|
|
T v1 = field.Get(ei.Vertex1);
|
|
T v2 = field.Get(ei.Vertex2);
|
|
field.Set(this->EdgePointsOffset + workIndex,
|
|
static_cast<T>(MIRinternal::Scale(T(v1 - v2), ei.Weight) + v2));
|
|
if (ei.Weight > vtkm::Float64(1) || ei.Weight < vtkm::Float64(0))
|
|
{
|
|
this->RaiseError("Error in edge weight, assigned value not it interval [0,1].");
|
|
}
|
|
}
|
|
|
|
private:
|
|
vtkm::Id EdgePointsOffset;
|
|
};
|
|
|
|
class PerformInCellInterpolations : public vtkm::worklet::WorkletReduceByKey
|
|
{
|
|
public:
|
|
using ControlSignature = void(KeysIn keys, ValuesIn toReduce, ReducedValuesOut centroid);
|
|
|
|
using ExecutionSignature = void(_2, _3);
|
|
|
|
template <typename MappedValueVecType, typename MappedValueType>
|
|
VTKM_EXEC void operator()(const MappedValueVecType& toReduce, MappedValueType& centroid) const
|
|
{
|
|
vtkm::IdComponent numValues = toReduce.GetNumberOfComponents();
|
|
MappedValueType sum = toReduce[0];
|
|
for (vtkm::IdComponent i = 1; i < numValues; i++)
|
|
{
|
|
MappedValueType value = toReduce[i];
|
|
// static_cast is for when MappedValueType is a small int that gets promoted to int32.
|
|
sum = static_cast<MappedValueType>(sum + value);
|
|
}
|
|
centroid = MIRinternal::Scale(sum, 1. / static_cast<vtkm::Float64>(numValues));
|
|
}
|
|
};
|
|
|
|
template <typename Storage>
|
|
VTKM_CONT void operator()(const vtkm::cont::ArrayHandle<ValueType, Storage>& field) const
|
|
{
|
|
vtkm::worklet::Keys<vtkm::Id> interpolationKeys(InCellInterpolationKeys);
|
|
|
|
vtkm::Id numberOfOriginalValues = field.GetNumberOfValues();
|
|
vtkm::Id numberOfEdgePoints = EdgeInterpolationArray.GetNumberOfValues();
|
|
vtkm::Id numberOfInCellPoints = interpolationKeys.GetUniqueKeys().GetNumberOfValues();
|
|
|
|
ArrayHandleType result;
|
|
result.Allocate(numberOfOriginalValues + numberOfEdgePoints + numberOfInCellPoints);
|
|
vtkm::cont::Algorithm::CopySubRange(field, 0, numberOfOriginalValues, result);
|
|
|
|
PerformEdgeInterpolations edgeInterpWorklet(numberOfOriginalValues);
|
|
vtkm::worklet::DispatcherMapField<PerformEdgeInterpolations> edgeInterpDispatcher(
|
|
edgeInterpWorklet);
|
|
edgeInterpDispatcher.Invoke(this->EdgeInterpolationArray, result);
|
|
|
|
// Perform a gather on output to get all required values for calculation of
|
|
// centroids using the interpolation info array.
|
|
using IdHandle = vtkm::cont::ArrayHandle<vtkm::Id>;
|
|
using ValueHandle = vtkm::cont::ArrayHandle<ValueType>;
|
|
vtkm::cont::ArrayHandlePermutation<IdHandle, ValueHandle> toReduceValues(
|
|
InCellInterpolationInfo, result);
|
|
|
|
vtkm::cont::ArrayHandle<ValueType> reducedValues;
|
|
vtkm::worklet::DispatcherReduceByKey<PerformInCellInterpolations>
|
|
inCellInterpolationDispatcher;
|
|
inCellInterpolationDispatcher.Invoke(interpolationKeys, toReduceValues, reducedValues);
|
|
vtkm::Id inCellPointsOffset = numberOfOriginalValues + numberOfEdgePoints;
|
|
vtkm::cont::Algorithm::CopySubRange(
|
|
reducedValues, 0, reducedValues.GetNumberOfValues(), result, inCellPointsOffset);
|
|
*(this->Output) = result;
|
|
}
|
|
|
|
private:
|
|
vtkm::cont::ArrayHandle<EdgeInterpolation> EdgeInterpolationArray;
|
|
vtkm::cont::ArrayHandle<vtkm::Id> InCellInterpolationKeys;
|
|
vtkm::cont::ArrayHandle<vtkm::Id> InCellInterpolationInfo;
|
|
vtkm::Id EdgePointsOffset;
|
|
vtkm::Id InCellPointsOffset;
|
|
ArrayHandleType* Output;
|
|
};
|
|
|
|
template <typename IDLen, typename IDPos, typename IDList, typename VFList>
|
|
class InterpolateMIRFields
|
|
{
|
|
public:
|
|
InterpolateMIRFields(vtkm::cont::ArrayHandle<EdgeInterpolation> edgeInterpolationArray,
|
|
vtkm::cont::ArrayHandle<vtkm::Id> inCellInterpolationKeys,
|
|
vtkm::cont::ArrayHandle<vtkm::Id> inCellInterpolationInfo,
|
|
vtkm::Id edgePointsOffset,
|
|
vtkm::Id inCellPointsOffset,
|
|
IDLen* output1,
|
|
IDPos* output2,
|
|
IDList* output3,
|
|
VFList* output4)
|
|
: EdgeInterpolationArray(edgeInterpolationArray)
|
|
, InCellInterpolationKeys(inCellInterpolationKeys)
|
|
, InCellInterpolationInfo(inCellInterpolationInfo)
|
|
, EdgePointsOffset(edgePointsOffset)
|
|
, InCellPointsOffset(inCellPointsOffset)
|
|
, LenOut(output1)
|
|
, PosOut(output2)
|
|
, IDOut(output3)
|
|
, VFOut(output4)
|
|
{
|
|
}
|
|
|
|
class PerformEdgeInterpolations : public vtkm::worklet::WorkletMapField
|
|
{
|
|
public:
|
|
PerformEdgeInterpolations(vtkm::Id edgePointsOffset)
|
|
: EdgePointsOffset(edgePointsOffset)
|
|
{
|
|
}
|
|
|
|
using ControlSignature = void(FieldIn edgeInterpolations,
|
|
WholeArrayIn lengths,
|
|
WholeArrayIn positions,
|
|
WholeArrayInOut ids,
|
|
WholeArrayInOut vfs);
|
|
|
|
using ExecutionSignature = void(_1, _2, _3, _4, _5, WorkIndex);
|
|
|
|
template <typename EdgeInterp, typename IDL, typename IDO, typename IdsVec, typename VfsVec>
|
|
VTKM_EXEC void operator()(const EdgeInterp& ei,
|
|
const IDL& lengths,
|
|
const IDO& positions,
|
|
IdsVec& ids,
|
|
VfsVec& vfs,
|
|
const vtkm::Id workIndex) const
|
|
{
|
|
vtkm::Vec<vtkm::Id, 2> idOff;
|
|
vtkm::Vec<vtkm::Id, 2> idLen;
|
|
vtkm::Vec<vtkm::Id, 2> idInd;
|
|
vtkm::Vec<vtkm::Float64, 2> multiplier;
|
|
multiplier[1] = vtkm::Float64(1.0) - ei.Weight;
|
|
multiplier[0] = ei.Weight;
|
|
vtkm::Id uniqueMats = vtkm::Id(0);
|
|
|
|
idOff[0] = vtkm::Id(0);
|
|
idOff[1] = idOff[0];
|
|
idInd[0] = positions.Get(ei.Vertex1);
|
|
idInd[1] = positions.Get(ei.Vertex2);
|
|
idLen[0] = lengths.Get(ei.Vertex1);
|
|
idLen[1] = lengths.Get(ei.Vertex2);
|
|
vtkm::IdComponent numberOfPoints = 2;
|
|
vtkm::UInt8 hasWork = vtkm::UInt8(1);
|
|
while (hasWork != vtkm::UInt8(0))
|
|
{
|
|
hasWork = vtkm::UInt8(0);
|
|
vtkm::Id lowest = vtkm::Id(-1);
|
|
for (vtkm::IdComponent i = 0; i < numberOfPoints; i++)
|
|
{
|
|
if (idOff[i] < idLen[i])
|
|
{
|
|
vtkm::Id tmp = ids.Get(idInd[i] + idOff[i]);
|
|
if (lowest == vtkm::Id(-1) || tmp < lowest)
|
|
{
|
|
lowest = tmp;
|
|
hasWork = vtkm::UInt8(1);
|
|
}
|
|
}
|
|
}
|
|
if (hasWork != vtkm::UInt8(0))
|
|
{
|
|
vtkm::Float64 vfVal = vtkm::Float64(0);
|
|
for (vtkm::IdComponent i = 0; i < numberOfPoints; i++)
|
|
{
|
|
if (idOff[i] < idLen[i])
|
|
{
|
|
vtkm::Id tmp = ids.Get(idInd[i] + idOff[i]);
|
|
if (lowest == tmp)
|
|
{
|
|
vfVal += multiplier[i] * vfs.Get(idInd[i] + idOff[i]);
|
|
idOff[i]++;
|
|
}
|
|
}
|
|
}
|
|
ids.Set(positions.Get(this->EdgePointsOffset + workIndex) + uniqueMats, lowest);
|
|
vfs.Set(positions.Get(this->EdgePointsOffset + workIndex) + uniqueMats, vfVal);
|
|
uniqueMats++;
|
|
}
|
|
}
|
|
}
|
|
|
|
private:
|
|
vtkm::Id EdgePointsOffset;
|
|
};
|
|
class PerformEdgeInterpolations_C : public vtkm::worklet::WorkletMapField
|
|
{
|
|
private:
|
|
vtkm::Id EdgePointsOffset;
|
|
|
|
public:
|
|
PerformEdgeInterpolations_C(vtkm::Id edgePointsOffset)
|
|
: EdgePointsOffset(edgePointsOffset)
|
|
{
|
|
}
|
|
using ControlSignature = void(FieldIn edgeInterpolations,
|
|
WholeArrayInOut IDLengths,
|
|
WholeArrayIn IDOffsets,
|
|
WholeArrayIn IDs,
|
|
FieldOut edgeLength);
|
|
using ExecutionSignature = void(_1, _2, _3, _4, WorkIndex, _5);
|
|
template <typename EdgeInterp, typename IDL, typename IDO, typename IdsVec, typename ELL>
|
|
VTKM_EXEC void operator()(const EdgeInterp& ei,
|
|
IDL& lengths,
|
|
const IDO& positions,
|
|
const IdsVec& ids,
|
|
const vtkm::Id workIndex,
|
|
ELL& edgelength) const
|
|
{
|
|
vtkm::Vec<vtkm::Id, 2> idOff;
|
|
vtkm::Vec<vtkm::Id, 2> idLen;
|
|
vtkm::Vec<vtkm::Id, 2> idInd;
|
|
vtkm::Id uniqueMats = vtkm::Id(0);
|
|
|
|
idOff[0] = vtkm::Id(0);
|
|
idOff[1] = idOff[0];
|
|
idInd[0] = positions.Get(ei.Vertex1);
|
|
idInd[1] = positions.Get(ei.Vertex2);
|
|
idLen[0] = lengths.Get(ei.Vertex1);
|
|
idLen[1] = lengths.Get(ei.Vertex2);
|
|
vtkm::IdComponent numberOfPoints = 2;
|
|
vtkm::UInt8 hasWork = vtkm::UInt8(1);
|
|
while (hasWork != vtkm::UInt8(0))
|
|
{
|
|
hasWork = vtkm::UInt8(0);
|
|
vtkm::Id lowest = vtkm::Id(-1);
|
|
for (vtkm::IdComponent i = 0; i < numberOfPoints; i++)
|
|
{
|
|
if (idOff[i] < idLen[i])
|
|
{
|
|
vtkm::Id tmp = ids.Get(idInd[i] + idOff[i]);
|
|
if (lowest == vtkm::Id(-1) || tmp < lowest)
|
|
{
|
|
lowest = tmp;
|
|
hasWork = vtkm::UInt8(1);
|
|
}
|
|
}
|
|
}
|
|
if (hasWork != vtkm::UInt8(0))
|
|
{
|
|
for (vtkm::IdComponent i = 0; i < numberOfPoints; i++)
|
|
{
|
|
if (idOff[i] < idLen[i])
|
|
{
|
|
vtkm::Id tmp = ids.Get(idInd[i] + idOff[i]);
|
|
if (lowest == tmp)
|
|
{
|
|
idOff[i]++;
|
|
}
|
|
}
|
|
}
|
|
uniqueMats++;
|
|
}
|
|
}
|
|
lengths.Set(this->EdgePointsOffset + workIndex, uniqueMats);
|
|
edgelength = uniqueMats;
|
|
}
|
|
};
|
|
|
|
class PerformInCellInterpolations_C : public vtkm::worklet::WorkletReduceByKey
|
|
{
|
|
public:
|
|
using ControlSignature = void(KeysIn keys,
|
|
ValuesIn toReduce,
|
|
WholeArrayIn IDLengths,
|
|
WholeArrayIn IDOffsets,
|
|
WholeArrayIn IDs,
|
|
ReducedValuesOut centroid);
|
|
|
|
using ExecutionSignature = void(_2, _3, _4, _5, _6);
|
|
|
|
template <typename MappedValueVecType,
|
|
typename MappedValueType,
|
|
typename IDArr,
|
|
typename IDOff,
|
|
typename IdsVec>
|
|
VTKM_EXEC void operator()(const MappedValueVecType& toReduce,
|
|
const IDArr& lengths,
|
|
const IDOff& positions,
|
|
const IdsVec& ids,
|
|
MappedValueType& numIdNeeded) const
|
|
{
|
|
vtkm::IdComponent numberOfPoints = toReduce.GetNumberOfComponents();
|
|
// ToReduce is simply the indexArray, giving us point information (since this is reduce by key)
|
|
// numIdNeeded is the output length of this key
|
|
using IdVec = vtkm::Vec<vtkm::Id, 8>;
|
|
IdVec idOff = vtkm::TypeTraits<IdVec>::ZeroInitialization();
|
|
IdVec idLen = vtkm::TypeTraits<IdVec>::ZeroInitialization();
|
|
IdVec idInd = vtkm::TypeTraits<IdVec>::ZeroInitialization();
|
|
vtkm::Id uniqueMats = vtkm::Id(0);
|
|
|
|
for (vtkm::IdComponent i = 0; i < numberOfPoints; i++)
|
|
{
|
|
idOff[i] = 0;
|
|
idLen[i] = lengths.Get(toReduce[i]);
|
|
idInd[i] = positions.Get(toReduce[i]);
|
|
}
|
|
|
|
vtkm::UInt8 hasWork = vtkm::UInt8(1);
|
|
while (hasWork != vtkm::UInt8(0))
|
|
{
|
|
hasWork = vtkm::UInt8(0);
|
|
vtkm::Id lowest = vtkm::Id(-1);
|
|
for (vtkm::IdComponent i = 0; i < numberOfPoints; i++)
|
|
{
|
|
if (idOff[i] < idLen[i])
|
|
{
|
|
vtkm::Id tmp = ids.Get(idInd[i] + idOff[i]);
|
|
if (lowest == vtkm::Id(-1) || tmp < lowest)
|
|
{
|
|
lowest = tmp;
|
|
hasWork = vtkm::UInt8(1);
|
|
}
|
|
}
|
|
}
|
|
if (hasWork != vtkm::UInt8(0))
|
|
{
|
|
for (vtkm::IdComponent i = 0; i < numberOfPoints; i++)
|
|
{
|
|
if (idOff[i] < idLen[i])
|
|
{
|
|
vtkm::Id tmp = ids.Get(idInd[i] + idOff[i]);
|
|
if (lowest == tmp)
|
|
{
|
|
idOff[i]++;
|
|
}
|
|
}
|
|
}
|
|
uniqueMats++;
|
|
}
|
|
}
|
|
numIdNeeded = uniqueMats;
|
|
}
|
|
};
|
|
|
|
class PerformInCellInterpolations : public vtkm::worklet::WorkletReduceByKey
|
|
{
|
|
private:
|
|
vtkm::Id offset;
|
|
|
|
public:
|
|
PerformInCellInterpolations(vtkm::Id outputOffsetForBookkeeping)
|
|
: offset(outputOffsetForBookkeeping)
|
|
{
|
|
}
|
|
using ControlSignature = void(KeysIn keys,
|
|
ValuesIn toReduce,
|
|
WholeArrayIn IDLengths,
|
|
WholeArrayIn IDOffsets,
|
|
WholeArrayIn IDs,
|
|
WholeArrayIn VFs,
|
|
ReducedValuesIn indexOff,
|
|
ReducedValuesOut reindexedOut,
|
|
WholeArrayOut outputIDs,
|
|
WholeArrayOut outputVFs);
|
|
|
|
using ExecutionSignature = void(_2, _3, _4, _5, _6, _7, _8, _9, _10);
|
|
|
|
template <typename MappedValueVecType,
|
|
typename IDArr,
|
|
typename IDOff,
|
|
typename IdsVec,
|
|
typename VfsVec,
|
|
typename IndexIn,
|
|
typename IndexOut,
|
|
typename OutID,
|
|
typename OutVF>
|
|
VTKM_EXEC void operator()(const MappedValueVecType& toReduce,
|
|
const IDArr& lengths,
|
|
const IDOff& positions,
|
|
const IdsVec& ids,
|
|
const VfsVec& vfs,
|
|
const IndexIn& localOffset,
|
|
IndexOut& globalOffset,
|
|
OutID& outIDs,
|
|
OutVF& outVFs) const
|
|
{
|
|
|
|
globalOffset = localOffset + this->offset;
|
|
vtkm::IdComponent numberOfPoints = toReduce.GetNumberOfComponents();
|
|
// ToReduce is simply the indexArray, giving us point information (since this is reduce by key)
|
|
|
|
// numIdNeeded is the output length of this key
|
|
using IdVec = vtkm::Vec<vtkm::Id, 8>;
|
|
IdVec idOff = vtkm::TypeTraits<IdVec>::ZeroInitialization();
|
|
IdVec idLen = vtkm::TypeTraits<IdVec>::ZeroInitialization();
|
|
IdVec idInd = vtkm::TypeTraits<IdVec>::ZeroInitialization();
|
|
vtkm::Id uniqueMats = vtkm::Id(0);
|
|
|
|
for (vtkm::IdComponent i = 0; i < numberOfPoints; i++)
|
|
{
|
|
idOff[i] = 0;
|
|
idLen[i] = lengths.Get(toReduce[i]);
|
|
idInd[i] = positions.Get(toReduce[i]);
|
|
}
|
|
|
|
vtkm::UInt8 hasWork = vtkm::UInt8(1);
|
|
while (hasWork != vtkm::UInt8(0))
|
|
{
|
|
hasWork = vtkm::UInt8(0);
|
|
vtkm::Id lowest = vtkm::Id(-1);
|
|
for (vtkm::IdComponent i = 0; i < numberOfPoints; i++)
|
|
{
|
|
if (idOff[i] < idLen[i])
|
|
{
|
|
vtkm::Id tmp = ids.Get(idInd[i] + idOff[i]);
|
|
if (lowest == vtkm::Id(-1) || tmp < lowest)
|
|
{
|
|
lowest = tmp;
|
|
hasWork = vtkm::UInt8(1);
|
|
}
|
|
}
|
|
}
|
|
if (hasWork != vtkm::UInt8(0))
|
|
{
|
|
vtkm::Float64 val = vtkm::Float64(0);
|
|
for (vtkm::IdComponent i = 0; i < numberOfPoints; i++)
|
|
{
|
|
if (idOff[i] < idLen[i])
|
|
{
|
|
vtkm::Id tmp = ids.Get(idInd[i] + idOff[i]);
|
|
if (lowest == tmp)
|
|
{
|
|
val += vfs.Get(idInd[i] + idOff[i]);
|
|
idOff[i]++;
|
|
}
|
|
}
|
|
}
|
|
outVFs.Set(localOffset + uniqueMats, val / vtkm::Float64(numberOfPoints));
|
|
outIDs.Set(localOffset + uniqueMats, lowest);
|
|
uniqueMats++;
|
|
}
|
|
}
|
|
}
|
|
};
|
|
|
|
VTKM_CONT void operator()(
|
|
const vtkm::cont::ArrayHandle<vtkm::Id, vtkm::cont::StorageTagBasic>& originalLen,
|
|
const vtkm::cont::ArrayHandle<vtkm::Id, vtkm::cont::StorageTagBasic>& originalPos,
|
|
const vtkm::cont::ArrayHandle<vtkm::Id, vtkm::cont::StorageTagBasic>& originalIDs,
|
|
const vtkm::cont::ArrayHandle<vtkm::Float64, vtkm::cont::StorageTagBasic>& originalVFs) const
|
|
{
|
|
vtkm::worklet::Keys<vtkm::Id> interpolationKeys(InCellInterpolationKeys);
|
|
vtkm::Id numberOfOriginalPos = originalLen.GetNumberOfValues();
|
|
vtkm::Id numberOfEdgePoints = EdgeInterpolationArray.GetNumberOfValues();
|
|
|
|
vtkm::cont::ArrayHandle<vtkm::Id> lengthArr;
|
|
vtkm::cont::ArrayHandle<vtkm::Id> posArr;
|
|
vtkm::cont::ArrayHandle<vtkm::Id> idArr;
|
|
vtkm::cont::ArrayHandle<vtkm::Float64> vfArr;
|
|
lengthArr.Allocate(numberOfOriginalPos + numberOfEdgePoints);
|
|
posArr.Allocate(numberOfOriginalPos + numberOfEdgePoints);
|
|
vtkm::cont::Algorithm::CopySubRange(originalLen, 0, numberOfOriginalPos, lengthArr);
|
|
vtkm::cont::Algorithm::CopySubRange(originalPos, 0, numberOfOriginalPos, posArr);
|
|
|
|
vtkm::cont::ArrayHandle<vtkm::Id> edgeLengths;
|
|
PerformEdgeInterpolations_C edgeCountWorklet(numberOfOriginalPos);
|
|
vtkm::worklet::DispatcherMapField<PerformEdgeInterpolations_C> edgeInterpDispatcher_C(
|
|
edgeCountWorklet);
|
|
edgeInterpDispatcher_C.Invoke(
|
|
this->EdgeInterpolationArray, lengthArr, posArr, originalIDs, edgeLengths);
|
|
|
|
vtkm::Id idLengthFromJustEdges = vtkm::cont::Algorithm::Reduce(edgeLengths, vtkm::Id(0));
|
|
|
|
idArr.Allocate(originalIDs.GetNumberOfValues() + idLengthFromJustEdges);
|
|
vfArr.Allocate(originalIDs.GetNumberOfValues() + idLengthFromJustEdges);
|
|
vtkm::cont::Algorithm::CopySubRange(originalIDs, 0, originalIDs.GetNumberOfValues(), idArr);
|
|
vtkm::cont::Algorithm::CopySubRange(originalVFs, 0, originalIDs.GetNumberOfValues(), vfArr);
|
|
vtkm::cont::Algorithm::ScanExclusive(lengthArr, posArr);
|
|
|
|
// Accept that you will have to copy data :| Maybe can speed this up with some special logic...
|
|
PerformEdgeInterpolations edgeInterpWorklet(numberOfOriginalPos);
|
|
vtkm::worklet::DispatcherMapField<PerformEdgeInterpolations> edgeInterpDispatcher(
|
|
edgeInterpWorklet);
|
|
edgeInterpDispatcher.Invoke(this->EdgeInterpolationArray, lengthArr, posArr, idArr, vfArr);
|
|
|
|
// Need to run actual edgeInterpDispatcher, we then reduce the values
|
|
|
|
|
|
vtkm::cont::ArrayHandleIndex pointArr(numberOfOriginalPos + numberOfEdgePoints);
|
|
vtkm::cont::ArrayHandle<vtkm::Id> pointArrCp;
|
|
vtkm::cont::ArrayCopy(pointArr, pointArrCp);
|
|
using IdHandle = vtkm::cont::ArrayHandle<vtkm::Id>;
|
|
using ValueHandle = vtkm::cont::ArrayHandle<vtkm::Id>;
|
|
vtkm::cont::ArrayHandlePermutation<IdHandle, ValueHandle> toReduceValues(
|
|
InCellInterpolationInfo, pointArrCp);
|
|
|
|
PerformInCellInterpolations_C incellCountWorklet;
|
|
vtkm::cont::ArrayHandle<vtkm::Id> reducedIDCounts;
|
|
vtkm::worklet::DispatcherReduceByKey<PerformInCellInterpolations_C> cellCountDispatcher(
|
|
incellCountWorklet);
|
|
cellCountDispatcher.Invoke(
|
|
interpolationKeys, toReduceValues, lengthArr, posArr, idArr, reducedIDCounts);
|
|
|
|
vtkm::cont::ArrayHandle<vtkm::Id> reducedIDOffsets;
|
|
vtkm::Id totalIDLen = vtkm::cont::Algorithm::ScanExclusive(reducedIDCounts, reducedIDOffsets);
|
|
|
|
PerformInCellInterpolations incellWorklet(originalIDs.GetNumberOfValues() +
|
|
idLengthFromJustEdges);
|
|
vtkm::cont::ArrayHandle<vtkm::Id> cellids, cellOffsets;
|
|
vtkm::cont::ArrayHandle<vtkm::Float64> cellvfs;
|
|
|
|
cellids.Allocate(totalIDLen);
|
|
cellvfs.Allocate(totalIDLen);
|
|
vtkm::worklet::DispatcherReduceByKey<PerformInCellInterpolations> cellInterpDispatcher(
|
|
incellWorklet);
|
|
cellInterpDispatcher.Invoke(interpolationKeys,
|
|
toReduceValues,
|
|
lengthArr,
|
|
posArr,
|
|
idArr,
|
|
vfArr,
|
|
reducedIDOffsets,
|
|
cellOffsets,
|
|
cellids,
|
|
cellvfs);
|
|
|
|
vtkm::Id inCellVFOffset = originalIDs.GetNumberOfValues() + idLengthFromJustEdges;
|
|
vtkm::cont::Algorithm::CopySubRange(cellids, 0, totalIDLen, idArr, inCellVFOffset);
|
|
vtkm::cont::Algorithm::CopySubRange(cellvfs, 0, totalIDLen, vfArr, inCellVFOffset);
|
|
vtkm::cont::Algorithm::CopySubRange(reducedIDCounts,
|
|
0,
|
|
reducedIDCounts.GetNumberOfValues(),
|
|
lengthArr,
|
|
numberOfOriginalPos + numberOfEdgePoints);
|
|
vtkm::cont::Algorithm::CopySubRange(cellOffsets,
|
|
0,
|
|
cellOffsets.GetNumberOfValues(),
|
|
posArr,
|
|
numberOfOriginalPos + numberOfEdgePoints);
|
|
|
|
*(this->LenOut) = lengthArr;
|
|
*(this->PosOut) = posArr;
|
|
*(this->IDOut) = idArr;
|
|
*(this->VFOut) = vfArr;
|
|
}
|
|
|
|
private:
|
|
vtkm::cont::ArrayHandle<EdgeInterpolation> EdgeInterpolationArray;
|
|
vtkm::cont::ArrayHandle<vtkm::Id> InCellInterpolationKeys;
|
|
vtkm::cont::ArrayHandle<vtkm::Id> InCellInterpolationInfo;
|
|
vtkm::Id EdgePointsOffset;
|
|
vtkm::Id InCellPointsOffset;
|
|
IDLen* LenOut;
|
|
IDPos* PosOut;
|
|
IDList* IDOut;
|
|
VFList* VFOut;
|
|
};
|
|
|
|
template <typename LookbackArr, typename WeightArr>
|
|
class InterpolateLookbackField
|
|
{
|
|
public:
|
|
InterpolateLookbackField(vtkm::cont::ArrayHandle<EdgeInterpolation> edgeInterpolationArray,
|
|
vtkm::cont::ArrayHandle<vtkm::Id> inCellInterpolationKeys,
|
|
vtkm::cont::ArrayHandle<vtkm::Id> inCellInterpolationInfo,
|
|
vtkm::Id edgePointsOffset,
|
|
vtkm::Id inCellPointsOffset,
|
|
LookbackArr* output,
|
|
WeightArr* output2)
|
|
: EdgeInterpolationArray(edgeInterpolationArray)
|
|
, InCellInterpolationKeys(inCellInterpolationKeys)
|
|
, InCellInterpolationInfo(inCellInterpolationInfo)
|
|
, EdgePointsOffset(edgePointsOffset)
|
|
, InCellPointsOffset(inCellPointsOffset)
|
|
, Output(output)
|
|
, Output2(output2)
|
|
{
|
|
}
|
|
class PerformEdgeInterpolations : public vtkm::worklet::WorkletMapField
|
|
{
|
|
public:
|
|
PerformEdgeInterpolations(vtkm::Id edgePointsOffset)
|
|
: EdgePointsOffset(edgePointsOffset)
|
|
{
|
|
}
|
|
|
|
using ControlSignature = void(FieldIn edgeInterpolations,
|
|
WholeArrayInOut inoutID,
|
|
WholeArrayInOut inoutWeights);
|
|
|
|
using ExecutionSignature = void(_1, _2, _3, WorkIndex);
|
|
|
|
template <typename EdgeInterp, typename InOutId, typename InOutWeight>
|
|
VTKM_EXEC void operator()(const EdgeInterp& ei,
|
|
InOutId& field,
|
|
InOutWeight& field1,
|
|
const vtkm::Id workIndex) const
|
|
{
|
|
|
|
vtkm::Vec<vtkm::IdComponent, 2> curOff;
|
|
vtkm::Vec<vtkm::Float64, 2> mult;
|
|
vtkm::Vec<vtkm::Id, 8> centroid;
|
|
vtkm::Vec<vtkm::Float64, 8> weight;
|
|
vtkm::Vec<vtkm::Vec<vtkm::Id, 8>, 2> keys;
|
|
vtkm::Vec<vtkm::Vec<vtkm::Float64, 8>, 2> weights;
|
|
keys[0] = field.Get(ei.Vertex1);
|
|
keys[1] = field.Get(ei.Vertex2);
|
|
weights[0] = field1.Get(ei.Vertex1);
|
|
weights[1] = field1.Get(ei.Vertex2);
|
|
for (vtkm::IdComponent i = 0; i < 8; i++)
|
|
{
|
|
weight[i] = 0;
|
|
centroid[i] = -1;
|
|
}
|
|
curOff[0] = 0;
|
|
curOff[1] = 0;
|
|
mult[0] = ei.Weight;
|
|
mult[1] = vtkm::Float64(1.0) - ei.Weight;
|
|
for (vtkm::IdComponent j = 0; j < 8; j++)
|
|
{
|
|
vtkm::Id lowestID = vtkm::Id(-1);
|
|
for (vtkm::IdComponent i = 0; i < 2; i++)
|
|
{
|
|
if (curOff[i] < 8 &&
|
|
(lowestID == vtkm::Id(-1) ||
|
|
((keys[i])[curOff[i]] != vtkm::Id(-1) && (keys[i])[curOff[i]] < lowestID)))
|
|
{
|
|
lowestID = (keys[i])[curOff[i]];
|
|
}
|
|
if (curOff[i] < 8 && (keys[i])[curOff[i]] == vtkm::Id(-1))
|
|
{
|
|
curOff[i] = 8;
|
|
}
|
|
}
|
|
if (lowestID == vtkm::Id(-1))
|
|
{
|
|
break;
|
|
}
|
|
centroid[j] = lowestID;
|
|
for (vtkm::IdComponent i = 0; i < 2; i++)
|
|
{
|
|
if (curOff[i] < 8 && lowestID == (keys[i])[curOff[i]])
|
|
{
|
|
weight[j] += mult[i] * weights[i][curOff[i]];
|
|
curOff[i]++;
|
|
}
|
|
}
|
|
}
|
|
|
|
field.Set(this->EdgePointsOffset + workIndex, centroid);
|
|
field1.Set(this->EdgePointsOffset + workIndex, weight);
|
|
}
|
|
|
|
private:
|
|
vtkm::Id EdgePointsOffset;
|
|
};
|
|
|
|
class PerformInCellInterpolations : public vtkm::worklet::WorkletReduceByKey
|
|
{
|
|
public:
|
|
using ControlSignature = void(KeysIn keys,
|
|
ValuesIn toReduceID,
|
|
WholeArrayIn Keys,
|
|
WholeArrayIn weights,
|
|
ReducedValuesOut id,
|
|
ReducedValuesOut weight);
|
|
|
|
using ExecutionSignature = void(_2, _3, _4, _5, _6);
|
|
|
|
template <typename IDs,
|
|
typename VecOfVecIDs,
|
|
typename VecOfVecWeights,
|
|
typename VecId,
|
|
typename VecWeight>
|
|
VTKM_EXEC void operator()(const IDs& ids,
|
|
const VecOfVecIDs& keysIn,
|
|
const VecOfVecWeights& weightsIn,
|
|
VecId& centroid,
|
|
VecWeight& weight) const
|
|
{
|
|
vtkm::IdComponent numValues = ids.GetNumberOfComponents();
|
|
vtkm::Vec<vtkm::IdComponent, 8> curOff;
|
|
vtkm::Vec<vtkm::Vec<vtkm::Id, 8>, 8> keys;
|
|
vtkm::Vec<vtkm::Vec<vtkm::Float64, 8>, 8> weights;
|
|
for (vtkm::IdComponent i = 0; i < 8; i++)
|
|
{
|
|
weight[i] = 0;
|
|
centroid[i] = -1;
|
|
curOff[i] = 0;
|
|
}
|
|
for (vtkm::IdComponent i = 0; i < numValues; i++)
|
|
{
|
|
keys[i] = keysIn.Get(ids[i]);
|
|
weights[i] = weightsIn.Get(ids[i]);
|
|
}
|
|
for (vtkm::IdComponent i = numValues; i < 8; i++)
|
|
{
|
|
curOff[i] = 8;
|
|
}
|
|
for (vtkm::IdComponent j = 0; j < 8; j++)
|
|
{
|
|
vtkm::Id lowestID = vtkm::Id(-1);
|
|
for (vtkm::IdComponent i = 0; i < numValues; i++)
|
|
{
|
|
auto tmp = keys[i];
|
|
if (curOff[i] < 8 &&
|
|
(lowestID == vtkm::Id(-1) ||
|
|
(tmp[curOff[i]] != vtkm::Id(-1) && tmp[curOff[i]] < lowestID)))
|
|
{
|
|
lowestID = tmp[curOff[i]];
|
|
}
|
|
|
|
if (curOff[i] < 8 && tmp[curOff[i]] == vtkm::Id(-1))
|
|
{
|
|
curOff[i] = 8;
|
|
}
|
|
}
|
|
if (lowestID == vtkm::Id(-1))
|
|
{
|
|
break;
|
|
}
|
|
centroid[j] = lowestID;
|
|
for (vtkm::IdComponent i = 0; i < numValues; i++)
|
|
{
|
|
auto tmp = keys[i];
|
|
if (curOff[i] < 8 && lowestID == tmp[curOff[i]])
|
|
{
|
|
auto w = weights[i];
|
|
weight[j] += w[curOff[i]];
|
|
curOff[i]++;
|
|
}
|
|
}
|
|
}
|
|
for (vtkm::IdComponent j = 0; j < 8; j++)
|
|
{
|
|
weight[j] *= 1. / static_cast<vtkm::Float64>(numValues);
|
|
VTKM_ASSERT(curOff[j] == 8);
|
|
}
|
|
}
|
|
};
|
|
|
|
template <typename ValueType, typename ValueType1, typename Storage, typename Storage2>
|
|
VTKM_CONT void operator()(
|
|
const vtkm::cont::ArrayHandle<ValueType, Storage>& fieldID,
|
|
const vtkm::cont::ArrayHandle<ValueType1, Storage2>& weightsField) const
|
|
{
|
|
vtkm::worklet::Keys<vtkm::Id> interpolationKeys(InCellInterpolationKeys);
|
|
|
|
vtkm::Id numberOfOriginalValues = fieldID.GetNumberOfValues();
|
|
vtkm::Id numberOfEdgePoints = EdgeInterpolationArray.GetNumberOfValues();
|
|
vtkm::Id numberOfInCellPoints = interpolationKeys.GetUniqueKeys().GetNumberOfValues();
|
|
LookbackArr result;
|
|
result.Allocate(numberOfOriginalValues + numberOfEdgePoints + numberOfInCellPoints);
|
|
vtkm::cont::Algorithm::CopySubRange(fieldID, 0, numberOfOriginalValues, result);
|
|
WeightArr result2;
|
|
result2.Allocate(numberOfOriginalValues + numberOfEdgePoints + numberOfInCellPoints);
|
|
vtkm::cont::Algorithm::CopySubRange(weightsField, 0, numberOfOriginalValues, result2);
|
|
|
|
PerformEdgeInterpolations edgeInterpWorklet(numberOfOriginalValues);
|
|
vtkm::worklet::DispatcherMapField<PerformEdgeInterpolations> edgeInterpDispatcher(
|
|
edgeInterpWorklet);
|
|
edgeInterpDispatcher.Invoke(this->EdgeInterpolationArray, result, result2);
|
|
|
|
// Perform a gather on output to get all required values for calculation of
|
|
// centroids using the interpolation info array.oi
|
|
vtkm::cont::ArrayHandleIndex nout(numberOfOriginalValues + numberOfEdgePoints);
|
|
auto toReduceValues = make_ArrayHandlePermutation(InCellInterpolationInfo, nout);
|
|
|
|
vtkm::cont::ArrayHandle<ValueType> reducedValues;
|
|
vtkm::cont::ArrayHandle<ValueType1> reducedWeights;
|
|
vtkm::worklet::DispatcherReduceByKey<PerformInCellInterpolations>
|
|
inCellInterpolationDispatcher;
|
|
inCellInterpolationDispatcher.Invoke(
|
|
interpolationKeys, toReduceValues, result, result2, reducedValues, reducedWeights);
|
|
vtkm::Id inCellPointsOffset = numberOfOriginalValues + numberOfEdgePoints;
|
|
vtkm::cont::Algorithm::CopySubRange(
|
|
reducedValues, 0, reducedValues.GetNumberOfValues(), result, inCellPointsOffset);
|
|
vtkm::cont::Algorithm::CopySubRange(
|
|
reducedWeights, 0, reducedWeights.GetNumberOfValues(), result2, inCellPointsOffset);
|
|
*(this->Output) = result;
|
|
*(this->Output2) = result2;
|
|
}
|
|
|
|
private:
|
|
vtkm::cont::ArrayHandle<EdgeInterpolation> EdgeInterpolationArray;
|
|
vtkm::cont::ArrayHandle<vtkm::Id> InCellInterpolationKeys;
|
|
vtkm::cont::ArrayHandle<vtkm::Id> InCellInterpolationInfo;
|
|
vtkm::Id EdgePointsOffset;
|
|
vtkm::Id InCellPointsOffset;
|
|
LookbackArr* Output;
|
|
WeightArr* Output2;
|
|
};
|
|
void ProcessSimpleMIRField(
|
|
const vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 8>, vtkm::cont::StorageTagBasic>& orLookback,
|
|
const vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float64, 8>, vtkm::cont::StorageTagBasic>&
|
|
orWeights,
|
|
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 8>, vtkm::cont::StorageTagBasic>& newLookback,
|
|
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float64, 8>, vtkm::cont::StorageTagBasic>& newweights)
|
|
const
|
|
{
|
|
auto worker = InterpolateLookbackField<vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id, 8>>,
|
|
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float64, 8>>>(
|
|
this->EdgePointsInterpolation,
|
|
this->InCellInterpolationKeys,
|
|
this->InCellInterpolationInfo,
|
|
this->EdgePointsOffset,
|
|
this->InCellPointsOffset,
|
|
&newLookback,
|
|
&newweights);
|
|
worker(orLookback, orWeights);
|
|
}
|
|
void ProcessMIRField(
|
|
const vtkm::cont::ArrayHandle<vtkm::Id, vtkm::cont::StorageTagBasic> orLen,
|
|
const vtkm::cont::ArrayHandle<vtkm::Id, vtkm::cont::StorageTagBasic> orPos,
|
|
const vtkm::cont::ArrayHandle<vtkm::Id, vtkm::cont::StorageTagBasic> orIDs,
|
|
const vtkm::cont::ArrayHandle<vtkm::Float64, vtkm::cont::StorageTagBasic> orVFs,
|
|
vtkm::cont::ArrayHandle<vtkm::Id, vtkm::cont::StorageTagBasic>& newLen,
|
|
vtkm::cont::ArrayHandle<vtkm::Id, vtkm::cont::StorageTagBasic>& newPos,
|
|
vtkm::cont::ArrayHandle<vtkm::Id, vtkm::cont::StorageTagBasic>& newIDs,
|
|
vtkm::cont::ArrayHandle<vtkm::Float64, vtkm::cont::StorageTagBasic>& newVFs) const
|
|
{
|
|
auto worker =
|
|
InterpolateMIRFields<vtkm::cont::ArrayHandle<vtkm::Id, vtkm::cont::StorageTagBasic>,
|
|
vtkm::cont::ArrayHandle<vtkm::Id, vtkm::cont::StorageTagBasic>,
|
|
vtkm::cont::ArrayHandle<vtkm::Id, vtkm::cont::StorageTagBasic>,
|
|
vtkm::cont::ArrayHandle<vtkm::Float64, vtkm::cont::StorageTagBasic>>(
|
|
this->EdgePointsInterpolation,
|
|
this->InCellInterpolationKeys,
|
|
this->InCellInterpolationInfo,
|
|
this->EdgePointsOffset,
|
|
this->InCellPointsOffset,
|
|
&newLen,
|
|
&newPos,
|
|
&newIDs,
|
|
&newVFs);
|
|
worker(orLen, orPos, orIDs, orVFs);
|
|
}
|
|
|
|
template <typename ValueType, typename StorageType>
|
|
vtkm::cont::ArrayHandle<ValueType> ProcessPointField(
|
|
const vtkm::cont::ArrayHandle<ValueType, StorageType>& fieldData) const
|
|
{
|
|
using ResultType = vtkm::cont::ArrayHandle<ValueType, vtkm::cont::StorageTagBasic>;
|
|
using Worker = InterpolateField<ResultType>;
|
|
ResultType output;
|
|
Worker worker = Worker(this->EdgePointsInterpolation,
|
|
this->InCellInterpolationKeys,
|
|
this->InCellInterpolationInfo,
|
|
this->EdgePointsOffset,
|
|
this->InCellPointsOffset,
|
|
&output);
|
|
worker(fieldData);
|
|
return output;
|
|
}
|
|
|
|
|
|
|
|
|
|
template <typename ValueType, typename StorageType>
|
|
vtkm::cont::ArrayHandle<ValueType> ProcessCellField(
|
|
const vtkm::cont::ArrayHandle<ValueType, StorageType>& fieldData) const
|
|
{
|
|
// Use a temporary permutation array to simplify the mapping:
|
|
auto tmp = vtkm::cont::make_ArrayHandlePermutation(this->CellMapOutputToInput, fieldData);
|
|
|
|
// Copy into an array with default storage:
|
|
vtkm::cont::ArrayHandle<ValueType> result;
|
|
vtkm::cont::ArrayCopy(tmp, result);
|
|
|
|
return result;
|
|
}
|
|
|
|
private:
|
|
MIRCases::MIRTables MIRTablesInstance;
|
|
vtkm::cont::ArrayHandle<EdgeInterpolation> EdgePointsInterpolation;
|
|
vtkm::cont::ArrayHandle<vtkm::Id> InCellInterpolationKeys;
|
|
vtkm::cont::ArrayHandle<vtkm::Id> InCellInterpolationInfo;
|
|
vtkm::cont::ArrayHandle<vtkm::Id> CellMapOutputToInput;
|
|
vtkm::Id EdgePointsOffset;
|
|
vtkm::Id InCellPointsOffset;
|
|
};
|
|
|
|
}
|
|
}
|
|
|
|
|
|
|
|
|
|
namespace vtkm
|
|
{
|
|
namespace worklet
|
|
{
|
|
template <typename IDType, typename FloatType>
|
|
struct MIRObject : public vtkm::cont::ExecutionAndControlObjectBase
|
|
{
|
|
public:
|
|
class MIRObjectPortal
|
|
{
|
|
public:
|
|
VTKM_EXEC FloatType GetVFForPoint(IDType point, IDType matID, IDType) const
|
|
{
|
|
IDType vfPos = BinarySearchForKey(
|
|
matID, this->PPos.Get(point), this->PPos.Get(point) + this->PLens.Get(point) - 1);
|
|
if (vfPos == IDType(-1))
|
|
{
|
|
for (vtkm::IdComponent s = vtkm::IdComponent(this->PPos.Get(point));
|
|
s < vtkm::IdComponent(this->PPos.Get(point) + this->PLens.Get(point));
|
|
s++)
|
|
{
|
|
if (this->PIDs.Get(s) == matID)
|
|
{
|
|
return FloatType(0);
|
|
}
|
|
}
|
|
return FloatType(0);
|
|
}
|
|
else
|
|
{
|
|
return this->PVFs.Get(vfPos);
|
|
}
|
|
}
|
|
|
|
private:
|
|
VTKM_EXEC IDType BinarySearchForKey(IDType matID, IDType low, IDType high) const
|
|
{
|
|
if (high < low)
|
|
{
|
|
return IDType(-1);
|
|
}
|
|
IDType mid = (low + high) / 2;
|
|
|
|
IDType midpoint = this->PIDs.Get(mid);
|
|
if (matID == midpoint)
|
|
{
|
|
return mid;
|
|
}
|
|
if (matID > midpoint)
|
|
{
|
|
return BinarySearchForKey(matID, (mid + 1), high);
|
|
}
|
|
else
|
|
{
|
|
return BinarySearchForKey(matID, low, mid - 1);
|
|
}
|
|
}
|
|
|
|
typename vtkm::cont::ArrayHandle<IDType, vtkm::cont::StorageTagBasic>::ReadPortalType PLens;
|
|
typename vtkm::cont::ArrayHandle<IDType, vtkm::cont::StorageTagBasic>::ReadPortalType PPos;
|
|
typename vtkm::cont::ArrayHandle<IDType, vtkm::cont::StorageTagBasic>::ReadPortalType PIDs;
|
|
typename vtkm::cont::ArrayHandle<FloatType, vtkm::cont::StorageTagBasic>::ReadPortalType PVFs;
|
|
friend struct MIRObject;
|
|
};
|
|
|
|
VTKM_CONT vtkm::cont::ArrayHandle<IDType> getPointLenArr() { return this->pointLen; }
|
|
VTKM_CONT vtkm::cont::ArrayHandle<IDType> getPointPosArr() { return this->pointPos; }
|
|
VTKM_CONT vtkm::cont::ArrayHandle<IDType> getPointIDArr() { return this->pointIDs; }
|
|
VTKM_CONT vtkm::cont::ArrayHandle<FloatType> getPointVFArr() { return this->pointVFs; }
|
|
|
|
template <typename IDInput, typename FloatInput>
|
|
MIRObject(const IDInput len, const IDInput pos, const IDInput ids, const FloatInput floats)
|
|
{
|
|
vtkm::cont::ArrayCopy(len, pointLen);
|
|
vtkm::cont::ArrayCopy(pos, pointPos);
|
|
vtkm::cont::ArrayCopy(ids, pointIDs);
|
|
vtkm::cont::ArrayCopy(floats, pointVFs);
|
|
}
|
|
|
|
MIRObjectPortal PrepareForExecution(vtkm::cont::DeviceAdapterId device, vtkm::cont::Token& token)
|
|
{
|
|
MIRObjectPortal portal;
|
|
portal.PLens = this->pointLen.PrepareForInput(device, token);
|
|
portal.PPos = this->pointPos.PrepareForInput(device, token);
|
|
portal.PIDs = this->pointIDs.PrepareForInput(device, token);
|
|
portal.PVFs = this->pointVFs.PrepareForInput(device, token);
|
|
return portal;
|
|
}
|
|
|
|
private:
|
|
vtkm::cont::ArrayHandle<IDType> pointLen, pointPos, pointIDs;
|
|
vtkm::cont::ArrayHandle<FloatType> pointVFs;
|
|
};
|
|
|
|
struct CombineVFsForPoints_C : public vtkm::worklet::WorkletVisitPointsWithCells
|
|
{
|
|
public:
|
|
using ControlSignature = void(CellSetIn cellSet,
|
|
FieldInCell lens,
|
|
FieldInCell pos,
|
|
WholeArrayIn ids,
|
|
FieldOutPoint idcount);
|
|
using ExecutionSignature = void(CellCount, _2, _3, _4, _5);
|
|
using InputDomain = _1;
|
|
|
|
template <typename LenVec, typename PosVec, typename IdsVec, typename OutVec>
|
|
VTKM_EXEC void operator()(vtkm::IdComponent numCells,
|
|
const LenVec& len,
|
|
const PosVec& pos,
|
|
const IdsVec& ids,
|
|
OutVec& outlength) const
|
|
{
|
|
|
|
// This is for the number of VFs in the surrounding cells...
|
|
// We assume that the ids are sorted.
|
|
outlength = vtkm::Id(0);
|
|
|
|
vtkm::Id uniqueMats = vtkm::Id(0);
|
|
|
|
using ida = vtkm::Id;
|
|
|
|
ida lowest = ids.Get(pos[0]);
|
|
ida prevLowest = ida(-1);
|
|
ida largest = ida(-1);
|
|
|
|
for (vtkm::IdComponent ci = 0; ci < numCells; ci++)
|
|
{
|
|
vtkm::IdComponent l = vtkm::IdComponent(pos[ci] + len[ci]);
|
|
for (vtkm::IdComponent idi = vtkm::IdComponent(pos[ci]); idi < l; idi++)
|
|
{
|
|
ida tmp = ids.Get(idi);
|
|
largest = vtkm::Maximum()(tmp, largest);
|
|
}
|
|
}
|
|
|
|
while (prevLowest != lowest)
|
|
{
|
|
for (vtkm::IdComponent ci = 0; ci < numCells; ci++)
|
|
{
|
|
vtkm::IdComponent l = vtkm::IdComponent(pos[ci] + len[ci]);
|
|
for (vtkm::IdComponent idi = vtkm::IdComponent(pos[ci]); idi < l; idi++)
|
|
{
|
|
ida tmp = ids.Get(idi);
|
|
if (tmp < lowest && tmp > prevLowest)
|
|
{
|
|
lowest = tmp;
|
|
}
|
|
}
|
|
}
|
|
uniqueMats++;
|
|
prevLowest = ida(lowest);
|
|
lowest = ida(largest);
|
|
}
|
|
outlength = uniqueMats;
|
|
}
|
|
};
|
|
struct CombineVFsForPoints : public vtkm::worklet::WorkletVisitPointsWithCells
|
|
{
|
|
public:
|
|
using ControlSignature = void(CellSetIn cellSet,
|
|
FieldInCell lens,
|
|
FieldInCell pos,
|
|
WholeArrayIn ids,
|
|
WholeArrayIn vfs,
|
|
FieldInPoint actpos,
|
|
WholeArrayOut idx,
|
|
WholeArrayOut vfx);
|
|
using ExecutionSignature = void(CellCount, _2, _3, _4, _5, _6, _7, _8);
|
|
using InputDomain = _1;
|
|
template <typename LenVec,
|
|
typename PosVec,
|
|
typename IdsVec,
|
|
typename VfsVec,
|
|
typename PosVec2,
|
|
typename OutVec,
|
|
typename OutVec2>
|
|
VTKM_EXEC void operator()(vtkm::IdComponent numCells,
|
|
const LenVec& len,
|
|
const PosVec& pos,
|
|
const IdsVec& ids,
|
|
const VfsVec& vfs,
|
|
const PosVec2& posit,
|
|
OutVec& outid,
|
|
OutVec2& outvf) const
|
|
{
|
|
|
|
// This is for the number of VFs in the surrounding cells...
|
|
// We assume that the ids are sorted.
|
|
|
|
|
|
vtkm::Id uniqueMats = vtkm::Id(0);
|
|
|
|
using ida = vtkm::Id;
|
|
|
|
ida lowest = ids.Get(pos[0]);
|
|
ida prevLowest = ida(-1);
|
|
ida largest = ida(-1);
|
|
|
|
for (vtkm::IdComponent ci = 0; ci < numCells; ci++)
|
|
{
|
|
vtkm::IdComponent l = vtkm::IdComponent(pos[ci] + len[ci]);
|
|
for (vtkm::IdComponent idi = vtkm::IdComponent(pos[ci]); idi < l; idi++)
|
|
{
|
|
ida tmp = ids.Get(idi);
|
|
largest = vtkm::Maximum()(tmp, largest);
|
|
}
|
|
}
|
|
|
|
while (prevLowest != lowest)
|
|
{
|
|
for (vtkm::IdComponent ci = 0; ci < numCells; ci++)
|
|
{
|
|
vtkm::IdComponent l = vtkm::IdComponent(pos[ci] + len[ci]);
|
|
for (vtkm::IdComponent idi = vtkm::IdComponent(pos[ci]); idi < l; idi++)
|
|
{
|
|
ida tmp = ids.Get(idi);
|
|
if (tmp < lowest && tmp > prevLowest)
|
|
{
|
|
lowest = tmp;
|
|
}
|
|
}
|
|
}
|
|
outid.Set(posit + uniqueMats, lowest);
|
|
vtkm::Float64 avg = vtkm::Float64(0);
|
|
for (vtkm::IdComponent ci = 0; ci < numCells; ci++)
|
|
{
|
|
vtkm::IdComponent l = vtkm::IdComponent(pos[ci] + len[ci]);
|
|
for (vtkm::IdComponent idi = vtkm::IdComponent(pos[ci]); idi < l; idi++)
|
|
{
|
|
ida tmp = ids.Get(idi);
|
|
if (tmp == lowest)
|
|
{
|
|
avg += vfs.Get(idi);
|
|
}
|
|
}
|
|
}
|
|
outvf.Set(posit + uniqueMats, avg / vtkm::Float64(numCells));
|
|
uniqueMats++;
|
|
prevLowest = ida(lowest);
|
|
lowest = ida(largest);
|
|
}
|
|
}
|
|
};
|
|
struct ExtractVFsForMIR_C : public vtkm::worklet::WorkletVisitCellsWithPoints
|
|
{
|
|
public:
|
|
using ControlSignature = void(CellSetIn cellSet, FieldOutCell numPointsCount);
|
|
using ExecutionSignature = void(PointCount, _2);
|
|
using InputDomain = _1;
|
|
template <typename OutVec>
|
|
VTKM_EXEC void operator()(vtkm::IdComponent numPoints, OutVec& outlength) const
|
|
{
|
|
outlength = numPoints;
|
|
}
|
|
};
|
|
struct ExtractVFsForMIR : public vtkm::worklet::WorkletVisitCellsWithPoints
|
|
{
|
|
public:
|
|
using ControlSignature = void(CellSetIn cellSet,
|
|
ExecObject mir_obj,
|
|
FieldInCell previousMatID,
|
|
FieldOutCell curvfVals,
|
|
FieldOutCell prevvfVals);
|
|
|
|
using ExecutionSignature = void(PointCount, VisitIndex, PointIndices, _2, _3, _4, _5);
|
|
using InputDomain = _1;
|
|
using ScatterType = vtkm::worklet::ScatterCounting;
|
|
template <typename CountArrayType>
|
|
VTKM_CONT static ScatterType MakeScatter(const CountArrayType& countArray)
|
|
{
|
|
VTKM_IS_ARRAY_HANDLE(CountArrayType);
|
|
return ScatterType(countArray);
|
|
}
|
|
template <typename DA, typename prevID, typename OutVec, typename OutVec2, typename pointVec>
|
|
VTKM_EXEC void operator()(vtkm::IdComponent numPoints,
|
|
vtkm::IdComponent index,
|
|
pointVec& pointIDs,
|
|
const DA& mirobj,
|
|
const prevID& previousID,
|
|
OutVec& outVF,
|
|
OutVec2& prevOutVF) const
|
|
{
|
|
(void)numPoints;
|
|
outVF = OutVec(0);
|
|
prevOutVF = OutVec2(0);
|
|
outVF = mirobj.GetVFForPoint(pointIDs[index], this->target, 0);
|
|
if (previousID == 0)
|
|
{
|
|
prevOutVF = 0;
|
|
}
|
|
else
|
|
{
|
|
prevOutVF = mirobj.GetVFForPoint(pointIDs[index], previousID, 0);
|
|
}
|
|
}
|
|
ExtractVFsForMIR(vtkm::Id targetMat)
|
|
: target(targetMat)
|
|
{
|
|
}
|
|
|
|
private:
|
|
vtkm::Id target;
|
|
};
|
|
struct CalcVol : public vtkm::worklet::WorkletVisitCellsWithPoints
|
|
{
|
|
public:
|
|
using ControlSignature = void(CellSetIn cellSet,
|
|
ExecObject mirTables,
|
|
FieldInPoint verts,
|
|
FieldOutCell vol);
|
|
using ExecutionSignature = void(PointCount, CellShape, _2, _3, _4);
|
|
template <typename Arrout, typename PointListIn, typename Dev, typename CellShape>
|
|
VTKM_EXEC void operator()(const vtkm::IdComponent pointCount,
|
|
const CellShape& cellShape,
|
|
const Dev& mirTable,
|
|
const PointListIn& vertPos,
|
|
Arrout& volumeOut) const
|
|
{
|
|
vtkm::IdComponent numFaces = mirTable.GetNumberOfFaces(static_cast<vtkm::Id>(cellShape.Id));
|
|
|
|
vtkm::Float64 totVol = vtkm::Float64(0);
|
|
vtkm::IdComponent offset = mirTable.GetFaceOffset(static_cast<vtkm::Id>(cellShape.Id));
|
|
//VTKM_LOG_S(vtkm::cont::LogLevel::Info, "CELL! " << offset << " " << numFaces);
|
|
auto av1 = vertPos[0];
|
|
for (vtkm::IdComponent i = 1; i < pointCount; i++)
|
|
{
|
|
av1 += vertPos[i];
|
|
}
|
|
auto av = av1 * (vtkm::Float64(1.0) / vtkm::Float64(pointCount));
|
|
|
|
for (vtkm::IdComponent i = 0; i < numFaces; i++)
|
|
{
|
|
vtkm::UInt8 p1 = mirTable.GetPoint(offset++);
|
|
vtkm::UInt8 p2 = mirTable.GetPoint(offset++);
|
|
vtkm::UInt8 p3 = mirTable.GetPoint(offset++);
|
|
//VTKM_LOG_S(vtkm::cont::LogLevel::Info, (i+1) << "/" << numFaces << " " << static_cast<vtkm::Id>(p1)<<" " << static_cast<vtkm::Id>(p2)<<" " << static_cast<vtkm::Id>(p3) << " " << offset);
|
|
auto v1 = vertPos[p1];
|
|
auto v2 = vertPos[p2];
|
|
auto v3 = vertPos[p3];
|
|
|
|
auto v4 = v1 - av;
|
|
auto v5 = v2 - av;
|
|
auto v6 = v3 - av;
|
|
totVol += vtkm::Abs(vtkm::Dot(v4, vtkm::Cross(v5, v6))) / 6;
|
|
}
|
|
volumeOut = totVol;
|
|
}
|
|
};
|
|
struct CalcError_C : public vtkm::worklet::WorkletReduceByKey
|
|
{
|
|
public:
|
|
using ControlSignature = void(KeysIn cellID,
|
|
ValuesIn cellColors,
|
|
WholeArrayIn origLen,
|
|
WholeArrayIn origPos,
|
|
WholeArrayIn origIDs,
|
|
WholeArrayOut newlengthsOut);
|
|
using ExecutionSignature = void(ValueCount, _1, _2, _3, _4, _5, _6);
|
|
using InputDomain = _1;
|
|
template <typename Colors, typename ORL, typename ORP, typename ORID, typename NLO>
|
|
VTKM_EXEC void operator()(const vtkm::IdComponent numCells,
|
|
const vtkm::Id cellID,
|
|
const Colors& cellCol,
|
|
const ORL& orgLen,
|
|
const ORP& orgPos,
|
|
const ORID& orgID,
|
|
NLO& outputLen) const
|
|
{
|
|
// Although I don't doubt for a minute that keys is sorted and hence the output would be too, but this ensures I don't deal with a headache if they change that.
|
|
// The orgLen and orgPos are the true, original cell IDs and VFs
|
|
// Luckily indexing into cellID should be quick compared to orgLen...
|
|
vtkm::Id lowest = orgID.Get(orgPos.Get(0));
|
|
vtkm::Id originalInd = 0;
|
|
vtkm::Id orgLen1 = orgLen.Get(cellID);
|
|
vtkm::Id orgPos1 = orgPos.Get(cellID);
|
|
vtkm::Id uniqueMats = 0;
|
|
vtkm::Id largest = orgID.Get(orgLen1 + orgPos1 - 1);
|
|
for (vtkm::IdComponent i = 0; i < numCells; i++)
|
|
{
|
|
vtkm::Id tmp = cellCol[i];
|
|
largest = vtkm::Maximum()(tmp, largest);
|
|
}
|
|
vtkm::Id prevLowest = vtkm::Id(-1);
|
|
lowest = vtkm::Id(0);
|
|
while (prevLowest != largest)
|
|
{
|
|
if (originalInd < orgLen1)
|
|
{
|
|
lowest = orgID.Get(orgPos1 + originalInd);
|
|
}
|
|
for (vtkm::IdComponent i = 0; i < numCells; i++)
|
|
{
|
|
vtkm::Id tmp = cellCol[i];
|
|
if (tmp > prevLowest)
|
|
{
|
|
lowest = vtkm::Minimum()(tmp, lowest);
|
|
}
|
|
}
|
|
if (originalInd < orgLen1)
|
|
{
|
|
if (orgID.Get(orgPos1 + originalInd) == lowest)
|
|
{
|
|
originalInd++;
|
|
}
|
|
}
|
|
uniqueMats++;
|
|
|
|
prevLowest = lowest;
|
|
lowest = largest;
|
|
}
|
|
outputLen.Set(cellID, uniqueMats);
|
|
}
|
|
};
|
|
struct CalcError : public vtkm::worklet::WorkletReduceByKey
|
|
{
|
|
private:
|
|
vtkm::Float64 lerping;
|
|
|
|
public:
|
|
CalcError(vtkm::Float64 errorLerp)
|
|
: lerping(errorLerp)
|
|
{
|
|
}
|
|
using ControlSignature = void(KeysIn cellID,
|
|
ValuesIn cellColors,
|
|
ValuesIn cellVols,
|
|
WholeArrayIn origLen,
|
|
WholeArrayIn origPos,
|
|
WholeArrayIn origIDs,
|
|
WholeArrayIn origVFs,
|
|
WholeArrayIn curLen,
|
|
WholeArrayIn curPos,
|
|
WholeArrayIn curID,
|
|
WholeArrayIn curVF,
|
|
WholeArrayIn newLength,
|
|
WholeArrayIn newposition,
|
|
WholeArrayOut newIDs,
|
|
WholeArrayOut newVFs,
|
|
WholeArrayIn origVols,
|
|
ReducedValuesOut totalErr);
|
|
using ExecutionSignature =
|
|
void(ValueCount, _1, _2, _3, _4, _5, _6, _7, _8, _9, _10, _11, _12, _13, _14, _15, _16, _17);
|
|
using InputDomain = _1;
|
|
template <typename Colors,
|
|
typename ORL,
|
|
typename ORP,
|
|
typename ORID,
|
|
typename NLO,
|
|
typename ORVF,
|
|
typename NID,
|
|
typename NVF,
|
|
typename Vols,
|
|
typename TEO,
|
|
typename CPos,
|
|
typename CLen,
|
|
typename CID,
|
|
typename CVF,
|
|
typename NLen,
|
|
typename OVols>
|
|
VTKM_EXEC void operator()(const vtkm::IdComponent numCells,
|
|
const vtkm::Id cellID,
|
|
const Colors& cellCol,
|
|
const Vols& cellVolumes,
|
|
const ORL& orgLen,
|
|
const ORP& orgPos,
|
|
const ORID& orgID,
|
|
const ORVF& orgVF,
|
|
const CLen& curLen,
|
|
const CPos& curPos,
|
|
const CID& curID,
|
|
const CVF& curVF,
|
|
const NLen&,
|
|
const NLO& inputPos,
|
|
NID& inputIDs,
|
|
NVF& inputVFs,
|
|
const OVols& orgVols,
|
|
TEO& totalErrorOut) const
|
|
{
|
|
// Although I don't doubt for a minute that keys is sorted and hence the output would be too, but this ensures I don't deal with a headache if they change that.
|
|
// The orgLen and orgPos are the true, original cell IDs and VFs
|
|
// Luckily indexing into cellID should be quick compared to orgLen...
|
|
//{
|
|
vtkm::Id lowest = orgID.Get(orgPos.Get(0));
|
|
vtkm::Id originalInd = 0;
|
|
vtkm::Id orgLen1 = orgLen.Get(cellID);
|
|
vtkm::Id orgPos1 = orgPos.Get(cellID);
|
|
vtkm::Id uniqueMats = 0;
|
|
vtkm::Id largest = orgID.Get(orgLen1 + orgPos1 - 1);
|
|
|
|
//vtkm::Id canConnect = vtkm::Id(0);
|
|
for (vtkm::IdComponent i = 0; i < numCells; i++)
|
|
{
|
|
vtkm::Id tmp = cellCol[i];
|
|
largest = vtkm::Maximum()(tmp, largest);
|
|
}
|
|
vtkm::Id prevLowest = vtkm::Id(-1);
|
|
|
|
vtkm::Id currentIndex = curPos.Get(cellID);
|
|
vtkm::Id currentLens = curLen.Get(cellID) + currentIndex;
|
|
//}
|
|
|
|
vtkm::Float64 totalError = vtkm::Float64(0);
|
|
while (prevLowest != largest)
|
|
{
|
|
if (originalInd < orgLen1)
|
|
{
|
|
lowest = orgID.Get(orgPos1 + originalInd);
|
|
}
|
|
for (vtkm::IdComponent i = 0; i < numCells; i++)
|
|
{
|
|
vtkm::Id tmp = cellCol[i];
|
|
if (tmp > prevLowest)
|
|
{
|
|
lowest = vtkm::Minimum()(tmp, lowest);
|
|
}
|
|
}
|
|
vtkm::Float64 totalVolForColor = vtkm::Float64(0);
|
|
for (vtkm::IdComponent i = 0; i < numCells; i++)
|
|
{
|
|
vtkm::Id tmp = cellCol[i];
|
|
if (tmp == lowest)
|
|
{
|
|
totalVolForColor += cellVolumes[i];
|
|
}
|
|
}
|
|
if (originalInd < orgLen1)
|
|
{
|
|
if (orgID.Get(orgPos1 + originalInd) == lowest)
|
|
{
|
|
totalVolForColor -= orgVF.Get(orgPos1 + originalInd) * orgVols.Get(cellID);
|
|
originalInd++;
|
|
}
|
|
}
|
|
|
|
vtkm::Float64 prevTarget = vtkm::Float64(0);
|
|
if (currentIndex < currentLens)
|
|
{
|
|
|
|
if (curID.Get(currentIndex) == lowest)
|
|
{
|
|
prevTarget = curVF.Get(currentIndex);
|
|
currentIndex++;
|
|
}
|
|
}
|
|
//vtkm::Float64 tmp = prevTarget;
|
|
prevTarget += this->lerping * (-totalVolForColor) / orgVols.Get(cellID);
|
|
totalError += vtkm::Abs(totalVolForColor);
|
|
//VTKM_LOG_S(vtkm::cont::LogLevel::Warn, "Lerping " << tmp << " -> " << prevTarget << " :| " << totalVolForColor);
|
|
//VTKM_LOG_S(vtkm::cont::LogLevel::Info, cellID << ": " << uniqueMats << " PT: " << tmp << " AVPTR "
|
|
// << totalVolForColor << " L: " << this->lerping << " and " << prevTarget << " / " << totalError
|
|
// << "\n" << inputPos.Get(cellID));
|
|
inputIDs.Set(inputPos.Get(cellID) + uniqueMats, lowest);
|
|
inputVFs.Set(inputPos.Get(cellID) + uniqueMats, prevTarget);
|
|
uniqueMats++;
|
|
|
|
prevLowest = lowest;
|
|
lowest = largest;
|
|
}
|
|
totalErrorOut = totalError;
|
|
}
|
|
};
|
|
struct CheckFor2D : public vtkm::worklet::WorkletVisitCellsWithPoints
|
|
{
|
|
using ControlSignature = void(CellSetIn,
|
|
FieldOutCell is2D,
|
|
FieldOutCell is3D,
|
|
FieldOutCell isOther);
|
|
using ExecutionSignature = void(CellShape, _2, _3, _4);
|
|
using InputDomain = _1;
|
|
template <typename OO, typename OP, typename OQ, typename SHAPE>
|
|
VTKM_EXEC void operator()(const SHAPE shape, OO& is2D, OP& is3D, OQ& isOther) const
|
|
{
|
|
is2D = vtkm::Id(0);
|
|
is3D = vtkm::Id(0);
|
|
|
|
if (shape.Id == vtkm::CellShapeIdToTag<vtkm::CELL_SHAPE_TRIANGLE>::Tag().Id ||
|
|
shape.Id == vtkm::CellShapeIdToTag<vtkm::CELL_SHAPE_POLYGON>::Tag().Id ||
|
|
shape.Id == vtkm::CellShapeIdToTag<vtkm::CELL_SHAPE_QUAD>::Tag().Id ||
|
|
shape.Id == vtkm::Id(6) // Tri strip?
|
|
|| shape.Id == vtkm::Id(8) /* Pixel? */)
|
|
{
|
|
is2D = vtkm::Id(1);
|
|
}
|
|
else if (shape.Id == vtkm::Id(0) /*Empty*/
|
|
|| shape.Id == vtkm::CellShapeIdToTag<vtkm::CELL_SHAPE_LINE>::Tag().Id ||
|
|
shape.Id == vtkm::CellShapeIdToTag<vtkm::CELL_SHAPE_POLY_LINE>::Tag().Id ||
|
|
shape.Id == vtkm::CellShapeIdToTag<vtkm::CELL_SHAPE_VERTEX>::Tag().Id ||
|
|
shape.Id == vtkm::Id(2) /* Poly Vertex? */)
|
|
{
|
|
isOther = vtkm::Id(1);
|
|
}
|
|
else if (shape.Id == vtkm::CellShapeIdToTag<vtkm::CELL_SHAPE_TETRA>::Tag().Id ||
|
|
shape.Id == vtkm::CellShapeIdToTag<vtkm::CELL_SHAPE_HEXAHEDRON>::Tag().Id ||
|
|
shape.Id == vtkm::CellShapeIdToTag<vtkm::CELL_SHAPE_WEDGE>::Tag().Id ||
|
|
shape.Id == vtkm::CellShapeIdToTag<vtkm::CELL_SHAPE_PYRAMID>::Tag().Id ||
|
|
shape.Id == vtkm::Id(11) /* Voxel? */)
|
|
{
|
|
is3D = vtkm::Id(1);
|
|
}
|
|
else
|
|
{
|
|
// Truly is other
|
|
isOther = vtkm::Id(1);
|
|
}
|
|
}
|
|
};
|
|
struct ConstructCellWeightList : public vtkm::worklet::WorkletMapField
|
|
{
|
|
using ControlSignature = void(FieldIn pointIDs, FieldOut VecLookback, FieldOut VecWeights);
|
|
using ExecutionSignature = void(InputIndex, _2, _3);
|
|
using InputDomain = _1;
|
|
template <typename VO1, typename VO2>
|
|
VTKM_EXEC void operator()(vtkm::Id& in, VO1& lookback, VO2& weights) const
|
|
{
|
|
for (vtkm::IdComponent i = 0; i < 8; i++)
|
|
{
|
|
lookback[i] = vtkm::Id(-1);
|
|
weights[i] = vtkm::Float64(0);
|
|
}
|
|
lookback[0] = in;
|
|
weights[0] = vtkm::Float64(1);
|
|
}
|
|
};
|
|
struct DestructPointWeightList : public vtkm::worklet::WorkletMapField
|
|
{
|
|
using ControlSignature = void(FieldIn pointIDs,
|
|
FieldIn pointWeights,
|
|
WholeArrayIn originalVals,
|
|
FieldOut newVals);
|
|
using ExecutionSignature = void(_1, _2, _3, _4);
|
|
using InputDomain = _1;
|
|
template <typename PID, typename PW, typename OV, typename NV>
|
|
VTKM_EXEC void operator()(const PID& pids, const PW& pws, const OV& ov, NV& newVals) const
|
|
{
|
|
VTKM_ASSERT(pids[0] != -1);
|
|
newVals = static_cast<NV>(ov.Get(pids[0]) * pws[0]);
|
|
for (vtkm::IdComponent i = 1; i < 8; i++)
|
|
{
|
|
if (pids[i] == vtkm::Id(-1))
|
|
{
|
|
break;
|
|
}
|
|
else
|
|
{
|
|
newVals += static_cast<NV>(ov.Get(pids[i]) * pws[i]);
|
|
}
|
|
}
|
|
}
|
|
};
|
|
}
|
|
|
|
}
|
|
|
|
#endif
|