Filter for adding ghost zones.

This commit is contained in:
Dave Pugmire 2019-02-07 15:26:20 -05:00
parent 4e1fd2892c
commit bebbd92e32
7 changed files with 593 additions and 0 deletions

@ -47,6 +47,7 @@ set(headers
Pair.h
Range.h
RangeId.h
RangeId2.h
RangeId3.h
StaticAssert.h
Swap.h

166
vtkm/RangeId2.h Normal file

@ -0,0 +1,166 @@
//============================================================================
// 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_RangeId2_h
#define vtk_m_RangeId2_h
#include <vtkm/RangeId.h>
namespace vtkm
{
/// \brief Represent 2D integer range.
///
/// \c vtkm::RangeId2 is a helper class for representing a 2D range of integer
/// values. The typical use of this class is to express a box of indices
/// in the x, y, and z directions.
///
/// \c RangeId2 also contains several helper functions for computing and
/// maintaining the range.
///
struct RangeId2
{
vtkm::RangeId X;
vtkm::RangeId Y;
RangeId2() = default;
VTKM_EXEC_CONT
RangeId2(const vtkm::RangeId& xrange, const vtkm::RangeId& yrange)
: X(xrange)
, Y(yrange)
{
}
VTKM_EXEC_CONT
RangeId2(vtkm::Id minX, vtkm::Id maxX, vtkm::Id minY, vtkm::Id maxY)
: X(vtkm::RangeId(minX, maxX))
, Y(vtkm::RangeId(minY, maxY))
{
}
/// Initialize range with an array of 6 values in the order xmin, xmax,
/// ymin, ymax, zmin, zmax.
///
VTKM_EXEC_CONT
explicit RangeId2(const vtkm::Id range[4])
: X(vtkm::RangeId(range[0], range[1]))
, Y(vtkm::RangeId(range[2], range[3]))
{
}
/// Initialize range with the minimum and the maximum corners
///
VTKM_EXEC_CONT
RangeId2(const vtkm::Id2& min, const vtkm::Id2& max)
: X(vtkm::RangeId(min[0], max[0]))
, Y(vtkm::RangeId(min[1], max[1]))
{
}
/// \b Determine if the range is non-empty.
///
/// \c IsNonEmpty returns true if the range is non-empty.
///
VTKM_EXEC_CONT
bool IsNonEmpty() const { return (this->X.IsNonEmpty() && this->Y.IsNonEmpty()); }
/// \b Determines if an Id2 value is within the range.
///
VTKM_EXEC_CONT
bool Contains(const vtkm::Id2& val) const
{
return (this->X.Contains(val[0]) && this->Y.Contains(val[1]));
}
/// \b Returns the center of the range.
///
/// \c Center computes the middle of the range.
///
VTKM_EXEC_CONT
vtkm::Id2 Center() const { return vtkm::Id2(this->X.Center(), this->Y.Center()); }
VTKM_EXEC_CONT
vtkm::Id2 Dimensions() const { return vtkm::Id2(this->X.Length(), this->Y.Length()); }
/// \b Expand range to include a value.
///
/// This version of \c Include expands the range just enough to include the
/// given value. If the range already include this value, then
/// nothing is done.
///
template <typename T>
VTKM_EXEC_CONT void Include(const vtkm::Vec<T, 2>& point)
{
this->X.Include(point[0]);
this->Y.Include(point[1]);
}
/// \b Expand range to include other range.
///
/// This version of \c Include expands the range just enough to include
/// the other range. Essentially it is the union of the two ranges.
///
VTKM_EXEC_CONT
void Include(const vtkm::RangeId2& range)
{
this->X.Include(range.X);
this->Y.Include(range.Y);
}
/// \b Return the union of this and another range.
///
/// This is a nondestructive form of \c Include.
///
VTKM_EXEC_CONT
vtkm::RangeId2 Union(const vtkm::RangeId2& other) const
{
vtkm::RangeId2 unionRangeId2(*this);
unionRangeId2.Include(other);
return unionRangeId2;
}
/// \b Operator for union
///
VTKM_EXEC_CONT
vtkm::RangeId2 operator+(const vtkm::RangeId2& other) const { return this->Union(other); }
VTKM_EXEC_CONT
bool operator==(const vtkm::RangeId2& range) const
{
return ((this->X == range.X) && (this->Y == range.Y));
}
VTKM_EXEC_CONT
bool operator!=(const vtkm::RangeId2& range) const
{
return ((this->X != range.X) || (this->Y != range.Y));
}
};
} // namespace vtkm
/// Helper function for printing range during testing
///
static inline VTKM_CONT std::ostream& operator<<(std::ostream& stream, const vtkm::RangeId2& range)
{
return stream << "{ X:" << range.X << ", Y:" << range.Y << " }";
}
#endif //vtk_m_RangeId2_h

@ -0,0 +1,58 @@
//============================================================================
// 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 2016 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
// Copyright 2016 UT-Battelle, LLC.
// Copyright 2016 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_filter_AddGhostZone_h
#define vtk_m_filter_AddGhostZone_h
#include <vtkm/filter/FilterDataSet.h>
namespace vtkm
{
namespace filter
{
struct AddGhostZonePolicy : vtkm::filter::PolicyBase<AddGhostZonePolicy>
{
using FieldTypeList = vtkm::ListTagBase<vtkm::UInt8>;
};
class AddGhostZone : public vtkm::filter::FilterDataSet<AddGhostZone>
{
public:
VTKM_CONT
AddGhostZone();
template <typename Policy>
VTKM_CONT vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet& inData,
vtkm::filter::PolicyBase<Policy> policy);
template <typename ValueType, typename Storage, typename Policy>
VTKM_CONT bool DoMapField(vtkm::cont::DataSet& result,
const vtkm::cont::ArrayHandle<ValueType, Storage>& input,
const vtkm::filter::FieldMetadata& fieldMeta,
vtkm::filter::PolicyBase<Policy>);
private:
};
}
} // namespace vtkm::filter
#include <vtkm/filter/AddGhostZone.hxx>
#endif //vtk_m_filter_AddGhostZone_h

