Fix CellSetExplicit CellToPoint
The computation of CellToPoint was not correctly handling the case where some points are not part of any cell.
This commit is contained in:
Sujin Philip
parent
c077f16fb9
commit
46aa74519a
@ -22,16 +22,11 @@
|
||||
|
||||
#include <vtkm/CellShape.h>
|
||||
#include <vtkm/TopologyElementTag.h>
|
||||
#include <vtkm/cont/ArrayHandleConstant.h>
|
||||
#include <vtkm/cont/CellSet.h>
|
||||
#include <vtkm/cont/internal/ConnectivityExplicitInternals.h>
|
||||
#include <vtkm/exec/ConnectivityExplicit.h>
|
||||
|
||||
#include <vtkm/cont/ArrayHandleConstant.h>
|
||||
#include <vtkm/cont/ArrayHandleCounting.h>
|
||||
#include <vtkm/exec/ExecutionWholeArray.h>
|
||||
#include <vtkm/worklet/DispatcherMapField.h>
|
||||
#include <vtkm/worklet/WorkletMapField.h>
|
||||
|
||||
#include <map>
|
||||
#include <utility>
|
||||
|
||||
@ -348,98 +343,14 @@ public:
|
||||
// nothing to do
|
||||
}
|
||||
|
||||
// Worklet to expand the PointToCell numIndices array by repeating cell index
|
||||
class ExpandIndices : public vtkm::worklet::WorkletMapField
|
||||
{
|
||||
public:
|
||||
typedef void ControlSignature(FieldIn<> cellIndex,
|
||||
FieldIn<> offset,
|
||||
FieldIn<> numIndices,
|
||||
WholeArrayOut<> cellIndices);
|
||||
typedef void ExecutionSignature(_1, _2, _3, _4);
|
||||
typedef _1 InputDomain;
|
||||
|
||||
VTKM_CONT
|
||||
ExpandIndices() {}
|
||||
|
||||
template <typename PortalType>
|
||||
VTKM_EXEC void operator()(const vtkm::Id& cellIndex,
|
||||
const vtkm::Id& offset,
|
||||
const vtkm::Id& numIndices,
|
||||
const PortalType& cellIndices) const
|
||||
{
|
||||
VTKM_ASSERT(cellIndices.GetNumberOfValues() >= offset + numIndices);
|
||||
vtkm::Id startIndex = offset;
|
||||
for (vtkm::Id i = 0; i < numIndices; i++)
|
||||
{
|
||||
cellIndices.Set(startIndex++, cellIndex);
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
template <typename Device>
|
||||
VTKM_CONT void CreateConnectivity(Device,
|
||||
vtkm::TopologyElementTagCell,
|
||||
vtkm::TopologyElementTagPoint)
|
||||
{
|
||||
// PointToCell connectivity array (point indices) will be
|
||||
// transformed into the CellToPoint numIndices array using reduction
|
||||
//
|
||||
// PointToCell numIndices array using expansion will be
|
||||
// transformed into the CellToPoint connectivity array
|
||||
|
||||
if (this->CellToPoint.ElementsValid)
|
||||
{
|
||||
return;
|
||||
}
|
||||
|
||||
typedef vtkm::cont::DeviceAdapterAlgorithm<Device> Algorithm;
|
||||
|
||||
// Sizes of the PointToCell information
|
||||
vtkm::Id numberOfCells = this->GetNumberOfCells();
|
||||
vtkm::Id connectivityLength = this->PointToCell.Connectivity.GetNumberOfValues();
|
||||
|
||||
// PointToCell connectivity will be basis of CellToPoint numIndices
|
||||
vtkm::cont::ArrayHandle<vtkm::Id> pointIndices;
|
||||
Algorithm::Copy(this->PointToCell.Connectivity, pointIndices);
|
||||
|
||||
// PointToCell numIndices will be basis of CellToPoint connectivity
|
||||
|
||||
this->CellToPoint.Connectivity.Allocate(connectivityLength);
|
||||
vtkm::cont::ArrayHandleCounting<vtkm::Id> index(0, 1, numberOfCells);
|
||||
|
||||
this->PointToCell.BuildIndexOffsets(Device());
|
||||
vtkm::worklet::DispatcherMapField<ExpandIndices, Device> expandDispatcher;
|
||||
expandDispatcher.Invoke(index,
|
||||
this->PointToCell.IndexOffsets,
|
||||
this->PointToCell.NumIndices,
|
||||
this->CellToPoint.Connectivity);
|
||||
|
||||
// SortByKey where key is PointToCell connectivity and value is the expanded cellIndex
|
||||
Algorithm::SortByKey(pointIndices, this->CellToPoint.Connectivity);
|
||||
|
||||
if (this->GetNumberOfPoints() <= 0 && connectivityLength > 0)
|
||||
{
|
||||
this->NumberOfPoints = pointIndices.GetPortalControl().Get(connectivityLength - 1) + 1;
|
||||
}
|
||||
vtkm::Id numberOfPoints = this->GetNumberOfPoints();
|
||||
|
||||
// CellToPoint numIndices from the now sorted PointToCell connectivity
|
||||
vtkm::cont::ArrayHandleConstant<vtkm::Id> numArray(1, connectivityLength);
|
||||
vtkm::cont::ArrayHandle<vtkm::Id> uniquePoints;
|
||||
vtkm::cont::ArrayHandle<vtkm::Id> numIndices;
|
||||
uniquePoints.Allocate(numberOfPoints);
|
||||
numIndices.Allocate(numberOfPoints);
|
||||
|
||||
Algorithm::ReduceByKey(pointIndices, numArray, uniquePoints, numIndices, vtkm::Add());
|
||||
|
||||
// Set the CellToPoint information
|
||||
this->CellToPoint.Shapes = vtkm::cont::make_ArrayHandleConstant(
|
||||
static_cast<vtkm::UInt8>(CELL_SHAPE_VERTEX), numberOfPoints);
|
||||
Algorithm::Copy(numIndices, this->CellToPoint.NumIndices);
|
||||
|
||||
this->CellToPoint.ElementsValid = true;
|
||||
this->CellToPoint.IndexOffsetsValid = false;
|
||||
internal::ComputeCellToPointConnectivity(
|
||||
this->CellToPoint, this->PointToCell, this->GetNumberOfPoints(), Device());
|
||||
}
|
||||
|
||||
void PrintSummary(std::ostream& out) const VTKM_OVERRIDE
|
||||
|
||||
@ -20,10 +20,16 @@
|
||||
#ifndef vtk_m_cont_internal_ConnectivityExplicitInternals_h
|
||||
#define vtk_m_cont_internal_ConnectivityExplicitInternals_h
|
||||
|
||||
#include <vtkm/CellShape.h>
|
||||
#include <vtkm/cont/ArrayHandle.h>
|
||||
#include <vtkm/cont/ArrayHandleCast.h>
|
||||
#include <vtkm/cont/ArrayHandleConstant.h>
|
||||
#include <vtkm/cont/ArrayHandleCounting.h>
|
||||
#include <vtkm/cont/DeviceAdapterAlgorithm.h>
|
||||
#include <vtkm/cont/internal/DeviceAdapterError.h>
|
||||
#include <vtkm/exec/ExecutionWholeArray.h>
|
||||
#include <vtkm/worklet/DispatcherMapField.h>
|
||||
#include <vtkm/worklet/WorkletMapField.h>
|
||||
|
||||
namespace vtkm
|
||||
{
|
||||
@ -164,6 +170,115 @@ struct ConnectivityExplicitInternals
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
// Worklet to expand the PointToCell numIndices array by repeating cell index
|
||||
class ExpandIndices : public vtkm::worklet::WorkletMapField
|
||||
{
|
||||
public:
|
||||
typedef void ControlSignature(FieldIn<> cellIndex,
|
||||
FieldIn<> offset,
|
||||
FieldIn<> numIndices,
|
||||
WholeArrayOut<> cellIndices);
|
||||
typedef void ExecutionSignature(_1, _2, _3, _4);
|
||||
using InputDomain = _1;
|
||||
|
||||
VTKM_CONT
|
||||
ExpandIndices() {}
|
||||
|
||||
template <typename PortalType>
|
||||
VTKM_EXEC void operator()(const vtkm::Id& cellIndex,
|
||||
const vtkm::Id& offset,
|
||||
const vtkm::Id& numIndices,
|
||||
const PortalType& cellIndices) const
|
||||
{
|
||||
VTKM_ASSERT(cellIndices.GetNumberOfValues() >= offset + numIndices);
|
||||
vtkm::Id startIndex = offset;
|
||||
for (vtkm::Id i = 0; i < numIndices; i++)
|
||||
{
|
||||
cellIndices.Set(startIndex++, cellIndex);
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
class ScatterValues : public vtkm::worklet::WorkletMapField
|
||||
{
|
||||
public:
|
||||
typedef void ControlSignature(FieldIn<> index, FieldIn<> value, WholeArrayOut<> output);
|
||||
typedef void ExecutionSignature(_1, _2, _3);
|
||||
using InputDomain = _1;
|
||||
|
||||
template <typename T, typename PortalType>
|
||||
VTKM_EXEC void operator()(const vtkm::Id& index, const T& value, const PortalType& output) const
|
||||
{
|
||||
output.Set(index, value);
|
||||
}
|
||||
};
|
||||
|
||||
template <typename PointToCell, typename C2PShapeStorageTag, typename Device>
|
||||
void ComputeCellToPointConnectivity(ConnectivityExplicitInternals<C2PShapeStorageTag>& cell2Point,
|
||||
const PointToCell& point2Cell,
|
||||
vtkm::Id numberOfPoints,
|
||||
Device)
|
||||
{
|
||||
// PointToCell connectivity array (point indices) will be
|
||||
// transformed into the CellToPoint numIndices array using reduction
|
||||
//
|
||||
// PointToCell numIndices array using expansion will be
|
||||
// transformed into the CellToPoint connectivity array
|
||||
|
||||
if (cell2Point.ElementsValid)
|
||||
{
|
||||
return;
|
||||
}
|
||||
|
||||
using Algorithm = vtkm::cont::DeviceAdapterAlgorithm<Device>;
|
||||
|
||||
// Sizes of the PointToCell information
|
||||
vtkm::Id numberOfCells = point2Cell.NumIndices.GetNumberOfValues();
|
||||
vtkm::Id connectivityLength = point2Cell.Connectivity.GetNumberOfValues();
|
||||
|
||||
// PointToCell connectivity will be basis of CellToPoint numIndices
|
||||
vtkm::cont::ArrayHandle<vtkm::Id> pointIndices;
|
||||
Algorithm::Copy(point2Cell.Connectivity, pointIndices);
|
||||
|
||||
// PointToCell numIndices will be basis of CellToPoint connectivity
|
||||
|
||||
cell2Point.Connectivity.Allocate(connectivityLength);
|
||||
vtkm::cont::ArrayHandleCounting<vtkm::Id> index(0, 1, numberOfCells);
|
||||
|
||||
vtkm::worklet::DispatcherMapField<ExpandIndices, Device> expandDispatcher;
|
||||
expandDispatcher.Invoke(
|
||||
index, point2Cell.IndexOffsets, point2Cell.NumIndices, cell2Point.Connectivity);
|
||||
|
||||
// SortByKey where key is PointToCell connectivity and value is the expanded cellIndex
|
||||
Algorithm::SortByKey(pointIndices, cell2Point.Connectivity);
|
||||
|
||||
// CellToPoint numIndices from the now sorted PointToCell connectivity
|
||||
vtkm::cont::ArrayHandleConstant<vtkm::IdComponent> numArray(1, connectivityLength);
|
||||
vtkm::cont::ArrayHandle<vtkm::Id> uniquePoints;
|
||||
vtkm::cont::ArrayHandle<vtkm::IdComponent> numIndices;
|
||||
Algorithm::ReduceByKey(pointIndices, numArray, uniquePoints, numIndices, vtkm::Add());
|
||||
|
||||
// if not all the points have a cell
|
||||
if (uniquePoints.GetNumberOfValues() < numberOfPoints)
|
||||
{
|
||||
vtkm::cont::ArrayHandle<vtkm::IdComponent> fullNumIndices;
|
||||
Algorithm::Copy(vtkm::cont::ArrayHandleConstant<vtkm::IdComponent>(0, numberOfPoints),
|
||||
fullNumIndices);
|
||||
vtkm::worklet::DispatcherMapField<ScatterValues, Device>().Invoke(
|
||||
uniquePoints, numIndices, fullNumIndices);
|
||||
numIndices = fullNumIndices;
|
||||
}
|
||||
|
||||
// Set the CellToPoint information
|
||||
cell2Point.Shapes = vtkm::cont::make_ArrayHandleConstant(
|
||||
static_cast<vtkm::UInt8>(CELL_SHAPE_VERTEX), numberOfPoints);
|
||||
cell2Point.NumIndices = numIndices;
|
||||
|
||||
cell2Point.ElementsValid = true;
|
||||
cell2Point.IndexOffsetsValid = false;
|
||||
}
|
||||
}
|
||||
}
|
||||
} // namespace vtkm::cont::internal
|
||||
|
||||
@ -47,6 +47,7 @@ set(unit_tests
|
||||
UnitTestArrayHandleUniformPointCoordinates.cxx
|
||||
UnitTestArrayHandleConcatenate.cxx
|
||||
UnitTestArrayPortalToIterators.cxx
|
||||
UnitTestCellSetExplicit.cxx
|
||||
UnitTestContTesting.cxx
|
||||
UnitTestDataSetBuilderExplicit.cxx
|
||||
UnitTestDataSetBuilderRectilinear.cxx
|
||||
|
||||
161
vtkm/cont/testing/UnitTestCellSetExplicit.cxx
Normal file
161
vtkm/cont/testing/UnitTestCellSetExplicit.cxx
Normal file
@ -0,0 +1,161 @@
|
||||
//============================================================================
|
||||
// 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 2017 Sandia Corporation.
|
||||
// Copyright 2017 UT-Battelle, LLC.
|
||||
// Copyright 2017 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.
|
||||
//============================================================================
|
||||
#include <vtkm/cont/CellSetExplicit.h>
|
||||
|
||||
#include <vtkm/cont/ArrayHandle.h>
|
||||
#include <vtkm/cont/testing/Testing.h>
|
||||
#include <vtkm/worklet/DispatcherMapTopology.h>
|
||||
#include <vtkm/worklet/WorkletMapTopology.h>
|
||||
|
||||
namespace
|
||||
{
|
||||
|
||||
const vtkm::Id numberOfPoints = 11;
|
||||
|
||||
vtkm::UInt8 shapes[] = { static_cast<vtkm::UInt8>(vtkm::CELL_SHAPE_HEXAHEDRON),
|
||||
static_cast<vtkm::UInt8>(vtkm::CELL_SHAPE_PYRAMID),
|
||||
static_cast<vtkm::UInt8>(vtkm::CELL_SHAPE_TETRA),
|
||||
static_cast<vtkm::UInt8>(vtkm::CELL_SHAPE_WEDGE) };
|
||||
|
||||
vtkm::IdComponent numIndices[] = { 8, 5, 4, 6 };
|
||||
|
||||
vtkm::Id indexOffset[] = { 0, 8, 13, 17 };
|
||||
|
||||
vtkm::Id connectivity[] = { 0, 1, 5, 4, 3, 2, 6, 7, 1, 5, 6, 2, 8, 5, 8, 10, 6, 4, 7, 9, 5, 6, 10 };
|
||||
|
||||
template <typename T, std::size_t Length>
|
||||
vtkm::Id ArrayLength(const T (&)[Length])
|
||||
{
|
||||
return static_cast<vtkm::Id>(Length);
|
||||
}
|
||||
|
||||
// all points are part of atleast 1 cell
|
||||
vtkm::cont::CellSetExplicit<> MakeTestCellSet1()
|
||||
{
|
||||
vtkm::cont::CellSetExplicit<> cs;
|
||||
cs.Fill(numberOfPoints,
|
||||
vtkm::cont::make_ArrayHandle(shapes, 4),
|
||||
vtkm::cont::make_ArrayHandle(numIndices, 4),
|
||||
vtkm::cont::make_ArrayHandle(connectivity, ArrayLength(connectivity)),
|
||||
vtkm::cont::make_ArrayHandle(indexOffset, 4));
|
||||
return cs;
|
||||
}
|
||||
|
||||
// some points are not part of any cell
|
||||
vtkm::cont::CellSetExplicit<> MakeTestCellSet2()
|
||||
{
|
||||
vtkm::cont::CellSetExplicit<> cs;
|
||||
cs.Fill(
|
||||
numberOfPoints,
|
||||
vtkm::cont::make_ArrayHandle(shapes + 1, 2),
|
||||
vtkm::cont::make_ArrayHandle(numIndices + 1, 2),
|
||||
vtkm::cont::make_ArrayHandle(connectivity + indexOffset[1], indexOffset[3] - indexOffset[1]));
|
||||
return cs;
|
||||
}
|
||||
|
||||
struct WorkletPointToCell : public vtkm::worklet::WorkletMapPointToCell
|
||||
{
|
||||
typedef void ControlSignature(CellSetIn cellset, FieldOutCell<IdType> numPoints);
|
||||
typedef void ExecutionSignature(PointIndices, _2);
|
||||
using InputDomain = _1;
|
||||
|
||||
template <typename PointIndicesType>
|
||||
VTKM_EXEC void operator()(const PointIndicesType& pointIndices, vtkm::Id& numPoints) const
|
||||
{
|
||||
numPoints = pointIndices.GetNumberOfComponents();
|
||||
}
|
||||
};
|
||||
|
||||
struct WorkletCellToPoint : public vtkm::worklet::WorkletMapCellToPoint
|
||||
{
|
||||
typedef void ControlSignature(CellSetIn cellset, FieldOutPoint<IdType> numCells);
|
||||
typedef void ExecutionSignature(CellIndices, _2);
|
||||
using InputDomain = _1;
|
||||
|
||||
template <typename CellIndicesType>
|
||||
VTKM_EXEC void operator()(const CellIndicesType& cellIndices, vtkm::Id& numCells) const
|
||||
{
|
||||
numCells = cellIndices.GetNumberOfComponents();
|
||||
}
|
||||
};
|
||||
|
||||
void TestCellSetExplicit()
|
||||
{
|
||||
vtkm::cont::CellSetExplicit<> cellset;
|
||||
vtkm::cont::ArrayHandle<vtkm::Id> result;
|
||||
|
||||
std::cout << "----------------------------------------------------\n";
|
||||
std::cout << "Testing Case 1 (all points are part of atleast 1 cell): \n";
|
||||
cellset = MakeTestCellSet1();
|
||||
|
||||
std::cout << "\tTesting PointToCell\n";
|
||||
vtkm::worklet::DispatcherMapTopology<WorkletPointToCell>().Invoke(cellset, result);
|
||||
|
||||
VTKM_TEST_ASSERT(result.GetNumberOfValues() == cellset.GetNumberOfCells(),
|
||||
"result length not equal to number of cells");
|
||||
for (vtkm::Id i = 0; i < result.GetNumberOfValues(); ++i)
|
||||
{
|
||||
VTKM_TEST_ASSERT(result.GetPortalConstControl().Get(i) == numIndices[i], "incorrect result");
|
||||
}
|
||||
|
||||
std::cout << "\tTesting CellToPoint\n";
|
||||
vtkm::worklet::DispatcherMapTopology<WorkletCellToPoint>().Invoke(cellset, result);
|
||||
|
||||
VTKM_TEST_ASSERT(result.GetNumberOfValues() == cellset.GetNumberOfPoints(),
|
||||
"result length not equal to number of points");
|
||||
|
||||
vtkm::Id expected1[] = { 1, 2, 2, 1, 2, 4, 4, 2, 2, 1, 2 };
|
||||
for (vtkm::Id i = 0; i < result.GetNumberOfValues(); ++i)
|
||||
{
|
||||
VTKM_TEST_ASSERT(result.GetPortalConstControl().Get(i) == expected1[i], "incorrect result");
|
||||
}
|
||||
|
||||
std::cout << "----------------------------------------------------\n";
|
||||
std::cout << "Testing Case 2 (some points are not part of any cell): \n";
|
||||
cellset = MakeTestCellSet2();
|
||||
|
||||
std::cout << "\tTesting PointToCell\n";
|
||||
vtkm::worklet::DispatcherMapTopology<WorkletPointToCell>().Invoke(cellset, result);
|
||||
|
||||
VTKM_TEST_ASSERT(result.GetNumberOfValues() == cellset.GetNumberOfCells(),
|
||||
"result length not equal to number of cells");
|
||||
VTKM_TEST_ASSERT(result.GetPortalConstControl().Get(0) == numIndices[1] &&
|
||||
result.GetPortalConstControl().Get(1) == numIndices[2],
|
||||
"incorrect result");
|
||||
|
||||
std::cout << "\tTesting CellToPoint\n";
|
||||
vtkm::worklet::DispatcherMapTopology<WorkletCellToPoint>().Invoke(cellset, result);
|
||||
|
||||
VTKM_TEST_ASSERT(result.GetNumberOfValues() == cellset.GetNumberOfPoints(),
|
||||
"result length not equal to number of points");
|
||||
|
||||
vtkm::Id expected2[] = { 0, 1, 1, 0, 0, 2, 2, 0, 2, 0, 1 };
|
||||
for (vtkm::Id i = 0; i < result.GetNumberOfValues(); ++i)
|
||||
{
|
||||
VTKM_TEST_ASSERT(result.GetPortalConstControl().Get(i) == expected2[i], "incorrect result");
|
||||
}
|
||||
}
|
||||
|
||||
} // anonymous namespace
|
||||
|
||||
int UnitTestCellSetExplicit(int, char* [])
|
||||
{
|
||||
return vtkm::cont::testing::Testing::Run(TestCellSetExplicit);
|
||||
}
|
||||
Loading…
Reference in New Issue
Block a user