Merge topic 'cell-locator'
41679cb5 Add a CellLocator 02f48cfa Fix multiple declaration of DistributeCellData 9e0650ad Update Newton's Method to return solution status Acked-by: Kitware Robot <kwrobot@kitware.com> Merge-request: !957
This commit is contained in:
commit
4253d12062
@ -26,6 +26,14 @@
|
||||
namespace vtkm
|
||||
{
|
||||
|
||||
template <typename ScalarType, vtkm::IdComponent Size>
|
||||
struct NewtonsMethodResult
|
||||
{
|
||||
bool Valid;
|
||||
bool Converged;
|
||||
vtkm::Vec<ScalarType, Size> Solution;
|
||||
};
|
||||
|
||||
/// Uses Newton's method (a.k.a. Newton-Raphson method) to solve a nonlinear
|
||||
/// system of equations. This function assumes that the number of variables
|
||||
/// equals the number of equations. Newton's method operates on an iterative
|
||||
@ -41,7 +49,7 @@ template <typename ScalarType,
|
||||
vtkm::IdComponent Size,
|
||||
typename JacobianFunctor,
|
||||
typename FunctionFunctor>
|
||||
VTKM_EXEC_CONT vtkm::Vec<ScalarType, Size> NewtonsMethod(
|
||||
VTKM_EXEC_CONT NewtonsMethodResult<ScalarType, Size> NewtonsMethod(
|
||||
JacobianFunctor jacobianEvaluator,
|
||||
FunctionFunctor functionEvaluator,
|
||||
vtkm::Vec<ScalarType, Size> desiredFunctionOutput,
|
||||
@ -54,6 +62,7 @@ VTKM_EXEC_CONT vtkm::Vec<ScalarType, Size> NewtonsMethod(
|
||||
|
||||
VectorType x = initialGuess;
|
||||
|
||||
bool valid = false;
|
||||
bool converged = false;
|
||||
for (vtkm::IdComponent iteration = 0; !converged && (iteration < maxIterations); iteration++)
|
||||
{
|
||||
@ -69,9 +78,12 @@ VTKM_EXEC_CONT vtkm::Vec<ScalarType, Size> NewtonsMethod(
|
||||
MatrixType jacobian = jacobianEvaluator(x);
|
||||
VectorType currentFunctionOutput = functionEvaluator(x);
|
||||
|
||||
bool valid; // Ignored.
|
||||
VectorType deltaX =
|
||||
vtkm::SolveLinearSystem(jacobian, currentFunctionOutput - desiredFunctionOutput, valid);
|
||||
if (!valid)
|
||||
{
|
||||
break;
|
||||
}
|
||||
|
||||
x = x - deltaX;
|
||||
|
||||
@ -83,7 +95,7 @@ VTKM_EXEC_CONT vtkm::Vec<ScalarType, Size> NewtonsMethod(
|
||||
}
|
||||
|
||||
// Not checking whether converged.
|
||||
return x;
|
||||
return { valid, converged, x };
|
||||
}
|
||||
|
||||
} // namespace vtkm
|
||||
|
@ -44,6 +44,7 @@ set(headers
|
||||
ArrayHandleConcatenate.h
|
||||
ArrayRangeCompute.h
|
||||
ArrayRangeCompute.hxx
|
||||
CellLocatorTwoLevelUniformGrid.h
|
||||
CellSet.h
|
||||
CellSetExplicit.h
|
||||
CellSetListTag.h
|
||||
|
747
vtkm/cont/CellLocatorTwoLevelUniformGrid.h
Normal file
747
vtkm/cont/CellLocatorTwoLevelUniformGrid.h
Normal file
@ -0,0 +1,747 @@
|
||||
//============================================================================
|
||||
// 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 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
|
||||
// Copyright 2017 UT-Battelle, LLC.
|
||||
// Copyright 2017 Los Alamos National Security.
|
||||
//
|
||||
// Under the terms of Contract DE-NA0003525 with NTESS,
|
||||
// the U.S. Government retains certain rights in this software.
|
||||
//
|
||||
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
|
||||
// Laboratory (LANL), the U.S. Government retains certain rights in
|
||||
// this software.
|
||||
//============================================================================
|
||||
#ifndef vtk_m_cont_CellLocatorTwoLevelUniformGrid_h
|
||||
#define vtk_m_cont_CellLocatorTwoLevelUniformGrid_h
|
||||
|
||||
#include <vtkm/cont/ArrayHandle.h>
|
||||
#include <vtkm/cont/ArrayHandleConstant.h>
|
||||
#include <vtkm/cont/ArrayHandleTransform.h>
|
||||
#include <vtkm/cont/CoordinateSystem.h>
|
||||
#include <vtkm/cont/DeviceAdapterAlgorithm.h>
|
||||
#include <vtkm/cont/DynamicCellSet.h>
|
||||
|
||||
#include <vtkm/exec/CellInside.h>
|
||||
#include <vtkm/exec/ParametricCoordinates.h>
|
||||
|
||||
#include <vtkm/worklet/DispatcherMapField.h>
|
||||
#include <vtkm/worklet/DispatcherMapTopology.h>
|
||||
#include <vtkm/worklet/WorkletMapField.h>
|
||||
#include <vtkm/worklet/WorkletMapTopology.h>
|
||||
|
||||
#include <vtkm/Math.h>
|
||||
#include <vtkm/Types.h>
|
||||
#include <vtkm/VecFromPortalPermute.h>
|
||||
#include <vtkm/VecTraits.h>
|
||||
|
||||
namespace vtkm
|
||||
{
|
||||
namespace cont
|
||||
{
|
||||
|
||||
class CellLocatorTwoLevelUniformGrid
|
||||
{
|
||||
public:
|
||||
CellLocatorTwoLevelUniformGrid()
|
||||
: DensityL1(32.0f)
|
||||
, DensityL2(2.0f)
|
||||
{
|
||||
}
|
||||
|
||||
/// Get/Set the desired approximate number of cells per level 1 bin
|
||||
///
|
||||
void SetDensityL1(vtkm::FloatDefault val) { this->DensityL1 = val; }
|
||||
vtkm::FloatDefault GetDensityL1() const { return this->DensityL1; }
|
||||
|
||||
/// Get/Set the desired approximate number of cells per level 1 bin
|
||||
///
|
||||
void SetDensityL2(vtkm::FloatDefault val) { this->DensityL2 = val; }
|
||||
vtkm::FloatDefault GetDensityL2() const { return this->DensityL2; }
|
||||
|
||||
void SetCellSet(const vtkm::cont::DynamicCellSet& cellset) { this->CellSet = cellset; }
|
||||
|
||||
const vtkm::cont::DynamicCellSet& GetCellSet() const { return this->CellSet; }
|
||||
|
||||
void SetCoordinates(const vtkm::cont::CoordinateSystem& coords) { this->Coordinates = coords; }
|
||||
|
||||
const vtkm::cont::CoordinateSystem& GetCoordinates() const { return this->Coordinates; }
|
||||
|
||||
void PrintSummary(std::ostream& out) const
|
||||
{
|
||||
out << "DensityL1: " << this->DensityL1 << "\n";
|
||||
out << "DensityL2: " << this->DensityL2 << "\n";
|
||||
out << "Input CellSet: \n";
|
||||
this->CellSet.PrintSummary(out);
|
||||
out << "Input Coordinates: \n";
|
||||
this->Coordinates.PrintSummary(out);
|
||||
out << "LookupStructure:\n";
|
||||
out << " TopLevelGrid\n";
|
||||
out << " Dimensions: " << this->LookupStructure.TopLevel.Dimensions << "\n";
|
||||
out << " Origin: " << this->LookupStructure.TopLevel.Origin << "\n";
|
||||
out << " BinSize: " << this->LookupStructure.TopLevel.BinSize << "\n";
|
||||
out << " LeafDimensions:\n";
|
||||
vtkm::cont::printSummary_ArrayHandle(this->LookupStructure.LeafDimensions, out);
|
||||
out << " LeafStartIndex:\n";
|
||||
vtkm::cont::printSummary_ArrayHandle(this->LookupStructure.LeafStartIndex, out);
|
||||
out << " CellStartIndex:\n";
|
||||
vtkm::cont::printSummary_ArrayHandle(this->LookupStructure.CellStartIndex, out);
|
||||
out << " CellCount:\n";
|
||||
vtkm::cont::printSummary_ArrayHandle(this->LookupStructure.CellCount, out);
|
||||
out << " CellIds:\n";
|
||||
vtkm::cont::printSummary_ArrayHandle(this->LookupStructure.CellIds, out);
|
||||
}
|
||||
|
||||
private:
|
||||
using DimensionType = vtkm::Int16;
|
||||
using DimVec3 = vtkm::Vec<DimensionType, 3>;
|
||||
using FloatVec3 = vtkm::Vec<vtkm::FloatDefault, 3>;
|
||||
|
||||
struct BinsBBox
|
||||
{
|
||||
DimVec3 Min;
|
||||
DimVec3 Max;
|
||||
};
|
||||
|
||||
|
||||
|
||||
struct Bounds
|
||||
{
|
||||
FloatVec3 Min;
|
||||
FloatVec3 Max;
|
||||
};
|
||||
|
||||
struct Grid
|
||||
{
|
||||
DimVec3 Dimensions;
|
||||
FloatVec3 Origin;
|
||||
FloatVec3 BinSize;
|
||||
};
|
||||
|
||||
struct TwoLevelUniformGrid
|
||||
{
|
||||
Grid TopLevel;
|
||||
|
||||
vtkm::cont::ArrayHandle<DimVec3> LeafDimensions;
|
||||
vtkm::cont::ArrayHandle<vtkm::Id> LeafStartIndex;
|
||||
|
||||
vtkm::cont::ArrayHandle<vtkm::Id> CellStartIndex;
|
||||
vtkm::cont::ArrayHandle<vtkm::Id> CellCount;
|
||||
vtkm::cont::ArrayHandle<vtkm::Id> CellIds;
|
||||
};
|
||||
|
||||
VTKM_EXEC_CONT static DimVec3 ComputeGridDimension(vtkm::Id numberOfCells,
|
||||
const FloatVec3& size,
|
||||
vtkm::FloatDefault density)
|
||||
{
|
||||
vtkm::FloatDefault nsides = 1.0f;
|
||||
vtkm::FloatDefault volume = 1.0f;
|
||||
vtkm::FloatDefault maxside = vtkm::Max(size[0], vtkm::Max(size[1], size[2]));
|
||||
for (int i = 0; i < 3; ++i)
|
||||
{
|
||||
if (size[i] / maxside >= 1e-4f)
|
||||
{
|
||||
nsides += 1.0f;
|
||||
volume *= size[i];
|
||||
}
|
||||
}
|
||||
|
||||
auto r = vtkm::Pow((static_cast<vtkm::FloatDefault>(numberOfCells) / (volume * density)),
|
||||
1.0f / nsides);
|
||||
return vtkm::Max(DimVec3(1),
|
||||
DimVec3(static_cast<DimensionType>(size[0] * r),
|
||||
static_cast<DimensionType>(size[1] * r),
|
||||
static_cast<DimensionType>(size[2] * r)));
|
||||
}
|
||||
|
||||
template <typename PointsVecType>
|
||||
VTKM_EXEC static Bounds ComputeCellBounds(const PointsVecType& points)
|
||||
{
|
||||
using CoordsType = typename vtkm::VecTraits<PointsVecType>::ComponentType;
|
||||
auto numPoints = vtkm::VecTraits<PointsVecType>::GetNumberOfComponents(points);
|
||||
|
||||
CoordsType minp = points[0], maxp = points[0];
|
||||
for (vtkm::IdComponent i = 1; i < numPoints; ++i)
|
||||
{
|
||||
minp = vtkm::Min(minp, points[i]);
|
||||
maxp = vtkm::Max(maxp, points[i]);
|
||||
}
|
||||
|
||||
return { FloatVec3(minp), FloatVec3(maxp) };
|
||||
}
|
||||
|
||||
// TODO: This function may return false positives for non 3D cells as the
|
||||
// tests are done on the projection of the point on the cell. Extra checks
|
||||
// should be added to test if the point actually falls on the cell.
|
||||
template <typename CellShapeTag, typename CoordsType>
|
||||
VTKM_EXEC static bool PointInsideCell(FloatVec3 point,
|
||||
CellShapeTag cellShape,
|
||||
CoordsType cellPoints,
|
||||
const vtkm::exec::FunctorBase& worklet,
|
||||
FloatVec3& parametricCoordinates)
|
||||
{
|
||||
auto bounds = ComputeCellBounds(cellPoints);
|
||||
if (point[0] >= bounds.Min[0] && point[0] <= bounds.Max[0] && point[1] >= bounds.Min[1] &&
|
||||
point[1] <= bounds.Max[1] && point[2] >= bounds.Min[2] && point[2] <= bounds.Max[2])
|
||||
{
|
||||
bool success = false;
|
||||
parametricCoordinates = vtkm::exec::WorldCoordinatesToParametricCoordinates(
|
||||
cellPoints, point, cellShape, success, worklet);
|
||||
return success && vtkm::exec::CellInside(parametricCoordinates, cellShape);
|
||||
}
|
||||
return false;
|
||||
}
|
||||
|
||||
VTKM_EXEC static BinsBBox ComputeIntersectingBins(const Bounds cellBounds, const Grid& grid)
|
||||
{
|
||||
auto minb = (cellBounds.Min - grid.Origin) / grid.BinSize;
|
||||
auto maxb = (cellBounds.Max - grid.Origin) / grid.BinSize;
|
||||
return { vtkm::Max(DimVec3(0), DimVec3(minb)),
|
||||
vtkm::Min(grid.Dimensions - DimVec3(1), DimVec3(maxb)) };
|
||||
}
|
||||
|
||||
VTKM_EXEC static vtkm::Id GetNumberOfBins(const BinsBBox& binsBBox)
|
||||
{
|
||||
return (binsBBox.Max[0] - binsBBox.Min[0] + 1) * (binsBBox.Max[1] - binsBBox.Min[1] + 1) *
|
||||
(binsBBox.Max[2] - binsBBox.Min[2] + 1);
|
||||
}
|
||||
|
||||
class BBoxIterator
|
||||
{
|
||||
public:
|
||||
VTKM_EXEC_CONT BBoxIterator(const BinsBBox& bbox, const DimVec3& dim)
|
||||
: BBox(bbox)
|
||||
, Dim(dim)
|
||||
, Idx(bbox.Min)
|
||||
, Step(1,
|
||||
dim[0] - (bbox.Max[0] - bbox.Min[0] + 1),
|
||||
(dim[0] * dim[1]) - ((bbox.Max[1] - bbox.Min[1] + 1) * dim[0]))
|
||||
, FlatIdx(bbox.Min[0] + (dim[0] * (bbox.Min[1] + (dim[1] * bbox.Min[2]))))
|
||||
{
|
||||
}
|
||||
|
||||
VTKM_EXEC_CONT void Init()
|
||||
{
|
||||
this->Idx = this->BBox.Min;
|
||||
this->FlatIdx =
|
||||
this->Idx[0] + (this->Dim[0] * (this->Idx[1] + (this->Dim[1] * this->Idx[2])));
|
||||
}
|
||||
|
||||
VTKM_EXEC_CONT bool Done() const { return this->Idx[2] > this->BBox.Max[2]; }
|
||||
|
||||
VTKM_EXEC_CONT void Next()
|
||||
{
|
||||
++this->Idx[0];
|
||||
this->FlatIdx += this->Step[0];
|
||||
if (this->Idx[0] > this->BBox.Max[0])
|
||||
{
|
||||
this->Idx[0] = this->BBox.Min[0];
|
||||
++this->Idx[1];
|
||||
this->FlatIdx += this->Step[1];
|
||||
if (this->Idx[1] > this->BBox.Max[1])
|
||||
{
|
||||
this->Idx[1] = this->BBox.Min[1];
|
||||
++this->Idx[2];
|
||||
this->FlatIdx += this->Step[2];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
VTKM_EXEC_CONT const DimVec3& GetIdx() const { return this->Idx; }
|
||||
|
||||
VTKM_EXEC_CONT vtkm::Id GetFlatIdx() const { return this->FlatIdx; }
|
||||
|
||||
private:
|
||||
BinsBBox BBox;
|
||||
DimVec3 Dim;
|
||||
DimVec3 Idx;
|
||||
vtkm::Id3 Step;
|
||||
vtkm::Id FlatIdx;
|
||||
};
|
||||
|
||||
public:
|
||||
class CountBinsL1 : public vtkm::worklet::WorkletMapPointToCell
|
||||
{
|
||||
public:
|
||||
typedef void ControlSignature(CellSetIn cellset,
|
||||
FieldInPoint<Vec3> coords,
|
||||
FieldOutCell<IdType> bincount);
|
||||
typedef void ExecutionSignature(_2, _3);
|
||||
|
||||
CountBinsL1(const Grid& grid)
|
||||
: L1Grid(grid)
|
||||
{
|
||||
}
|
||||
|
||||
template <typename PointsVecType>
|
||||
VTKM_EXEC void operator()(const PointsVecType& points, vtkm::Id& numBins) const
|
||||
{
|
||||
auto cellBounds = ComputeCellBounds(points);
|
||||
auto binsBBox = ComputeIntersectingBins(cellBounds, this->L1Grid);
|
||||
numBins = GetNumberOfBins(binsBBox);
|
||||
}
|
||||
|
||||
private:
|
||||
Grid L1Grid;
|
||||
};
|
||||
|
||||
class FindBinsL1 : public vtkm::worklet::WorkletMapPointToCell
|
||||
{
|
||||
public:
|
||||
typedef void ControlSignature(CellSetIn cellset,
|
||||
FieldInPoint<Vec3> coords,
|
||||
FieldInCell<IdType> offsets,
|
||||
WholeArrayOut<IdType> binIds);
|
||||
typedef void ExecutionSignature(_2, _3, _4);
|
||||
|
||||
FindBinsL1(const Grid& grid)
|
||||
: L1Grid(grid)
|
||||
{
|
||||
}
|
||||
|
||||
template <typename PointsVecType, typename BinIdsPortalType>
|
||||
VTKM_EXEC void operator()(const PointsVecType& points,
|
||||
vtkm::Id offset,
|
||||
BinIdsPortalType& binIds) const
|
||||
{
|
||||
auto cellBounds = ComputeCellBounds(points);
|
||||
auto binsBBox = ComputeIntersectingBins(cellBounds, this->L1Grid);
|
||||
|
||||
for (BBoxIterator i(binsBBox, this->L1Grid.Dimensions); !i.Done(); i.Next())
|
||||
{
|
||||
binIds.Set(offset, i.GetFlatIdx());
|
||||
++offset;
|
||||
}
|
||||
}
|
||||
|
||||
private:
|
||||
Grid L1Grid;
|
||||
};
|
||||
|
||||
class GenerateBinsL1 : public vtkm::worklet::WorkletMapField
|
||||
{
|
||||
public:
|
||||
typedef void ControlSignature(FieldIn<IdType> binIds,
|
||||
FieldIn<IdType> cellCounts,
|
||||
WholeArrayOut<vtkm::ListTagBase<DimVec3>> dimensions);
|
||||
typedef void ExecutionSignature(_1, _2, _3);
|
||||
|
||||
using InputDomain = _1;
|
||||
|
||||
GenerateBinsL1(FloatVec3 size, vtkm::FloatDefault density)
|
||||
: Size(size)
|
||||
, Density(density)
|
||||
{
|
||||
}
|
||||
|
||||
template <typename OutputDimensionsPortal>
|
||||
VTKM_EXEC void operator()(vtkm::Id binId,
|
||||
vtkm::Id numCells,
|
||||
OutputDimensionsPortal& dimensions) const
|
||||
{
|
||||
dimensions.Set(binId, ComputeGridDimension(numCells, this->Size, this->Density));
|
||||
}
|
||||
|
||||
private:
|
||||
FloatVec3 Size;
|
||||
vtkm::FloatDefault Density;
|
||||
};
|
||||
|
||||
class CountBinsL2 : public vtkm::worklet::WorkletMapPointToCell
|
||||
{
|
||||
public:
|
||||
typedef void ControlSignature(CellSetIn cellset,
|
||||
FieldInPoint<Vec3> coords,
|
||||
WholeArrayIn<vtkm::ListTagBase<DimVec3>> binDimensions,
|
||||
FieldOutCell<IdType> bincount);
|
||||
typedef void ExecutionSignature(_2, _3, _4);
|
||||
|
||||
CountBinsL2(const Grid& grid)
|
||||
: L1Grid(grid)
|
||||
{
|
||||
}
|
||||
|
||||
template <typename PointsVecType, typename BinDimensionsPortalType>
|
||||
VTKM_EXEC void operator()(const PointsVecType& points,
|
||||
const BinDimensionsPortalType& binDimensions,
|
||||
vtkm::Id& numBins) const
|
||||
{
|
||||
auto cellBounds = ComputeCellBounds(points);
|
||||
auto binsBBox = ComputeIntersectingBins(cellBounds, this->L1Grid);
|
||||
|
||||
numBins = 0;
|
||||
for (BBoxIterator i(binsBBox, this->L1Grid.Dimensions); !i.Done(); i.Next())
|
||||
{
|
||||
Grid leaf;
|
||||
leaf.Dimensions = binDimensions.Get(i.GetFlatIdx());
|
||||
leaf.Origin =
|
||||
this->L1Grid.Origin + (static_cast<FloatVec3>(i.GetIdx()) * this->L1Grid.BinSize);
|
||||
leaf.BinSize = this->L1Grid.BinSize / static_cast<FloatVec3>(leaf.Dimensions);
|
||||
|
||||
auto binsBBoxL2 = ComputeIntersectingBins(cellBounds, leaf);
|
||||
numBins += GetNumberOfBins(binsBBoxL2);
|
||||
}
|
||||
}
|
||||
|
||||
private:
|
||||
Grid L1Grid;
|
||||
};
|
||||
|
||||
class FindBinsL2 : public vtkm::worklet::WorkletMapPointToCell
|
||||
{
|
||||
public:
|
||||
typedef void ControlSignature(CellSetIn cellset,
|
||||
FieldInPoint<Vec3> coords,
|
||||
WholeArrayIn<vtkm::ListTagBase<DimVec3>> binDimensions,
|
||||
WholeArrayIn<IdType> binStarts,
|
||||
FieldInCell<IdType> offsets,
|
||||
WholeArrayOut<IdType> binIds,
|
||||
WholeArrayOut<IdType> cellIds);
|
||||
typedef void ExecutionSignature(InputIndex, _2, _3, _4, _5, _6, _7);
|
||||
|
||||
FindBinsL2(const Grid& grid)
|
||||
: L1Grid(grid)
|
||||
{
|
||||
}
|
||||
|
||||
template <typename PointsVecType,
|
||||
typename BinDimensionsPortalType,
|
||||
typename BinStartsPortalType,
|
||||
typename BinIdsPortalType,
|
||||
typename CellIdsPortalType>
|
||||
VTKM_EXEC void operator()(vtkm::Id cellId,
|
||||
const PointsVecType& points,
|
||||
const BinDimensionsPortalType& binDimensions,
|
||||
const BinStartsPortalType& binStarts,
|
||||
vtkm::Id offset,
|
||||
BinIdsPortalType binIds,
|
||||
CellIdsPortalType cellIds) const
|
||||
{
|
||||
auto cellBounds = ComputeCellBounds(points);
|
||||
auto binsBBox = ComputeIntersectingBins(cellBounds, this->L1Grid);
|
||||
|
||||
for (BBoxIterator i(binsBBox, this->L1Grid.Dimensions); !i.Done(); i.Next())
|
||||
{
|
||||
Grid leaf;
|
||||
leaf.Dimensions = binDimensions.Get(i.GetFlatIdx());
|
||||
leaf.Origin =
|
||||
this->L1Grid.Origin + (static_cast<FloatVec3>(i.GetIdx()) * this->L1Grid.BinSize);
|
||||
leaf.BinSize = this->L1Grid.BinSize / static_cast<FloatVec3>(leaf.Dimensions);
|
||||
|
||||
auto binsBBoxL2 = ComputeIntersectingBins(cellBounds, leaf);
|
||||
vtkm::Id leafStart = binStarts.Get(i.GetFlatIdx());
|
||||
for (BBoxIterator j(binsBBoxL2, leaf.Dimensions); !j.Done(); j.Next())
|
||||
{
|
||||
binIds.Set(offset, leafStart + j.GetFlatIdx());
|
||||
cellIds.Set(offset, cellId);
|
||||
++offset;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
private:
|
||||
Grid L1Grid;
|
||||
};
|
||||
|
||||
class GenerateBinsL2 : public vtkm::worklet::WorkletMapField
|
||||
{
|
||||
public:
|
||||
typedef void ControlSignature(FieldIn<IdType> binIds,
|
||||
FieldIn<IdType> startsIn,
|
||||
FieldIn<IdType> countsIn,
|
||||
WholeArrayOut<IdType> startsOut,
|
||||
WholeArrayOut<IdType> countsOut);
|
||||
typedef void ExecutionSignature(_1, _2, _3, _4, _5);
|
||||
|
||||
using InputDomain = _1;
|
||||
|
||||
template <typename CellStartsPortalType, typename CellCountsPortalType>
|
||||
VTKM_EXEC void operator()(vtkm::Id binIndex,
|
||||
vtkm::Id start,
|
||||
vtkm::Id count,
|
||||
CellStartsPortalType& cellStarts,
|
||||
CellCountsPortalType& cellCounts) const
|
||||
{
|
||||
cellStarts.Set(binIndex, start);
|
||||
cellCounts.Set(binIndex, count);
|
||||
}
|
||||
};
|
||||
|
||||
struct DimensionsToCount
|
||||
{
|
||||
VTKM_EXEC vtkm::Id operator()(const DimVec3& dim) const { return dim[0] * dim[1] * dim[2]; }
|
||||
};
|
||||
|
||||
/// Builds the cell locator lookup structure
|
||||
///
|
||||
template <typename DeviceAdapter,
|
||||
typename CellSetList = VTKM_DEFAULT_CELL_SET_LIST_TAG,
|
||||
typename CoordsTypeList = VTKM_DEFAULT_COORDINATE_SYSTEM_TYPE_LIST_TAG,
|
||||
typename CoordsStorageList = VTKM_DEFAULT_COORDINATE_SYSTEM_STORAGE_LIST_TAG>
|
||||
void Build(DeviceAdapter,
|
||||
CellSetList cellSetTypes = CellSetList(),
|
||||
CoordsTypeList coordsValueTypes = CoordsTypeList(),
|
||||
CoordsStorageList coordsStorageType = CoordsStorageList())
|
||||
{
|
||||
using Algorithm = vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapter>;
|
||||
|
||||
auto cellset = this->CellSet.ResetCellSetList(cellSetTypes);
|
||||
auto points =
|
||||
this->Coordinates.GetData().ResetTypeAndStorageLists(coordsValueTypes, coordsStorageType);
|
||||
TwoLevelUniformGrid ls;
|
||||
|
||||
// 1: Compute the top level grid
|
||||
auto bounds = this->Coordinates.GetBounds(coordsValueTypes, coordsStorageType);
|
||||
FloatVec3 bmin(static_cast<vtkm::FloatDefault>(bounds.X.Min),
|
||||
static_cast<vtkm::FloatDefault>(bounds.Y.Min),
|
||||
static_cast<vtkm::FloatDefault>(bounds.Z.Min));
|
||||
FloatVec3 bmax(static_cast<vtkm::FloatDefault>(bounds.X.Max),
|
||||
static_cast<vtkm::FloatDefault>(bounds.Y.Max),
|
||||
static_cast<vtkm::FloatDefault>(bounds.Z.Max));
|
||||
auto size = bmax - bmin;
|
||||
auto fudge = vtkm::Max(FloatVec3(1e-6f), size * 1e-3f);
|
||||
size += 2.0f * fudge;
|
||||
|
||||
ls.TopLevel.Dimensions =
|
||||
ComputeGridDimension(cellset.GetNumberOfCells(), size, this->DensityL1);
|
||||
ls.TopLevel.Origin = bmin - fudge;
|
||||
ls.TopLevel.BinSize = size / static_cast<FloatVec3>(ls.TopLevel.Dimensions);
|
||||
|
||||
// 2: For each cell, find the number of top level bins they intersect
|
||||
vtkm::cont::ArrayHandle<vtkm::Id> binCounts;
|
||||
CountBinsL1 countL1(ls.TopLevel);
|
||||
vtkm::worklet::DispatcherMapTopology<CountBinsL1, DeviceAdapter>(countL1).Invoke(
|
||||
cellset, points, binCounts);
|
||||
|
||||
// 3: Total number of unique (cell, bin) pairs (for pre-allocating arrays)
|
||||
vtkm::Id numPairsL1 = Algorithm::ScanExclusive(binCounts, binCounts);
|
||||
|
||||
// 4: For each cell find the top level bins that intersect it
|
||||
vtkm::cont::ArrayHandle<vtkm::Id> binIds;
|
||||
binIds.Allocate(numPairsL1);
|
||||
FindBinsL1 findL1(ls.TopLevel);
|
||||
vtkm::worklet::DispatcherMapTopology<FindBinsL1, DeviceAdapter>(findL1).Invoke(
|
||||
cellset, points, binCounts, binIds);
|
||||
binCounts.ReleaseResources();
|
||||
|
||||
// 5: From above, find the number of cells that intersect each top level bin
|
||||
Algorithm::Sort(binIds);
|
||||
vtkm::cont::ArrayHandle<vtkm::Id> bins;
|
||||
vtkm::cont::ArrayHandle<vtkm::Id> cellsPerBin;
|
||||
Algorithm::ReduceByKey(binIds,
|
||||
vtkm::cont::make_ArrayHandleConstant(vtkm::Id(1), numPairsL1),
|
||||
bins,
|
||||
cellsPerBin,
|
||||
vtkm::Sum());
|
||||
binIds.ReleaseResources();
|
||||
|
||||
// 6: Compute level-2 dimensions
|
||||
vtkm::Id numberOfBins =
|
||||
ls.TopLevel.Dimensions[0] * ls.TopLevel.Dimensions[1] * ls.TopLevel.Dimensions[2];
|
||||
Algorithm::Copy(vtkm::cont::make_ArrayHandleConstant(DimVec3(0), numberOfBins),
|
||||
ls.LeafDimensions);
|
||||
GenerateBinsL1 generateL1(ls.TopLevel.BinSize, this->DensityL2);
|
||||
vtkm::worklet::DispatcherMapField<GenerateBinsL1, DeviceAdapter>(generateL1)
|
||||
.Invoke(bins, cellsPerBin, ls.LeafDimensions);
|
||||
bins.ReleaseResources();
|
||||
cellsPerBin.ReleaseResources();
|
||||
|
||||
// 7: Compute number of level-2 bins
|
||||
vtkm::Id numberOfLeaves = Algorithm::ScanExclusive(
|
||||
vtkm::cont::make_ArrayHandleTransform(ls.LeafDimensions, DimensionsToCount()),
|
||||
ls.LeafStartIndex);
|
||||
|
||||
|
||||
// 8: For each cell, find the number of l2 bins they intersect
|
||||
CountBinsL2 countL2(ls.TopLevel);
|
||||
vtkm::worklet::DispatcherMapTopology<CountBinsL2, DeviceAdapter>(countL2).Invoke(
|
||||
cellset, points, ls.LeafDimensions, binCounts);
|
||||
|
||||
// 9: Total number of unique (cell, bin) pairs (for pre-allocating arrays)
|
||||
vtkm::Id numPairsL2 = Algorithm::ScanExclusive(binCounts, binCounts);
|
||||
|
||||
// 10: For each cell, find the l2 bins they intersect
|
||||
binIds.Allocate(numPairsL2);
|
||||
ls.CellIds.Allocate(numPairsL2);
|
||||
FindBinsL2 findL2(ls.TopLevel);
|
||||
vtkm::worklet::DispatcherMapTopology<FindBinsL2, DeviceAdapter>(findL2).Invoke(
|
||||
cellset, points, ls.LeafDimensions, ls.LeafStartIndex, binCounts, binIds, ls.CellIds);
|
||||
binCounts.ReleaseResources();
|
||||
|
||||
// 11: From above, find the cells that each l2 bin intersects
|
||||
Algorithm::SortByKey(binIds, ls.CellIds);
|
||||
Algorithm::ReduceByKey(binIds,
|
||||
vtkm::cont::make_ArrayHandleConstant(vtkm::Id(1), numPairsL2),
|
||||
bins,
|
||||
cellsPerBin,
|
||||
vtkm::Sum());
|
||||
binIds.ReleaseResources();
|
||||
|
||||
// 12: Generate the leaf bin arrays
|
||||
vtkm::cont::ArrayHandle<vtkm::Id> cellsStart;
|
||||
Algorithm::ScanExclusive(cellsPerBin, cellsStart);
|
||||
|
||||
Algorithm::Copy(vtkm::cont::ArrayHandleConstant<vtkm::Id>(0, numberOfLeaves),
|
||||
ls.CellStartIndex);
|
||||
Algorithm::Copy(vtkm::cont::ArrayHandleConstant<vtkm::Id>(0, numberOfLeaves), ls.CellCount);
|
||||
vtkm::worklet::DispatcherMapField<GenerateBinsL2, DeviceAdapter>().Invoke(
|
||||
bins, cellsStart, cellsPerBin, ls.CellStartIndex, ls.CellCount);
|
||||
|
||||
std::swap(this->LookupStructure, ls);
|
||||
}
|
||||
|
||||
template <typename DeviceAdapter>
|
||||
struct TwoLevelUniformGridExecution : public vtkm::exec::ExecutionObjectBase
|
||||
{
|
||||
template <typename T>
|
||||
using ArrayPortalConst =
|
||||
typename vtkm::cont::ArrayHandle<T>::template ExecutionTypes<DeviceAdapter>::PortalConst;
|
||||
|
||||
Grid TopLevel;
|
||||
|
||||
ArrayPortalConst<DimVec3> LeafDimensions;
|
||||
ArrayPortalConst<vtkm::Id> LeafStartIndex;
|
||||
|
||||
ArrayPortalConst<vtkm::Id> CellStartIndex;
|
||||
ArrayPortalConst<vtkm::Id> CellCount;
|
||||
ArrayPortalConst<vtkm::Id> CellIds;
|
||||
};
|
||||
|
||||
class FindCellWorklet : public vtkm::worklet::WorkletMapField
|
||||
{
|
||||
public:
|
||||
typedef void ControlSignature(FieldIn<Vec3> points,
|
||||
WholeCellSetIn<> cellSet,
|
||||
WholeArrayIn<Vec3> coordinates,
|
||||
ExecObject lookupStruct,
|
||||
FieldOut<IdType> cellIds,
|
||||
FieldOut<Vec3> parametricCoordinates);
|
||||
typedef void ExecutionSignature(_1, _2, _3, _4, _5, _6);
|
||||
|
||||
using InputDomain = _1;
|
||||
|
||||
template <typename PointType,
|
||||
typename CellSetType,
|
||||
typename CoordsPortalType,
|
||||
typename LookupStructureType>
|
||||
VTKM_EXEC void operator()(const PointType& point,
|
||||
const CellSetType& cellSet,
|
||||
const CoordsPortalType& coords,
|
||||
const LookupStructureType& lookupStruct,
|
||||
vtkm::Id& cellId,
|
||||
FloatVec3& parametricCoordinates) const
|
||||
{
|
||||
cellId = -1;
|
||||
FloatVec3 p(static_cast<FloatDefault>(point[0]),
|
||||
static_cast<FloatDefault>(point[1]),
|
||||
static_cast<FloatDefault>(point[2]));
|
||||
|
||||
const Grid& topLevelGrid = lookupStruct.TopLevel;
|
||||
|
||||
DimVec3 binId3 = static_cast<DimVec3>((p - topLevelGrid.Origin) / topLevelGrid.BinSize);
|
||||
if (binId3[0] >= 0 && binId3[0] < topLevelGrid.Dimensions[0] && binId3[1] >= 0 &&
|
||||
binId3[1] < topLevelGrid.Dimensions[1] && binId3[2] >= 0 &&
|
||||
binId3[2] < topLevelGrid.Dimensions[2])
|
||||
{
|
||||
vtkm::Id binId = binId3[0] +
|
||||
(topLevelGrid.Dimensions[0] * (binId3[1] + (topLevelGrid.Dimensions[1] * binId3[2])));
|
||||
|
||||
Grid leafGrid;
|
||||
leafGrid.Dimensions = lookupStruct.LeafDimensions.Get(binId);
|
||||
if (leafGrid.Dimensions[0] + leafGrid.Dimensions[1] + leafGrid.Dimensions[2] == 0)
|
||||
{
|
||||
return;
|
||||
}
|
||||
|
||||
leafGrid.Origin =
|
||||
topLevelGrid.Origin + static_cast<FloatVec3>(binId3) * topLevelGrid.BinSize;
|
||||
leafGrid.BinSize = topLevelGrid.BinSize / static_cast<FloatVec3>(leafGrid.Dimensions);
|
||||
|
||||
vtkm::Id leafStart = lookupStruct.LeafStartIndex.Get(binId);
|
||||
DimVec3 leafId3 = static_cast<DimVec3>((p - leafGrid.Origin) / leafGrid.BinSize);
|
||||
vtkm::Id leafId = leafStart +
|
||||
(leafId3[0] +
|
||||
(leafGrid.Dimensions[0] * (leafId3[1] + (leafGrid.Dimensions[1] * leafId3[2]))));
|
||||
VTKM_ASSERT(leafId3[0] >= 0 && leafId3[0] < leafGrid.Dimensions[0] && leafId3[1] >= 0 &&
|
||||
leafId3[1] < leafGrid.Dimensions[1] && leafId3[2] >= 0 &&
|
||||
leafId3[2] < leafGrid.Dimensions[2]);
|
||||
|
||||
vtkm::Id start = lookupStruct.CellStartIndex.Get(leafId);
|
||||
vtkm::Id end = start + lookupStruct.CellCount.Get(leafId);
|
||||
for (vtkm::Id i = start; i < end; ++i)
|
||||
{
|
||||
vtkm::Id cid = lookupStruct.CellIds.Get(i);
|
||||
auto indices = cellSet.GetIndices(cid);
|
||||
vtkm::VecFromPortalPermute<decltype(indices), CoordsPortalType> pts(&indices, coords);
|
||||
FloatVec3 pc;
|
||||
if (PointInsideCell(p, cellSet.GetCellShape(cid), pts, *this, pc))
|
||||
{
|
||||
cellId = cid;
|
||||
parametricCoordinates = pc;
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
/// Finds the containing cells for the given array of points. Returns the cell ids
|
||||
/// in the `cellIds` arrays. If a cell could not be found due to the point being
|
||||
/// outside all the cells or due to numerical errors, the cell id is set to -1.
|
||||
/// Parametric coordinates of the point inside the cell is returned in the
|
||||
/// `parametricCoords` array.
|
||||
///
|
||||
template <typename PointComponentType,
|
||||
typename PointStorageType,
|
||||
typename DeviceAdapter,
|
||||
typename CellSetList = VTKM_DEFAULT_CELL_SET_LIST_TAG,
|
||||
typename CoordsTypeList = VTKM_DEFAULT_COORDINATE_SYSTEM_TYPE_LIST_TAG,
|
||||
typename CoordsStorageList = VTKM_DEFAULT_COORDINATE_SYSTEM_STORAGE_LIST_TAG>
|
||||
void FindCells(
|
||||
const vtkm::cont::ArrayHandle<vtkm::Vec<PointComponentType, 3>, PointStorageType>& points,
|
||||
vtkm::cont::ArrayHandle<vtkm::Id>& cellIds,
|
||||
vtkm::cont::ArrayHandle<FloatVec3>& parametricCoords,
|
||||
DeviceAdapter device,
|
||||
CellSetList cellSetTypes = CellSetList(),
|
||||
CoordsTypeList coordsValueTypes = CoordsTypeList(),
|
||||
CoordsStorageList coordsStorageType = CoordsStorageList()) const
|
||||
{
|
||||
vtkm::worklet::DispatcherMapField<FindCellWorklet, DeviceAdapter>().Invoke(
|
||||
points,
|
||||
this->CellSet.ResetCellSetList(cellSetTypes),
|
||||
this->Coordinates.GetData().ResetTypeAndStorageLists(coordsValueTypes, coordsStorageType),
|
||||
this->PrepareForDevice(device),
|
||||
cellIds,
|
||||
parametricCoords);
|
||||
}
|
||||
|
||||
private:
|
||||
template <typename DeviceAdapter>
|
||||
TwoLevelUniformGridExecution<DeviceAdapter> PrepareForDevice(DeviceAdapter device) const
|
||||
{
|
||||
TwoLevelUniformGridExecution<DeviceAdapter> deviceObject;
|
||||
deviceObject.TopLevel = this->LookupStructure.TopLevel;
|
||||
deviceObject.LeafDimensions = this->LookupStructure.LeafDimensions.PrepareForInput(device);
|
||||
deviceObject.LeafStartIndex = this->LookupStructure.LeafStartIndex.PrepareForInput(device);
|
||||
deviceObject.CellStartIndex = this->LookupStructure.CellStartIndex.PrepareForInput(device);
|
||||
deviceObject.CellCount = this->LookupStructure.CellCount.PrepareForInput(device);
|
||||
deviceObject.CellIds = this->LookupStructure.CellIds.PrepareForInput(device);
|
||||
return deviceObject;
|
||||
}
|
||||
|
||||
vtkm::FloatDefault DensityL1, DensityL2;
|
||||
|
||||
vtkm::cont::DynamicCellSet CellSet;
|
||||
vtkm::cont::CoordinateSystem Coordinates;
|
||||
|
||||
TwoLevelUniformGrid LookupStructure;
|
||||
};
|
||||
}
|
||||
}
|
||||
|
||||
#endif // vtk_m_cont_CellLocatorTwoLevelUniformGrid_h
|
@ -20,6 +20,7 @@
|
||||
|
||||
set(unit_tests
|
||||
UnitTestCudaArrayHandleFancy.cu
|
||||
UnitTestCudaCellLocatorTwoLevelUniformGrid.cu
|
||||
UnitTestCudaComputeRange.cu
|
||||
UnitTestCudaDataSetExplicit.cu
|
||||
UnitTestCudaDataSetSingleType.cu
|
||||
|
@ -0,0 +1,26 @@
|
||||
//============================================================================
|
||||
// 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 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
|
||||
// Copyright 2017 UT-Battelle, LLC.
|
||||
// Copyright 2017 Los Alamos National Security.
|
||||
//
|
||||
// Under the terms of Contract DE-NA0003525 with NTESS,
|
||||
// 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/testing/TestingCellLocatorTwoLevelUniformGrid.h>
|
||||
|
||||
int UnitTestCudaCellLocatorTwoLevelUniformGrid(int, char* [])
|
||||
{
|
||||
return vtkm::cont::testing::Testing::Run(
|
||||
TestingCellLocatorTwoLevelUniformGrid<vtkm::cont::DeviceAdapterTagCuda>);
|
||||
}
|
@ -21,6 +21,7 @@
|
||||
set(unit_tests
|
||||
UnitTestSerialArrayHandle.cxx
|
||||
UnitTestSerialArrayHandleFancy.cxx
|
||||
UnitTestSerialCellLocatorTwoLevelUniformGrid.cxx
|
||||
UnitTestSerialComputeRange.cxx
|
||||
UnitTestSerialDataSetExplicit.cxx
|
||||
UnitTestSerialDataSetSingleType.cxx
|
||||
|
@ -0,0 +1,26 @@
|
||||
//============================================================================
|
||||
// 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 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
|
||||
// Copyright 2017 UT-Battelle, LLC.
|
||||
// Copyright 2017 Los Alamos National Security.
|
||||
//
|
||||
// Under the terms of Contract DE-NA0003525 with NTESS,
|
||||
// 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/testing/TestingCellLocatorTwoLevelUniformGrid.h>
|
||||
|
||||
int UnitTestSerialCellLocatorTwoLevelUniformGrid(int, char* [])
|
||||
{
|
||||
return vtkm::cont::testing::Testing::Run(
|
||||
TestingCellLocatorTwoLevelUniformGrid<vtkm::cont::DeviceAdapterTagSerial>);
|
||||
}
|
@ -21,6 +21,7 @@
|
||||
set(unit_tests
|
||||
UnitTestTBBArrayHandle.cxx
|
||||
UnitTestTBBArrayHandleFancy.cxx
|
||||
UnitTestTBBCellLocatorTwoLevelUniformGrid.cxx
|
||||
UnitTestTBBComputeRange.cxx
|
||||
UnitTestTBBDataSetExplicit.cxx
|
||||
UnitTestTBBDataSetSingleType.cxx
|
||||
|
@ -0,0 +1,26 @@
|
||||
//============================================================================
|
||||
// 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 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
|
||||
// Copyright 2017 UT-Battelle, LLC.
|
||||
// Copyright 2017 Los Alamos National Security.
|
||||
//
|
||||
// Under the terms of Contract DE-NA0003525 with NTESS,
|
||||
// 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/testing/TestingCellLocatorTwoLevelUniformGrid.h>
|
||||
|
||||
int UnitTestTBBCellLocatorTwoLevelUniformGrid(int, char* [])
|
||||
{
|
||||
return vtkm::cont::testing::Testing::Run(
|
||||
TestingCellLocatorTwoLevelUniformGrid<vtkm::cont::DeviceAdapterTagTBB>);
|
||||
}
|
@ -23,6 +23,7 @@ set(headers
|
||||
MakeTestDataSet.h
|
||||
Testing.h
|
||||
TestingArrayHandles.h
|
||||
TestingCellLocatorTwoLevelUniformGrid.h
|
||||
TestingComputeRange.h
|
||||
TestingDeviceAdapter.h
|
||||
TestingDataSetExplicit.h
|
||||
|
248
vtkm/cont/testing/TestingCellLocatorTwoLevelUniformGrid.h
Normal file
248
vtkm/cont/testing/TestingCellLocatorTwoLevelUniformGrid.h
Normal file
@ -0,0 +1,248 @@
|
||||
//============================================================================
|
||||
// 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 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
|
||||
// Copyright 2017 UT-Battelle, LLC.
|
||||
// Copyright 2017 Los Alamos National Security.
|
||||
//
|
||||
// Under the terms of Contract DE-NA0003525 with NTESS,
|
||||
// the U.S. Government retains certain rights in this software.
|
||||
//
|
||||
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
|
||||
// Laboratory (LANL), the U.S. Government retains certain rights in
|
||||
// this software.
|
||||
//============================================================================
|
||||
#ifndef vtk_m_cont_testing_TestingCellLocatorTwoLevelUniformGrid_h
|
||||
#define vtk_m_cont_testing_TestingCellLocatorTwoLevelUniformGrid_h
|
||||
|
||||
#include <vtkm/cont/CellLocatorTwoLevelUniformGrid.h>
|
||||
#include <vtkm/cont/DataSetBuilderUniform.h>
|
||||
#include <vtkm/cont/testing/Testing.h>
|
||||
|
||||
#include <vtkm/exec/CellInterpolate.h>
|
||||
|
||||
#include <vtkm/worklet/DispatcherMapTopology.h>
|
||||
#include <vtkm/worklet/ScatterPermutation.h>
|
||||
#include <vtkm/worklet/Tetrahedralize.h>
|
||||
#include <vtkm/worklet/Triangulate.h>
|
||||
#include <vtkm/worklet/WorkletMapTopology.h>
|
||||
|
||||
#include <vtkm/CellShape.h>
|
||||
|
||||
#include <ctime>
|
||||
#include <random>
|
||||
|
||||
namespace
|
||||
{
|
||||
|
||||
using PointType = vtkm::Vec<vtkm::FloatDefault, 3>;
|
||||
|
||||
std::default_random_engine RandomGenerator;
|
||||
|
||||
class ParametricToWorldCoordinates : public vtkm::worklet::WorkletMapPointToCell
|
||||
{
|
||||
public:
|
||||
typedef void ControlSignature(CellSetIn cellset,
|
||||
FieldInPoint<Vec3> coords,
|
||||
FieldInOutCell<Vec3> pcs,
|
||||
FieldOutCell<Vec3> wcs);
|
||||
typedef void ExecutionSignature(CellShape, _2, _3, _4);
|
||||
|
||||
using ScatterType = vtkm::worklet::ScatterPermutation<>;
|
||||
|
||||
explicit ParametricToWorldCoordinates(const vtkm::cont::ArrayHandle<vtkm::Id>& cellIds)
|
||||
: Scatter(cellIds)
|
||||
{
|
||||
}
|
||||
|
||||
const ScatterType& GetScatter() const { return this->Scatter; }
|
||||
|
||||
template <typename CellShapeTagType, typename PointsVecType>
|
||||
VTKM_EXEC void operator()(CellShapeTagType cellShape,
|
||||
PointsVecType points,
|
||||
const PointType& pc,
|
||||
PointType& wc) const
|
||||
{
|
||||
wc = vtkm::exec::CellInterpolate(points, pc, cellShape, *this);
|
||||
}
|
||||
|
||||
private:
|
||||
ScatterType Scatter;
|
||||
};
|
||||
|
||||
template <vtkm::IdComponent DIMENSIONS, typename DeviceAdapter>
|
||||
vtkm::cont::DataSet MakeTestDataSet(const vtkm::Vec<vtkm::Id, DIMENSIONS>& dims,
|
||||
DeviceAdapter device)
|
||||
{
|
||||
using Algorithm = vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapter>;
|
||||
using Connectivity = vtkm::internal::ConnectivityStructuredInternals<DIMENSIONS>;
|
||||
|
||||
const vtkm::IdComponent PointsPerCell = 1 << DIMENSIONS;
|
||||
|
||||
auto uniformDs =
|
||||
vtkm::cont::DataSetBuilderUniform::Create(dims,
|
||||
vtkm::Vec<vtkm::FloatDefault, DIMENSIONS>(0.0f),
|
||||
vtkm::Vec<vtkm::FloatDefault, DIMENSIONS>(1.0f));
|
||||
|
||||
// copy points
|
||||
vtkm::cont::ArrayHandle<PointType> points;
|
||||
Algorithm::Copy(uniformDs.GetCoordinateSystem()
|
||||
.GetData()
|
||||
.template Cast<vtkm::cont::ArrayHandleUniformPointCoordinates>(),
|
||||
points);
|
||||
|
||||
vtkm::Id numberOfCells = uniformDs.GetCellSet().GetNumberOfCells();
|
||||
vtkm::Id numberOfIndices = numberOfCells * PointsPerCell;
|
||||
|
||||
Connectivity structured;
|
||||
structured.SetPointDimensions(dims);
|
||||
|
||||
// copy connectivity
|
||||
vtkm::cont::ArrayHandle<vtkm::Id> connectivity;
|
||||
connectivity.Allocate(numberOfIndices);
|
||||
for (vtkm::Id i = 0, idx = 0; i < numberOfCells; ++i)
|
||||
{
|
||||
auto ptids = structured.GetPointsOfCell(i);
|
||||
for (vtkm::IdComponent j = 0; j < PointsPerCell; ++j, ++idx)
|
||||
{
|
||||
connectivity.GetPortalControl().Set(idx, ptids[j]);
|
||||
}
|
||||
}
|
||||
|
||||
auto uniformCs =
|
||||
uniformDs.GetCellSet().template Cast<vtkm::cont::CellSetStructured<DIMENSIONS>>();
|
||||
vtkm::cont::CellSetSingleType<> cellset;
|
||||
|
||||
// triangulate the cellset
|
||||
switch (DIMENSIONS)
|
||||
{
|
||||
case 2:
|
||||
cellset = vtkm::worklet::Triangulate().Run(uniformCs, device);
|
||||
break;
|
||||
case 3:
|
||||
cellset = vtkm::worklet::Tetrahedralize().Run(uniformCs, device);
|
||||
break;
|
||||
default:
|
||||
VTKM_ASSERT(false);
|
||||
}
|
||||
|
||||
// It is posible that the warping will result in invalid cells. So use a
|
||||
// local random generator with a known seed that does not create invalid cells.
|
||||
std::default_random_engine rgen;
|
||||
|
||||
// Warp the coordinates
|
||||
std::uniform_real_distribution<vtkm::FloatDefault> warpFactor(-0.25f, 0.25f);
|
||||
auto pointsPortal = points.GetPortalControl();
|
||||
for (vtkm::Id i = 0; i < pointsPortal.GetNumberOfValues(); ++i)
|
||||
{
|
||||
PointType warpVec(0);
|
||||
for (vtkm::IdComponent c = 0; c < DIMENSIONS; ++c)
|
||||
{
|
||||
warpVec[c] = warpFactor(rgen);
|
||||
}
|
||||
pointsPortal.Set(i, pointsPortal.Get(i) + warpVec);
|
||||
}
|
||||
|
||||
// build dataset
|
||||
vtkm::cont::DataSet out;
|
||||
out.AddCoordinateSystem(vtkm::cont::CoordinateSystem("coords", points));
|
||||
out.AddCellSet(cellset);
|
||||
return out;
|
||||
}
|
||||
|
||||
template <vtkm::IdComponent DIMENSIONS>
|
||||
void GenerateRandomInput(const vtkm::cont::DataSet& ds,
|
||||
vtkm::Id count,
|
||||
vtkm::cont::ArrayHandle<vtkm::Id>& cellIds,
|
||||
vtkm::cont::ArrayHandle<PointType>& pcoords,
|
||||
vtkm::cont::ArrayHandle<PointType>& wcoords)
|
||||
{
|
||||
vtkm::Id numberOfCells = ds.GetCellSet().GetNumberOfCells();
|
||||
|
||||
std::uniform_int_distribution<vtkm::Id> cellIdGen(0, numberOfCells - 1);
|
||||
|
||||
cellIds.Allocate(count);
|
||||
pcoords.Allocate(count);
|
||||
wcoords.Allocate(count);
|
||||
|
||||
for (vtkm::Id i = 0; i < count; ++i)
|
||||
{
|
||||
cellIds.GetPortalControl().Set(i, cellIdGen(RandomGenerator));
|
||||
|
||||
PointType pc(0.0f);
|
||||
vtkm::FloatDefault minPc = 1e-3f;
|
||||
vtkm::FloatDefault sum = 0.0f;
|
||||
for (vtkm::IdComponent c = 0; c < DIMENSIONS; ++c)
|
||||
{
|
||||
vtkm::FloatDefault maxPc =
|
||||
1.0f - (static_cast<vtkm::FloatDefault>(DIMENSIONS - c) * minPc) - sum;
|
||||
std::uniform_real_distribution<vtkm::FloatDefault> pcoordGen(minPc, maxPc);
|
||||
pc[c] = pcoordGen(RandomGenerator);
|
||||
sum += pc[c];
|
||||
}
|
||||
pcoords.GetPortalControl().Set(i, pc);
|
||||
}
|
||||
|
||||
ParametricToWorldCoordinates pc2wc(cellIds);
|
||||
vtkm::worklet::DispatcherMapTopology<ParametricToWorldCoordinates>(pc2wc).Invoke(
|
||||
ds.GetCellSet(), ds.GetCoordinateSystem().GetData(), pcoords, wcoords);
|
||||
}
|
||||
|
||||
template <vtkm::IdComponent DIMENSIONS, typename DeviceAdapter>
|
||||
void TestCellLocator(const vtkm::Vec<vtkm::Id, DIMENSIONS>& dim,
|
||||
vtkm::Id numberOfPoints,
|
||||
DeviceAdapter device)
|
||||
{
|
||||
auto ds = MakeTestDataSet(dim, device);
|
||||
|
||||
std::cout << "Testing " << DIMENSIONS << "D dataset with " << ds.GetCellSet().GetNumberOfCells()
|
||||
<< " cells\n";
|
||||
|
||||
vtkm::cont::CellLocatorTwoLevelUniformGrid locator;
|
||||
locator.SetDensityL1(64.0f);
|
||||
locator.SetDensityL2(1.0f);
|
||||
locator.SetCellSet(ds.GetCellSet());
|
||||
locator.SetCoordinates(ds.GetCoordinateSystem());
|
||||
locator.Build(device);
|
||||
|
||||
vtkm::cont::ArrayHandle<vtkm::Id> expCellIds;
|
||||
vtkm::cont::ArrayHandle<PointType> expPCoords;
|
||||
vtkm::cont::ArrayHandle<PointType> points;
|
||||
GenerateRandomInput<DIMENSIONS>(ds, numberOfPoints, expCellIds, expPCoords, points);
|
||||
|
||||
std::cout << "Finding cells for " << numberOfPoints << " points\n";
|
||||
vtkm::cont::ArrayHandle<vtkm::Id> cellIds;
|
||||
vtkm::cont::ArrayHandle<PointType> pcoords;
|
||||
locator.FindCells(points, cellIds, pcoords, device);
|
||||
|
||||
for (vtkm::Id i = 0; i < numberOfPoints; ++i)
|
||||
{
|
||||
VTKM_TEST_ASSERT(cellIds.GetPortalConstControl().Get(i) ==
|
||||
expCellIds.GetPortalConstControl().Get(i),
|
||||
"Incorrect cell ids");
|
||||
VTKM_TEST_ASSERT(test_equal(pcoords.GetPortalConstControl().Get(i),
|
||||
expPCoords.GetPortalConstControl().Get(i),
|
||||
1e-3),
|
||||
"Incorrect parameteric coordinates");
|
||||
}
|
||||
}
|
||||
|
||||
} // anonymous
|
||||
|
||||
template <typename DeviceAdapter>
|
||||
void TestingCellLocatorTwoLevelUniformGrid()
|
||||
{
|
||||
vtkm::UInt32 seed = static_cast<vtkm::UInt32>(std::time(nullptr));
|
||||
std::cout << "Seed: " << seed << std::endl;
|
||||
RandomGenerator.seed(seed);
|
||||
|
||||
TestCellLocator(vtkm::Id3(8), 512, DeviceAdapter()); // 3D dataset
|
||||
TestCellLocator(vtkm::Id2(18), 512, DeviceAdapter()); // 2D dataset
|
||||
}
|
||||
|
||||
#endif // vtk_m_cont_testing_TestingCellLocatorTwoLevelUniformGrid_h
|
@ -23,6 +23,7 @@ set(headers
|
||||
CellDerivative.h
|
||||
CellEdge.h
|
||||
CellFace.h
|
||||
CellInside.h
|
||||
CellInterpolate.h
|
||||
ConnectivityExplicit.h
|
||||
ConnectivityPermuted.h
|
||||
|
118
vtkm/exec/CellInside.h
Normal file
118
vtkm/exec/CellInside.h
Normal file
@ -0,0 +1,118 @@
|
||||
//============================================================================
|
||||
// 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 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
|
||||
// Copyright 2017 UT-Battelle, LLC.
|
||||
// Copyright 2017 Los Alamos National Security.
|
||||
//
|
||||
// Under the terms of Contract DE-NA0003525 with NTESS,
|
||||
// the U.S. Government retains certain rights in this software.
|
||||
//
|
||||
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
|
||||
// Laboratory (LANL), the U.S. Government retains certain rights in
|
||||
// this software.
|
||||
//============================================================================
|
||||
#ifndef vtk_m_exec_CellInside_h
|
||||
#define vtk_m_exec_CellInside_h
|
||||
|
||||
#include <vtkm/CellShape.h>
|
||||
#include <vtkm/Types.h>
|
||||
|
||||
namespace vtkm
|
||||
{
|
||||
namespace exec
|
||||
{
|
||||
|
||||
template <typename T>
|
||||
static inline VTKM_EXEC bool CellInside(const vtkm::Vec<T, 3>&, vtkm::CellShapeTagEmpty)
|
||||
{
|
||||
return false;
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
static inline VTKM_EXEC bool CellInside(const vtkm::Vec<T, 3>& pcoords, vtkm::CellShapeTagVertex)
|
||||
{
|
||||
return pcoords[0] == T(0) && pcoords[1] == T(0) && pcoords[2] == T(0);
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
static inline VTKM_EXEC bool CellInside(const vtkm::Vec<T, 3>& pcoords, vtkm::CellShapeTagLine)
|
||||
{
|
||||
return pcoords[0] >= T(0) && pcoords[0] <= T(1);
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
static inline VTKM_EXEC bool CellInside(const vtkm::Vec<T, 3>& pcoords, vtkm::CellShapeTagTriangle)
|
||||
{
|
||||
return pcoords[0] >= T(0) && pcoords[1] >= T(0) && (pcoords[0] + pcoords[1] <= T(1));
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
static inline VTKM_EXEC bool CellInside(const vtkm::Vec<T, 3>& pcoords, vtkm::CellShapeTagPolygon)
|
||||
{
|
||||
return ((pcoords[0] - T(0.5)) * (pcoords[0] - T(0.5))) +
|
||||
((pcoords[1] - T(0.5)) * (pcoords[1] - T(0.5))) <=
|
||||
T(0.25);
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
static inline VTKM_EXEC bool CellInside(const vtkm::Vec<T, 3>& pcoords, vtkm::CellShapeTagQuad)
|
||||
{
|
||||
return pcoords[0] >= T(0) && pcoords[0] <= T(1) && pcoords[1] >= T(0) && pcoords[1] <= T(1);
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
static inline VTKM_EXEC bool CellInside(const vtkm::Vec<T, 3>& pcoords, vtkm::CellShapeTagTetra)
|
||||
{
|
||||
return pcoords[0] >= T(0) && pcoords[1] >= T(0) && pcoords[2] >= T(0) &&
|
||||
(pcoords[0] + pcoords[1] + pcoords[2] <= T(1));
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
static inline VTKM_EXEC bool CellInside(const vtkm::Vec<T, 3>& pcoords,
|
||||
vtkm::CellShapeTagHexahedron)
|
||||
{
|
||||
return pcoords[0] >= T(0) && pcoords[0] <= T(1) && pcoords[1] >= T(0) && pcoords[1] <= T(1) &&
|
||||
pcoords[2] >= T(0) && pcoords[2] <= T(1);
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
static inline VTKM_EXEC bool CellInside(const vtkm::Vec<T, 3>& pcoords, vtkm::CellShapeTagWedge)
|
||||
{
|
||||
return pcoords[0] >= T(0) && pcoords[1] >= T(0) && pcoords[2] >= T(0) && pcoords[2] <= T(1) &&
|
||||
(pcoords[0] + pcoords[1] <= T(1));
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
static inline VTKM_EXEC bool CellInside(const vtkm::Vec<T, 3>& pcoords, vtkm::CellShapeTagPyramid)
|
||||
{
|
||||
return pcoords[0] >= T(0) && pcoords[0] <= T(1) && pcoords[1] >= T(0) && pcoords[1] <= T(1) &&
|
||||
pcoords[2] >= T(0) && pcoords[2] <= T(1);
|
||||
}
|
||||
|
||||
/// Checks if the paramteric coordinates `pcoords` are on the inside for the
|
||||
/// specified cell type.
|
||||
///
|
||||
template <typename T>
|
||||
static inline VTKM_EXEC bool CellInside(const vtkm::Vec<T, 3>& pcoords,
|
||||
vtkm::CellShapeTagGeneric shape)
|
||||
{
|
||||
bool result = false;
|
||||
switch (shape.Id)
|
||||
{
|
||||
vtkmGenericCellShapeMacro(result = CellInside(pcoords, CellShapeTag()));
|
||||
default:
|
||||
break;
|
||||
}
|
||||
|
||||
return result;
|
||||
}
|
||||
}
|
||||
} // vtkm::exec
|
||||
|
||||
#endif // vtk_m_exec_CellInside_h
|
@ -711,13 +711,16 @@ static inline VTKM_EXEC typename WorldCoordVector::ComponentType
|
||||
WorldCoordinatesToParametricCoordinates3D(const WorldCoordVector& pointWCoords,
|
||||
const typename WorldCoordVector::ComponentType& wcoords,
|
||||
CellShapeTag,
|
||||
bool& success,
|
||||
const vtkm::exec::FunctorBase& worklet)
|
||||
{
|
||||
return vtkm::NewtonsMethod(
|
||||
auto result = vtkm::NewtonsMethod(
|
||||
JacobianFunctor3DCell<WorldCoordVector, CellShapeTag>(&pointWCoords),
|
||||
CoordinatesFunctor3DCell<WorldCoordVector, CellShapeTag>(&pointWCoords, &worklet),
|
||||
wcoords,
|
||||
typename WorldCoordVector::ComponentType(0.5f, 0.5f, 0.5f));
|
||||
success = result.Valid;
|
||||
return result.Solution;
|
||||
}
|
||||
|
||||
} // namespace detail
|
||||
@ -728,14 +731,16 @@ static inline VTKM_EXEC typename WorldCoordVector::ComponentType
|
||||
WorldCoordinatesToParametricCoordinates(const WorldCoordVector& pointWCoords,
|
||||
const typename WorldCoordVector::ComponentType& wcoords,
|
||||
vtkm::CellShapeTagGeneric shape,
|
||||
bool& success,
|
||||
const vtkm::exec::FunctorBase& worklet)
|
||||
{
|
||||
typename WorldCoordVector::ComponentType result;
|
||||
switch (shape.Id)
|
||||
{
|
||||
vtkmGenericCellShapeMacro(result = WorldCoordinatesToParametricCoordinates(
|
||||
pointWCoords, wcoords, CellShapeTag(), worklet));
|
||||
pointWCoords, wcoords, CellShapeTag(), success, worklet));
|
||||
default:
|
||||
success = false;
|
||||
worklet.RaiseError("Unknown cell shape sent to world 2 parametric.");
|
||||
return typename WorldCoordVector::ComponentType();
|
||||
}
|
||||
@ -748,9 +753,11 @@ static inline VTKM_EXEC typename WorldCoordVector::ComponentType
|
||||
WorldCoordinatesToParametricCoordinates(const WorldCoordVector&,
|
||||
const typename WorldCoordVector::ComponentType&,
|
||||
vtkm::CellShapeTagEmpty,
|
||||
bool& success,
|
||||
const vtkm::exec::FunctorBase& worklet)
|
||||
{
|
||||
worklet.RaiseError("Attempted to find point coordinates in empty cell.");
|
||||
success = false;
|
||||
return typename WorldCoordVector::ComponentType();
|
||||
}
|
||||
|
||||
@ -759,10 +766,12 @@ static inline VTKM_EXEC typename WorldCoordVector::ComponentType
|
||||
WorldCoordinatesToParametricCoordinates(const WorldCoordVector& pointWCoords,
|
||||
const typename WorldCoordVector::ComponentType&,
|
||||
vtkm::CellShapeTagVertex,
|
||||
bool& success,
|
||||
const vtkm::exec::FunctorBase& vtkmNotUsed(worklet))
|
||||
{
|
||||
(void)pointWCoords; // Silence compiler warnings.
|
||||
VTKM_ASSERT(pointWCoords.GetNumberOfComponents() == 1);
|
||||
success = true;
|
||||
return typename WorldCoordVector::ComponentType(0, 0, 0);
|
||||
}
|
||||
|
||||
@ -771,9 +780,11 @@ static inline VTKM_EXEC typename WorldCoordVector::ComponentType
|
||||
WorldCoordinatesToParametricCoordinates(const WorldCoordVector& pointWCoords,
|
||||
const typename WorldCoordVector::ComponentType& wcoords,
|
||||
vtkm::CellShapeTagLine,
|
||||
bool& success,
|
||||
const vtkm::exec::FunctorBase& vtkmNotUsed(worklet))
|
||||
{
|
||||
VTKM_ASSERT(pointWCoords.GetNumberOfComponents() == 2);
|
||||
success = true;
|
||||
|
||||
// Because this is a line, there is only one vaild parametric coordinate. Let
|
||||
// vec be the vector from the first point to the second point
|
||||
@ -798,8 +809,10 @@ static inline VTKM_EXEC vtkm::Vec<vtkm::FloatDefault, 3> WorldCoordinatesToParam
|
||||
const vtkm::VecAxisAlignedPointCoordinates<1>& pointWCoords,
|
||||
const vtkm::Vec<vtkm::FloatDefault, 3>& wcoords,
|
||||
vtkm::CellShapeTagLine,
|
||||
bool& success,
|
||||
const FunctorBase&)
|
||||
{
|
||||
success = true;
|
||||
return (wcoords - pointWCoords.GetOrigin()) / pointWCoords.GetSpacing();
|
||||
}
|
||||
|
||||
@ -808,8 +821,10 @@ static inline VTKM_EXEC typename WorldCoordVector::ComponentType
|
||||
WorldCoordinatesToParametricCoordinates(const WorldCoordVector& pointWCoords,
|
||||
const typename WorldCoordVector::ComponentType& wcoords,
|
||||
vtkm::CellShapeTagTriangle,
|
||||
bool& success,
|
||||
const vtkm::exec::FunctorBase& vtkmNotUsed(worklet))
|
||||
{
|
||||
success = true;
|
||||
return vtkm::exec::internal::ReverseInterpolateTriangle(pointWCoords, wcoords);
|
||||
}
|
||||
|
||||
@ -818,6 +833,7 @@ static inline VTKM_EXEC typename WorldCoordVector::ComponentType
|
||||
WorldCoordinatesToParametricCoordinates(const WorldCoordVector& pointWCoords,
|
||||
const typename WorldCoordVector::ComponentType& wcoords,
|
||||
vtkm::CellShapeTagPolygon,
|
||||
bool& success,
|
||||
const vtkm::exec::FunctorBase& worklet)
|
||||
{
|
||||
const vtkm::IdComponent numPoints = pointWCoords.GetNumberOfComponents();
|
||||
@ -826,16 +842,16 @@ WorldCoordinatesToParametricCoordinates(const WorldCoordVector& pointWCoords,
|
||||
{
|
||||
case 1:
|
||||
return WorldCoordinatesToParametricCoordinates(
|
||||
pointWCoords, wcoords, vtkm::CellShapeTagVertex(), worklet);
|
||||
pointWCoords, wcoords, vtkm::CellShapeTagVertex(), success, worklet);
|
||||
case 2:
|
||||
return WorldCoordinatesToParametricCoordinates(
|
||||
pointWCoords, wcoords, vtkm::CellShapeTagLine(), worklet);
|
||||
pointWCoords, wcoords, vtkm::CellShapeTagLine(), success, worklet);
|
||||
case 3:
|
||||
return WorldCoordinatesToParametricCoordinates(
|
||||
pointWCoords, wcoords, vtkm::CellShapeTagTriangle(), worklet);
|
||||
pointWCoords, wcoords, vtkm::CellShapeTagTriangle(), success, worklet);
|
||||
case 4:
|
||||
return WorldCoordinatesToParametricCoordinates(
|
||||
pointWCoords, wcoords, vtkm::CellShapeTagQuad(), worklet);
|
||||
pointWCoords, wcoords, vtkm::CellShapeTagQuad(), success, worklet);
|
||||
}
|
||||
|
||||
// If we are here, then there are 5 or more points on this polygon.
|
||||
@ -920,7 +936,7 @@ WorldCoordinatesToParametricCoordinates(const WorldCoordVector& pointWCoords,
|
||||
triangleWCoords[2] = pointWCoords[secondPointIndex];
|
||||
|
||||
WCoordType trianglePCoords = WorldCoordinatesToParametricCoordinates(
|
||||
triangleWCoords, wcoords, vtkm::CellShapeTagTriangle(), worklet);
|
||||
triangleWCoords, wcoords, vtkm::CellShapeTagTriangle(), success, worklet);
|
||||
|
||||
// trianglePCoords is in the triangle's parameter space rather than the
|
||||
// polygon's parameter space. We can find the polygon's parameter space by
|
||||
@ -940,6 +956,7 @@ static inline VTKM_EXEC typename WorldCoordVector::ComponentType
|
||||
WorldCoordinatesToParametricCoordinates(const WorldCoordVector& pointWCoords,
|
||||
const typename WorldCoordVector::ComponentType& wcoords,
|
||||
vtkm::CellShapeTagQuad,
|
||||
bool& success,
|
||||
const vtkm::exec::FunctorBase& worklet)
|
||||
{
|
||||
VTKM_ASSERT(pointWCoords.GetNumberOfComponents() == 4);
|
||||
@ -952,22 +969,25 @@ WorldCoordinatesToParametricCoordinates(const WorldCoordVector& pointWCoords,
|
||||
// plane that the polygon sits.
|
||||
vtkm::exec::internal::Space2D<T> space(pointWCoords[0], pointWCoords[1], pointWCoords[3]);
|
||||
|
||||
Vector2 pcoords = vtkm::NewtonsMethod(
|
||||
auto result = vtkm::NewtonsMethod(
|
||||
detail::JacobianFunctorQuad<WorldCoordVector, vtkm::CellShapeTagQuad>(&pointWCoords, &space),
|
||||
detail::CoordinatesFunctorQuad<WorldCoordVector, vtkm::CellShapeTagQuad>(
|
||||
&pointWCoords, &space, &worklet),
|
||||
space.ConvertCoordToSpace(wcoords),
|
||||
Vector2(0.5f, 0.5f));
|
||||
|
||||
return Vector3(pcoords[0], pcoords[1], 0);
|
||||
success = result.Valid;
|
||||
return Vector3(result.Solution[0], result.Solution[1], 0);
|
||||
}
|
||||
|
||||
static inline VTKM_EXEC vtkm::Vec<vtkm::FloatDefault, 3> WorldCoordinatesToParametricCoordinates(
|
||||
const vtkm::VecAxisAlignedPointCoordinates<2>& pointWCoords,
|
||||
const vtkm::Vec<vtkm::FloatDefault, 3>& wcoords,
|
||||
vtkm::CellShapeTagQuad,
|
||||
bool& success,
|
||||
const FunctorBase&)
|
||||
{
|
||||
success = true;
|
||||
return (wcoords - pointWCoords.GetOrigin()) / pointWCoords.GetSpacing();
|
||||
}
|
||||
|
||||
@ -976,9 +996,11 @@ static inline VTKM_EXEC typename WorldCoordVector::ComponentType
|
||||
WorldCoordinatesToParametricCoordinates(const WorldCoordVector& pointWCoords,
|
||||
const typename WorldCoordVector::ComponentType& wcoords,
|
||||
vtkm::CellShapeTagTetra,
|
||||
bool& success,
|
||||
const vtkm::exec::FunctorBase& vtkmNotUsed(worklet))
|
||||
{
|
||||
VTKM_ASSERT(pointWCoords.GetNumberOfComponents() == 4);
|
||||
success = true;
|
||||
|
||||
// We solve the world to parametric coordinates problem for tetrahedra
|
||||
// similarly to that for triangles. Before understanding this code, you
|
||||
@ -1029,20 +1051,23 @@ static inline VTKM_EXEC typename WorldCoordVector::ComponentType
|
||||
WorldCoordinatesToParametricCoordinates(const WorldCoordVector& pointWCoords,
|
||||
const typename WorldCoordVector::ComponentType& wcoords,
|
||||
vtkm::CellShapeTagHexahedron,
|
||||
bool& success,
|
||||
const vtkm::exec::FunctorBase& worklet)
|
||||
{
|
||||
VTKM_ASSERT(pointWCoords.GetNumberOfComponents() == 8);
|
||||
|
||||
return detail::WorldCoordinatesToParametricCoordinates3D(
|
||||
pointWCoords, wcoords, vtkm::CellShapeTagHexahedron(), worklet);
|
||||
pointWCoords, wcoords, vtkm::CellShapeTagHexahedron(), success, worklet);
|
||||
}
|
||||
|
||||
static inline VTKM_EXEC vtkm::Vec<vtkm::FloatDefault, 3> WorldCoordinatesToParametricCoordinates(
|
||||
const vtkm::VecAxisAlignedPointCoordinates<3>& pointWCoords,
|
||||
const vtkm::Vec<vtkm::FloatDefault, 3>& wcoords,
|
||||
vtkm::CellShapeTagHexahedron,
|
||||
bool& success,
|
||||
const FunctorBase&)
|
||||
{
|
||||
success = true;
|
||||
return (wcoords - pointWCoords.GetOrigin()) / pointWCoords.GetSpacing();
|
||||
}
|
||||
|
||||
@ -1051,12 +1076,13 @@ static inline VTKM_EXEC typename WorldCoordVector::ComponentType
|
||||
WorldCoordinatesToParametricCoordinates(const WorldCoordVector& pointWCoords,
|
||||
const typename WorldCoordVector::ComponentType& wcoords,
|
||||
vtkm::CellShapeTagWedge,
|
||||
bool& success,
|
||||
const vtkm::exec::FunctorBase& worklet)
|
||||
{
|
||||
VTKM_ASSERT(pointWCoords.GetNumberOfComponents() == 6);
|
||||
|
||||
return detail::WorldCoordinatesToParametricCoordinates3D(
|
||||
pointWCoords, wcoords, vtkm::CellShapeTagWedge(), worklet);
|
||||
pointWCoords, wcoords, vtkm::CellShapeTagWedge(), success, worklet);
|
||||
}
|
||||
|
||||
template <typename WorldCoordVector>
|
||||
@ -1064,12 +1090,13 @@ static inline VTKM_EXEC typename WorldCoordVector::ComponentType
|
||||
WorldCoordinatesToParametricCoordinates(const WorldCoordVector& pointWCoords,
|
||||
const typename WorldCoordVector::ComponentType& wcoords,
|
||||
vtkm::CellShapeTagPyramid,
|
||||
bool& success,
|
||||
const vtkm::exec::FunctorBase& worklet)
|
||||
{
|
||||
VTKM_ASSERT(pointWCoords.GetNumberOfComponents() == 5);
|
||||
|
||||
return detail::WorldCoordinatesToParametricCoordinates3D(
|
||||
pointWCoords, wcoords, vtkm::CellShapeTagPyramid(), worklet);
|
||||
pointWCoords, wcoords, vtkm::CellShapeTagPyramid(), success, worklet);
|
||||
}
|
||||
}
|
||||
} // namespace vtkm::exec
|
||||
|
@ -80,8 +80,9 @@ static void CompareCoordinates(const PointWCoordsType& pointWCoords,
|
||||
VTKM_TEST_ASSERT(test_equal(computedWCoords, trueWCoords, 0.01),
|
||||
"Computed wrong world coords from parametric coords.");
|
||||
|
||||
bool success = false;
|
||||
Vector3 computedPCoords = vtkm::exec::WorldCoordinatesToParametricCoordinates(
|
||||
pointWCoords, trueWCoords, shape, workletProxy);
|
||||
pointWCoords, trueWCoords, shape, success, workletProxy);
|
||||
VTKM_TEST_ASSERT(!errorMessage.IsErrorRaised(), messageBuffer);
|
||||
VTKM_TEST_ASSERT(test_equal(computedPCoords, truePCoords, 0.01),
|
||||
"Computed wrong parametric coords from world coords.");
|
||||
@ -166,8 +167,9 @@ void TestPCoordsSample(const PointWCoordsType& pointWCoords, CellShapeTag shape)
|
||||
Vector3 wcoords = vtkm::exec::ParametricCoordinatesToWorldCoordinates(
|
||||
pointWCoords, pcoords, shape, workletProxy);
|
||||
VTKM_TEST_ASSERT(!errorMessage.IsErrorRaised(), messageBuffer);
|
||||
bool success = false;
|
||||
Vector3 computedPCoords = vtkm::exec::WorldCoordinatesToParametricCoordinates(
|
||||
pointWCoords, wcoords, shape, workletProxy);
|
||||
pointWCoords, wcoords, shape, success, workletProxy);
|
||||
VTKM_TEST_ASSERT(!errorMessage.IsErrorRaised(), messageBuffer);
|
||||
|
||||
VTKM_TEST_ASSERT(test_equal(pcoords, computedPCoords, 0.05),
|
||||
|
@ -89,8 +89,9 @@ VTKM_EXEC_CONT inline bool Sample(const vtkm::Vec<vtkm::Vec<P, 3>, 8>& points,
|
||||
pointsVec.Append(points[i]);
|
||||
scalarVec.Append(scalars[i]);
|
||||
}
|
||||
bool success = false; // ignored
|
||||
vtkm::Vec<P, 3> pcoords = vtkm::exec::WorldCoordinatesToParametricCoordinates(
|
||||
pointsVec, sampleLocation, shapeTag, callingWorklet);
|
||||
pointsVec, sampleLocation, shapeTag, success, callingWorklet);
|
||||
P pmin, pmax;
|
||||
pmin = vtkm::Min(vtkm::Min(pcoords[0], pcoords[1]), pcoords[2]);
|
||||
pmax = vtkm::Max(vtkm::Max(pcoords[0], pcoords[1]), pcoords[2]);
|
||||
@ -112,8 +113,9 @@ VTKM_EXEC_CONT inline bool Sample(const vtkm::VecAxisAlignedPointCoordinates<3>&
|
||||
{
|
||||
|
||||
bool validSample = true;
|
||||
bool success;
|
||||
vtkm::Vec<P, 3> pcoords = vtkm::exec::WorldCoordinatesToParametricCoordinates(
|
||||
points, sampleLocation, vtkm::CellShapeTagHexahedron(), callingWorklet);
|
||||
points, sampleLocation, vtkm::CellShapeTagHexahedron(), success, callingWorklet);
|
||||
P pmin, pmax;
|
||||
pmin = vtkm::Min(vtkm::Min(pcoords[0], pcoords[1]), pcoords[2]);
|
||||
pmax = vtkm::Max(vtkm::Max(pcoords[0], pcoords[1]), pcoords[2]);
|
||||
|
@ -91,10 +91,11 @@ void TestNewtonsMethodTemplate()
|
||||
{
|
||||
std::cout << " " << initialGuess << std::endl;
|
||||
|
||||
Vector3 solution = vtkm::NewtonsMethod(
|
||||
auto result = vtkm::NewtonsMethod(
|
||||
EvaluateJacobian<T>(), EvaluateFunctions<T>(), desiredOutput, initialGuess, T(1e-6));
|
||||
|
||||
VTKM_TEST_ASSERT(test_equal(solution, expected1) || test_equal(solution, expected2),
|
||||
VTKM_TEST_ASSERT(test_equal(result.Solution, expected1) ||
|
||||
test_equal(result.Solution, expected2),
|
||||
"Newton's method did not converge to expected result.");
|
||||
}
|
||||
}
|
||||
|
@ -28,38 +28,38 @@ namespace vtkm
|
||||
namespace worklet
|
||||
{
|
||||
|
||||
//
|
||||
// Distribute multiple copies of cell data depending on cells create from original
|
||||
//
|
||||
struct DistributeCellData : public vtkm::worklet::WorkletMapField
|
||||
{
|
||||
typedef void ControlSignature(FieldIn<> inIndices, FieldOut<> outIndices);
|
||||
typedef void ExecutionSignature(_1, _2);
|
||||
|
||||
typedef vtkm::worklet::ScatterCounting ScatterType;
|
||||
|
||||
VTKM_CONT
|
||||
ScatterType GetScatter() const { return this->Scatter; }
|
||||
|
||||
template <typename CountArrayType, typename DeviceAdapter>
|
||||
VTKM_CONT DistributeCellData(const CountArrayType& countArray, DeviceAdapter device)
|
||||
: Scatter(countArray, device)
|
||||
{
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
VTKM_EXEC void operator()(T inputIndex, T& outputIndex) const
|
||||
{
|
||||
outputIndex = inputIndex;
|
||||
}
|
||||
|
||||
private:
|
||||
ScatterType Scatter;
|
||||
};
|
||||
|
||||
class Tetrahedralize
|
||||
{
|
||||
public:
|
||||
//
|
||||
// Distribute multiple copies of cell data depending on cells create from original
|
||||
//
|
||||
struct DistributeCellData : public vtkm::worklet::WorkletMapField
|
||||
{
|
||||
typedef void ControlSignature(FieldIn<> inIndices, FieldOut<> outIndices);
|
||||
typedef void ExecutionSignature(_1, _2);
|
||||
|
||||
typedef vtkm::worklet::ScatterCounting ScatterType;
|
||||
|
||||
VTKM_CONT
|
||||
ScatterType GetScatter() const { return this->Scatter; }
|
||||
|
||||
template <typename CountArrayType, typename DeviceAdapter>
|
||||
VTKM_CONT DistributeCellData(const CountArrayType& countArray, DeviceAdapter device)
|
||||
: Scatter(countArray, device)
|
||||
{
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
VTKM_EXEC void operator()(T inputIndex, T& outputIndex) const
|
||||
{
|
||||
outputIndex = inputIndex;
|
||||
}
|
||||
|
||||
private:
|
||||
ScatterType Scatter;
|
||||
};
|
||||
|
||||
Tetrahedralize()
|
||||
: OutCellsPerCell()
|
||||
{
|
||||
|
@ -28,38 +28,38 @@ namespace vtkm
|
||||
namespace worklet
|
||||
{
|
||||
|
||||
//
|
||||
// Distribute multiple copies of cell data depending on cells create from original
|
||||
//
|
||||
struct DistributeCellData : public vtkm::worklet::WorkletMapField
|
||||
{
|
||||
typedef void ControlSignature(FieldIn<> inIndices, FieldOut<> outIndices);
|
||||
typedef void ExecutionSignature(_1, _2);
|
||||
|
||||
typedef vtkm::worklet::ScatterCounting ScatterType;
|
||||
|
||||
VTKM_CONT
|
||||
ScatterType GetScatter() const { return this->Scatter; }
|
||||
|
||||
template <typename CountArrayType, typename DeviceAdapter>
|
||||
VTKM_CONT DistributeCellData(const CountArrayType& countArray, DeviceAdapter device)
|
||||
: Scatter(countArray, device)
|
||||
{
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
VTKM_EXEC void operator()(T inputIndex, T& outputIndex) const
|
||||
{
|
||||
outputIndex = inputIndex;
|
||||
}
|
||||
|
||||
private:
|
||||
ScatterType Scatter;
|
||||
};
|
||||
|
||||
class Triangulate
|
||||
{
|
||||
public:
|
||||
//
|
||||
// Distribute multiple copies of cell data depending on cells create from original
|
||||
//
|
||||
struct DistributeCellData : public vtkm::worklet::WorkletMapField
|
||||
{
|
||||
typedef void ControlSignature(FieldIn<> inIndices, FieldOut<> outIndices);
|
||||
typedef void ExecutionSignature(_1, _2);
|
||||
|
||||
typedef vtkm::worklet::ScatterCounting ScatterType;
|
||||
|
||||
VTKM_CONT
|
||||
ScatterType GetScatter() const { return this->Scatter; }
|
||||
|
||||
template <typename CountArrayType, typename DeviceAdapter>
|
||||
VTKM_CONT DistributeCellData(const CountArrayType& countArray, DeviceAdapter device)
|
||||
: Scatter(countArray, device)
|
||||
{
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
VTKM_EXEC void operator()(T inputIndex, T& outputIndex) const
|
||||
{
|
||||
outputIndex = inputIndex;
|
||||
}
|
||||
|
||||
private:
|
||||
ScatterType Scatter;
|
||||
};
|
||||
|
||||
Triangulate()
|
||||
: OutCellsPerCell()
|
||||
{
|
||||
|
Loading…
Reference in New Issue
Block a user