@ -0,0 +1,205 @@
//============================================================================
// 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 2016 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
// Copyright 2016 UT-Battelle, LLC.
// Copyright 2016 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/GhostCell.h>
#include <vtkm/RangeId.h>
#include <vtkm/RangeId2.h>
#include <vtkm/RangeId3.h>
#include <vtkm/Types.h>
#include <vtkm/cont/ArrayCopy.h>
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/ArrayHandleConstant.h>
#include <vtkm/filter/internal/CreateResult.h>
#include <vtkm/worklet/WorkletMapField.h>
namespace
{
struct TypeUInt8 : vtkm::ListTagBase<vtkm::UInt8>
{
};
class SetStructuredGhostZones1D : public vtkm::worklet::WorkletMapField
{
public:
SetStructuredGhostZones1D(const vtkm::Id& dim, const vtkm::Id& numLayers = 1)
: Dim(dim)
, NumLayers(numLayers)
, Range(NumLayers, Dim - NumLayers)
{
}
using ControlSignature = void(FieldIn, FieldOut);
using ExecutionSignature = void(_1, _2);
template <typename T>
VTKM_EXEC void operator()(const vtkm::Id& cellIndex, T& value) const
{
value = (Range.Contains(cellIndex) ? NormalCell : DuplicateCell);
}
private:
vtkm::Id Dim;
vtkm::Id NumLayers;
vtkm::RangeId Range;
static constexpr vtkm::UInt8 NormalCell =
static_cast<vtkm::UInt8>(vtkm::CellClassification::NORMAL);
static constexpr vtkm::UInt8 DuplicateCell =
static_cast<vtkm::UInt8>(vtkm::CellClassification::DUPLICATE);
};
class SetStructuredGhostZones2D : public vtkm::worklet::WorkletMapField
{
public:
SetStructuredGhostZones2D(const vtkm::Id2& dims, const vtkm::Id& numLayers = 1)
: Dims(dims)
, NumLayers(numLayers)
, Range(NumLayers, Dims[0] - NumLayers, NumLayers, Dims[1] - NumLayers)
{
}
using ControlSignature = void(FieldIn, FieldOut);
using ExecutionSignature = void(_1, _2);
template <typename T>
VTKM_EXEC void operator()(const vtkm::Id& cellIndex, T& value) const
{
vtkm::Id2 ij(cellIndex % Dims[0], cellIndex / Dims[0]);
value = (Range.Contains(ij) ? NormalCell : DuplicateCell);
}
private:
vtkm::Id2 Dims;
vtkm::Id NumLayers;
vtkm::RangeId2 Range;
static constexpr vtkm::UInt8 NormalCell =
static_cast<vtkm::UInt8>(vtkm::CellClassification::NORMAL);
static constexpr vtkm::UInt8 DuplicateCell =
static_cast<vtkm::UInt8>(vtkm::CellClassification::DUPLICATE);
};
class SetStructuredGhostZones3D : public vtkm::worklet::WorkletMapField
{
public:
SetStructuredGhostZones3D(const vtkm::Id3& dims, const vtkm::Id& numLayers = 1)
: Dims(dims)
, NumLayers(numLayers)
, Range(NumLayers,
Dims[0] - NumLayers,
NumLayers,
Dims[1] - NumLayers,
NumLayers,
Dims[2] - NumLayers)
{
}
using ControlSignature = void(FieldIn, FieldOut);
using ExecutionSignature = void(_1, _2);
template <typename T>
VTKM_EXEC void operator()(const vtkm::Id& cellIndex, T& value) const
{
vtkm::Id3 ijk(
cellIndex % Dims[0], (cellIndex / Dims[0]) % Dims[1], cellIndex / (Dims[0] * Dims[1]));
value = (Range.Contains(ijk) ? NormalCell : DuplicateCell);
}
private:
vtkm::Id3 Dims;
vtkm::Id NumLayers;
vtkm::RangeId3 Range;
static constexpr vtkm::UInt8 NormalCell =
static_cast<vtkm::UInt8>(vtkm::CellClassification::NORMAL);
static constexpr vtkm::UInt8 DuplicateCell =
static_cast<vtkm::UInt8>(vtkm::CellClassification::DUPLICATE);
};
};
namespace vtkm
{
namespace filter
{
inline VTKM_CONT AddGhostZone::AddGhostZone()
{
}
template <typename Policy>
inline VTKM_CONT vtkm::cont::DataSet AddGhostZone::DoExecute(const vtkm::cont::DataSet& input,
vtkm::filter::PolicyBase<Policy>)
{
const vtkm::cont::DynamicCellSet& cellset = input.GetCellSet(this->GetActiveCellSetIndex());
vtkm::Id numCells = cellset.GetNumberOfCells();
vtkm::cont::ArrayHandleIndex indexArray(numCells);
vtkm::cont::ArrayHandle<vtkm::UInt8> ghosts;
ghosts.Allocate(numCells);
//Structured cases are easy...
if (cellset.template IsType<vtkm::cont::CellSetStructured<1>>())
{
if (numCells <= 2)
throw vtkm::cont::ErrorFilterExecution("insufficient number of cells for AddGhostZone.");
vtkm::cont::CellSetStructured<1> cellset1d = cellset.Cast<vtkm::cont::CellSetStructured<1>>();
SetStructuredGhostZones1D structuredGhosts1D(cellset1d.GetCellDimensions());
vtkm::worklet::DispatcherMapField<SetStructuredGhostZones1D> dispatcher(structuredGhosts1D);
dispatcher.Invoke(indexArray, ghosts);
}
else if (cellset.template IsType<vtkm::cont::CellSetStructured<2>>())
{
if (numCells <= 4)
throw vtkm::cont::ErrorFilterExecution("insufficient number of cells for AddGhostZone.");
vtkm::cont::CellSetStructured<2> cellset2d = cellset.Cast<vtkm::cont::CellSetStructured<2>>();
SetStructuredGhostZones2D structuredGhosts2D(cellset2d.GetCellDimensions());
vtkm::worklet::DispatcherMapField<SetStructuredGhostZones2D> dispatcher(structuredGhosts2D);
dispatcher.Invoke(indexArray, ghosts);
}
else if (cellset.template IsType<vtkm::cont::CellSetStructured<3>>())
{
if (numCells <= 8)
throw vtkm::cont::ErrorFilterExecution("insufficient number of cells for AddGhostZone.");
vtkm::cont::CellSetStructured<3> cellset3d = cellset.Cast<vtkm::cont::CellSetStructured<3>>();
SetStructuredGhostZones3D structuredGhosts3D(cellset3d.GetCellDimensions());
vtkm::worklet::DispatcherMapField<SetStructuredGhostZones3D> dispatcher(structuredGhosts3D);
dispatcher.Invoke(indexArray, ghosts);
}
else
throw vtkm::cont::ErrorFilterExecution("Unsupported cellset type for AddGhostZone.");
vtkm::cont::DataSet output = internal::CreateResult(
input, ghosts, "vtkmGhostCells", vtkm::cont::Field::Association::CELL_SET, cellset.GetName());
return output;
}
template <typename ValueType, typename Storage, typename Policy>
inline VTKM_CONT bool AddGhostZone::DoMapField(vtkm::cont::DataSet&,
const vtkm::cont::ArrayHandle<ValueType, Storage>&,
const vtkm::filter::FieldMetadata&,
vtkm::filter::PolicyBase<Policy>)
{
return true;
}
}
}

@ -19,6 +19,7 @@
##============================================================================
set(headers
AddGhostZone.h
CellAverage.h
CellMeasures.h
CleanGrid.h
@ -79,6 +80,7 @@ set(headers
)
set(header_template_sources
AddGhostZone.hxx
CellAverage.hxx
CellMeasures.hxx
CleanGrid.hxx

@ -19,6 +19,7 @@
##============================================================================
set(unit_tests
UnitTestAddGhostZone.cxx
UnitTestCellAverageFilter.cxx
UnitTestCellMeasuresFilter.cxx
UnitTestCleanGrid.cxx

@ -0,0 +1,160 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//
// Copyright 2014 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
// Copyright 2014 UT-Battelle, LLC.
// Copyright 2014 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/DataSet.h>
#include <vtkm/cont/DataSetBuilderExplicit.h>
#include <vtkm/cont/DataSetBuilderRectilinear.h>
#include <vtkm/cont/DataSetBuilderUniform.h>
#include <vtkm/cont/DataSetFieldAdd.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/filter/AddGhostZone.h>
namespace
{
static vtkm::cont::DataSet MakeUniform(vtkm::Id numI, vtkm::Id numJ, vtkm::Id numK)
{
vtkm::cont::DataSetBuilderUniform dsb;
vtkm::cont::DataSet ds;
if (numJ == 0 && numK == 0)
ds = dsb.Create(numI + 1);
else if (numK == 0)
ds = dsb.Create(vtkm::Id2(numI + 1, numJ + 1));
else
ds = dsb.Create(vtkm::Id3(numI + 1, numJ + 1, numK + 1));
return ds;
}
static vtkm::cont::DataSet MakeRectilinear(vtkm::Id numI, vtkm::Id numJ, vtkm::Id numK)
{
vtkm::cont::DataSetBuilderRectilinear dsb;
vtkm::cont::DataSet ds;
std::size_t nx(static_cast<std::size_t>(numI + 1));
std::size_t ny(static_cast<std::size_t>(numJ + 1));
std::vector<float> x(nx), y(ny);
for (std::size_t i = 0; i < nx; i++)
x[i] = static_cast<float>(i);
for (std::size_t i = 0; i < ny; i++)
y[i] = static_cast<float>(i);
if (numK == 0)
ds = dsb.Create(x, y);
else
{
std::size_t nz(static_cast<std::size_t>(numK + 1));
std::vector<float> z(nz);
for (std::size_t i = 0; i < nz; i++)
z[i] = static_cast<float>(i);
ds = dsb.Create(x, y, z);
}
return ds;
}
void TestStructured()
{
std::cout << "Testing ghost cells for structured datasets." << std::endl;
// specify some 2d tests: {numI, numJ, numK, numGhostLayers}.
std::vector<std::vector<vtkm::Id>> tests2D = { { 8, 4, 0, 1 }, { 5, 5, 0, 1 }, { 10, 10, 0, 1 },
{ 10, 5, 0, 1 }, { 5, 10, 0, 1 }, { 20, 10, 0, 1 },
{ 10, 20, 0, 1 } };
std::vector<std::vector<vtkm::Id>> tests3D = { { 8, 8, 10, 1 }, { 5, 5, 5, 1 },
{ 10, 10, 10, 1 }, { 10, 5, 10, 1 },
{ 5, 10, 10, 1 }, { 20, 10, 10, 1 },
{ 10, 20, 10, 1 } };
std::vector<std::vector<vtkm::Id>> tests1D = {
{ 8, 0, 0, 1 }, { 5, 0, 0, 1 }, { 10, 0, 0, 1 }, { 20, 0, 0, 1 }
};
std::vector<std::vector<vtkm::Id>> tests;
tests.insert(tests.end(), tests1D.begin(), tests1D.end());
tests.insert(tests.end(), tests2D.begin(), tests2D.end());
tests.insert(tests.end(), tests3D.begin(), tests3D.end());
for (auto& t : tests)
{
vtkm::Id nx = t[0], ny = t[1], nz = t[2];
vtkm::Id nghost = t[3];
for (vtkm::Id layer = 1; layer <= nghost; layer++)
{
vtkm::cont::DataSet ds;
std::vector<std::string> dsTypes = { "uniform", "rectilinear" };
for (auto& dsType : dsTypes)
{
if (dsType == "uniform")
ds = MakeUniform(nx, ny, nz);
else if (dsType == "rectilinear")
ds = MakeRectilinear(nx, ny, nz);
vtkm::filter::AddGhostZone addGhost;
auto output = addGhost.Execute(ds, vtkm::filter::AddGhostZonePolicy());
//Validate the output.
VTKM_TEST_ASSERT(output.GetNumberOfCellSets() == 1, "Wrong number of cell sets in output");
VTKM_TEST_ASSERT(
output.HasField("vtkmGhostCells", vtkm::cont::Field::Association::CELL_SET),
"Ghost cells array not found in output");
vtkm::Id numCells = output.GetCellSet(0).GetNumberOfCells();
auto fieldArray =
output.GetField("vtkmGhostCells", vtkm::cont::Field::Association::CELL_SET).GetData();
VTKM_TEST_ASSERT(fieldArray.GetNumberOfValues() == numCells,
"Wrong number of values in ghost cell array");
//Check the number of normal cells.
vtkm::cont::ArrayHandle<vtkm::UInt8> ghostArray;
fieldArray.CopyTo(ghostArray);
vtkm::Id numNormalCells = 0;
auto portal = ghostArray.GetPortalConstControl();
const vtkm::UInt8 normalCell = static_cast<vtkm::UInt8>(vtkm::CellClassification::NORMAL);
for (vtkm::Id i = 0; i < numCells; i++)
if (portal.Get(i) == normalCell)
numNormalCells++;
vtkm::Id requiredNumNormalCells = (nx - 2 * layer);
if (ny > 0)
requiredNumNormalCells *= (ny - 2 * layer);
if (nz > 0)
requiredNumNormalCells *= (nz - 2 * layer);
VTKM_TEST_ASSERT(requiredNumNormalCells == numNormalCells, "Wrong number of normal cells");
}
}
}
}
void TestAddGhostZone()
{
TestStructured();
}
}
int UnitTestAddGhostZone(int argc, char* argv[])
{
return vtkm::cont::testing::Testing::Run(TestAddGhostZone, argc, argv);
}