Shuenhoy/vtk-m
1
0
Fork 0
mirror of https://gitlab.kitware.com/vtk/vtk-m synced 2024-09-16 17:22:55 +00:00
Code Issues 2 Actions Packages Projects Releases Wiki Activity

Add a split sharp edge worklet and filter

Browse Source 
It's a filter that Split sharp manifold edges where the feature angle
between the adjacent surfaces are larger than the threshold value.
When an edge is split, it would add a new point to the coordinates
and update the connectivity of an adjacent surface.
Ex. there are two adjacent triangles(0,1,2) and (2,1,3). Edge (1,2) needs
to be split. Two new points 4(duplication of point 1) an 5(duplication of point 2)
would be added and the later triangle's connectivity would be changed
to (5,4,3).
By default, all old point's fields would be copied to the new point.
Use with caution.
This commit is contained in:
Haocheng LIU 2018-09-07 16:49:51 -04:00
parent 6543942568
commit 246a58309c
14 changed files with 1314 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
docs/changelog/split-sharp-edges-filter.md Normal file
View File

@ -0,0 +1,12 @@
Add a split sharp edge filter
It's a filter that splits sharp manifold edges where the feature angle
between the adjacent surfaces are larger than the threshold value.
When an edge is split, it would add a new point to the coordinates
and update the connectivity of an adjacent surface.
Ex. There are two adjacent triangles(0,1,2) and (2,1,3). Edge (1,2) needs
to be split. Two new points 4(duplication of point 1) an 5(duplication of point 2)
would be added and the later triangle's connectivity would be changed
to (5,4,3).
By default, all old point's fields would be copied to the new point.
Use with caution.

70
vtkm/Bitset.h Normal file
View File

@ -0,0 +1,70 @@
//============================================================================
// 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_Bitset_h
#define vtk_m_Bitset_h
#include <assert.h>
#include <limits>
#include <vtkm/Types.h>
#include <vtkm/internal/ExportMacros.h>
namespace vtkm
{
/// \brief A bitmap to serve different needs.
/// Ex. Editing particular bits in a byte(s), checkint if particular bit values
/// are present or not. Once Cuda supports std::bitset, we should use the
/// standard one if possible
template <typename MaskType>
struct Bitset
{
VTKM_EXEC_CONT Bitset()
: Mask(0)
{
}
VTKM_EXEC_CONT void set(vtkm::Id bitIndex)
{
this->Mask = this->Mask | (static_cast<MaskType>(1) << bitIndex);
}
VTKM_EXEC_CONT void reset(vtkm::Id bitIndex)
{
this->Mask = this->Mask & ~(static_cast<MaskType>(1) << bitIndex);
}
VTKM_EXEC_CONT void toggle(vtkm::Id bitIndex)
{
this->Mask = this->Mask ^ (static_cast<MaskType>(0) << bitIndex);
}
VTKM_EXEC_CONT bool test(vtkm::Id bitIndex)
{
return ((this->Mask & (static_cast<MaskType>(1) << bitIndex)) != 0);
}
private:
MaskType Mask;
};
} // namespace vtkm
#endif //vtk_m_Bitset_h

1
vtkm/CMakeLists.txt
View File

@ -31,6 +31,7 @@ set(headers
BaseComponent.h
BinaryPredicates.h
BinaryOperators.h
Bitset.h
Bounds.h
CellShape.h
CellTraits.h

2
vtkm/filter/CMakeLists.txt
View File

@ -59,6 +59,7 @@ set(headers
PolicyBase.h
PolicyDefault.h
Probe.h
SplitSharpEdges.h
Streamline.h
SurfaceNormals.h
Tetrahedralize.h
@ -107,6 +108,7 @@ set(header_template_sources
PointElevation.hxx
PointTransform.hxx
Probe.hxx
SplitSharpEdges.hxx
Streamline.hxx
SurfaceNormals.hxx
Tetrahedralize.hxx

2
vtkm/filter/CoordinateSystemTransform.h
View File

@ -84,7 +84,7 @@ template <>
class FilterTraits<SphericalCoordinateTransform>
{
public:
//Point Elevation can only convert Float and Double Vec3 arrays
//CoordinateSystemTransformation can only convert Float and Double Vec3 arrays
using InputFieldTypeList = vtkm::TypeListTagFieldVec3;
};
}

2
vtkm/filter/PointTransform.h
View File

@ -75,7 +75,7 @@ template <typename S>
class FilterTraits<PointTransform<S>>
{
public:
//Point Elevation can only convert Float and Double Vec3 arrays
//PointTransformation can only convert Float and Double Vec3 arrays
using InputFieldTypeList = vtkm::TypeListTagFieldVec3;
};
}

88
vtkm/filter/SplitSharpEdges.h Normal file
View File

@ -0,0 +1,88 @@
//============================================================================
// 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.
//============================================================================
#ifndef vtk_m_filter_SplitSharpEdges_h
#define vtk_m_filter_SplitSharpEdges_h
#include <vtkm/filter/FilterDataSetWithField.h>
#include <vtkm/worklet/SplitSharpEdges.h>
namespace vtkm
{
namespace filter
{
/// \brief Split sharp manifold edges where the feature angle between the
/// adjacent surfaces are larger than the threshold value
///
/// Split sharp manifold edges where the feature angle between the adjacent
/// surfaces are larger than the threshold value. When an edge is split, it
/// would add a new point to the coordinates and update the connectivity of
/// an adjacent surface.
/// Ex. there are two adjacent triangles(0,1,2) and (2,1,3). Edge (1,2) needs
/// to be split. Two new points 4(duplication of point 1) an 5(duplication of point 2)
/// would be added and the later triangle's connectivity would be changed
/// to (5,4,3).
/// By default, all old point's fields would be copied to the new point.
/// Use with caution.
class SplitSharpEdges : public vtkm::filter::FilterDataSetWithField<SplitSharpEdges>
{
public:
VTKM_CONT
SplitSharpEdges();
VTKM_CONT
void SetFeatureAngle(vtkm::FloatDefault value) { this->FeatureAngle = value; }
VTKM_CONT
vtkm::FloatDefault GetFeatureAngle() const { return this->FeatureAngle; }
template <typename T, typename StorageType, typename DerivedPolicy, typename DeviceAdapter>
VTKM_CONT vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet& input,
const vtkm::cont::ArrayHandle<T, StorageType>& field,
const vtkm::filter::FieldMetadata& fieldMeta,
const vtkm::filter::PolicyBase<DerivedPolicy>& policy,
DeviceAdapter tag);
//Map a new field onto the resulting dataset after running the filter
template <typename T, typename StorageType, typename DerivedPolicy, typename DeviceAdapter>
VTKM_CONT bool DoMapField(vtkm::cont::DataSet& result,
const vtkm::cont::ArrayHandle<T, StorageType>& input,
const vtkm::filter::FieldMetadata& fieldMeta,
const vtkm::filter::PolicyBase<DerivedPolicy>& policy,
DeviceAdapter tag);
private:
vtkm::FloatDefault FeatureAngle;
vtkm::worklet::SplitSharpEdges Worklet;
};
template <>
class FilterTraits<SplitSharpEdges>
{ // SplitSharpEdges filter needs cell normals to decide split.
public:
using InputFieldTypeList = vtkm::TypeListTagFieldVec3;
};
}
} // namespace vtkm::filter
#include <vtkm/filter/SplitSharpEdges.hxx>
#endif // vtk_m_filter_SplitSharpEdges_h

96
vtkm/filter/SplitSharpEdges.hxx Normal file
View File

@ -0,0 +1,96 @@
//============================================================================
// 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/ArrayHandlePermutation.h>
#include <vtkm/cont/CellSetExplicit.h>
#include <vtkm/cont/CoordinateSystem.h>
#include <vtkm/cont/DynamicCellSet.h>
#include <vtkm/worklet/DispatcherMapTopology.h>
namespace vtkm
{
namespace filter
{
//-----------------------------------------------------------------------------
inline VTKM_CONT SplitSharpEdges::SplitSharpEdges()
: vtkm::filter::FilterDataSetWithField<SplitSharpEdges>()
, FeatureAngle(30.0)
{
}
//-----------------------------------------------------------------------------
template <typename T, typename StorageType, typename DerivedPolicy, typename DeviceAdapter>
inline VTKM_CONT vtkm::cont::DataSet SplitSharpEdges::DoExecute(
const vtkm::cont::DataSet& input,
const vtkm::cont::ArrayHandle<T, StorageType>& field,
const vtkm::filter::FieldMetadata& vtkmNotUsed(fieldMeta),
const vtkm::filter::PolicyBase<DerivedPolicy>& policy,
DeviceAdapter)
{
// Get the cells and coordinates of the dataset
const vtkm::cont::DynamicCellSet& cells = input.GetCellSet(this->GetActiveCellSetIndex());
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::FloatDefault, 3>> newCoords;
vtkm::cont::CellSetExplicit<> newCellset;
this->Worklet.Run(vtkm::filter::ApplyPolicy(cells, policy),
this->FeatureAngle,
field,
input.GetCoordinateSystem().GetData(),
newCoords,
newCellset,
DeviceAdapter());
vtkm::cont::DataSet output;
output.AddCellSet(newCellset);
output.AddCoordinateSystem(
vtkm::cont::CoordinateSystem(input.GetCoordinateSystem().GetName(), newCoords));
return output;
}
//-----------------------------------------------------------------------------
template <typename T, typename StorageType, typename DerivedPolicy, typename DeviceAdapter>
inline VTKM_CONT bool SplitSharpEdges::DoMapField(
vtkm::cont::DataSet& result,
const vtkm::cont::ArrayHandle<T, StorageType>& input,
const vtkm::filter::FieldMetadata& fieldMeta,
const vtkm::filter::PolicyBase<DerivedPolicy>&,
DeviceAdapter device)
{
if (fieldMeta.IsPointField())
{
// We copy the input handle to the result dataset, reusing the metadata
vtkm::cont::ArrayHandle<T> out = this->Worklet.ProcessPointField(input, device);
result.AddField(fieldMeta.AsField(out));
return true;
}
else if (fieldMeta.IsCellField())
{
result.AddField(fieldMeta.AsField(input));
return true;
}
else
{
return false;
}
}
}
}

1
vtkm/filter/testing/CMakeLists.txt
View File

@ -52,6 +52,7 @@ set(unit_tests
UnitTestPointElevationFilter.cxx
UnitTestPointTransform.cxx
UnitTestProbe.cxx
UnitTestSplitSharpEdgesFilter.cxx
UnitTestStreamlineFilter.cxx
UnitTestSurfaceNormalsFilter.cxx
UnitTestTetrahedralizeFilter.cxx

225
vtkm/filter/testing/UnitTestSplitSharpEdgesFilter.cxx Normal file
View File

@ -0,0 +1,225 @@
//============================================================================
// 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/filter/SplitSharpEdges.h>
#include <vtkm/filter/SurfaceNormals.h>
#include <vtkm/cont/testing/MakeTestDataSet.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/worklet/DispatcherMapTopology.h>
namespace
{
using NormalsArrayHandle = vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::FloatDefault, 3>>;
using DeviceAdapter = VTKM_DEFAULT_DEVICE_ADAPTER_TAG;
const vtkm::Vec<vtkm::FloatDefault, 3> expectedCoords[24] = {
{ 0.0, 0.0, 0.0 }, { 1.0, 0.0, 0.0 }, { 1.0, 0.0, 1.0 }, { 0.0, 0.0, 1.0 }, { 0.0, 1.0, 0.0 },
{ 1.0, 1.0, 0.0 }, { 1.0, 1.0, 1.0 }, { 0.0, 1.0, 1.0 }, { 0.0, 0.0, 0.0 }, { 0.0, 0.0, 0.0 },
{ 1.0, 0.0, 0.0 }, { 1.0, 0.0, 0.0 }, { 1.0, 0.0, 1.0 }, { 1.0, 0.0, 1.0 }, { 0.0, 0.0, 1.0 },
{ 0.0, 0.0, 1.0 }, { 0.0, 1.0, 0.0 }, { 0.0, 1.0, 0.0 }, { 1.0, 1.0, 0.0 }, { 1.0, 1.0, 0.0 },
{ 1.0, 1.0, 1.0 }, { 1.0, 1.0, 1.0 }, { 0.0, 1.0, 1.0 }, { 0.0, 1.0, 1.0 }
};
const std::vector<vtkm::Id> expectedConnectivityArray91{ 0, 1, 5, 4, 1, 2, 6, 5, 2, 3, 7, 6,
3, 0, 4, 7, 4, 5, 6, 7, 0, 3, 2, 1 };
const std::vector<vtkm::FloatDefault> expectedPointvar{ 10.1f, 20.1f, 30.2f, 40.2f, 50.3f, 60.3f,
70.3f, 80.3f, 10.1f, 10.1f, 20.1f, 20.1f,
30.2f, 30.2f, 40.2f, 40.2f, 50.3f, 50.3f,
60.3f, 60.3f, 70.3f, 70.3f, 80.3f, 80.3f };
vtkm::cont::DataSet Make3DExplicitSimpleCube()
{
vtkm::cont::DataSet dataSet;
vtkm::cont::DataSetBuilderExplicit dsb;
const int nVerts = 8;
const int nCells = 6;
using CoordType = vtkm::Vec<vtkm::FloatDefault, 3>;
std::vector<CoordType> coords = {
CoordType(0, 0, 0), // 0
CoordType(1, 0, 0), // 1
CoordType(1, 0, 1), // 2
CoordType(0, 0, 1), // 3
CoordType(0, 1, 0), // 4
CoordType(1, 1, 0), // 5
CoordType(1, 1, 1), // 6
CoordType(0, 1, 1) // 7
};
//Connectivity
std::vector<vtkm::UInt8> shapes;
std::vector<vtkm::IdComponent> numIndices;
for (size_t i = 0; i < 6; i++)
{
shapes.push_back(vtkm::CELL_SHAPE_QUAD);
numIndices.push_back(4);
}
std::vector<vtkm::Id> conn;
// Down face
conn.push_back(0);
conn.push_back(1);
conn.push_back(5);
conn.push_back(4);
// Right face
conn.push_back(1);
conn.push_back(2);
conn.push_back(6);
conn.push_back(5);
// Top face
conn.push_back(2);
conn.push_back(3);
conn.push_back(7);
conn.push_back(6);
// Left face
conn.push_back(3);
conn.push_back(0);
conn.push_back(4);
conn.push_back(7);
// Front face
conn.push_back(4);
conn.push_back(5);
conn.push_back(6);
conn.push_back(7);
// Back face
conn.push_back(0);
conn.push_back(3);
conn.push_back(2);
conn.push_back(1);
//Create the dataset.
dataSet = dsb.Create(coords, shapes, numIndices, conn, "coordinates", "cells");
vtkm::FloatDefault vars[nVerts] = { 10.1f, 20.1f, 30.2f, 40.2f, 50.3f, 60.3f, 70.3f, 80.3f };
vtkm::FloatDefault cellvar[nCells] = { 100.1f, 200.2f, 300.3f, 400.4f, 500.5f, 600.6f };
vtkm::cont::DataSetFieldAdd dsf;
dsf.AddPointField(dataSet, "pointvar", vars, nVerts);
dsf.AddCellField(dataSet, "cellvar", cellvar, nCells, "cells");
return dataSet;
}
void TestSplitSharpEdgesFilterSplitEveryEdge(vtkm::cont::DataSet& simpleCubeWithSN,
vtkm::filter::SplitSharpEdges& splitSharpEdgesFilter)
{
// Split every edge
vtkm::FloatDefault featureAngle = 89.0;
splitSharpEdgesFilter.SetFeatureAngle(featureAngle);
splitSharpEdgesFilter.SetActiveField("Normals", vtkm::cont::Field::Association::CELL_SET);
vtkm::cont::DataSet result = splitSharpEdgesFilter.Execute(simpleCubeWithSN);
auto newCoords = result.GetCoordinateSystem().GetData();
auto newCoordsP = newCoords.GetPortalConstControl();
vtkm::cont::ArrayHandle<vtkm::FloatDefault> newPointvarField;
result.GetField("pointvar").GetData().CopyTo(newPointvarField);
for (vtkm::IdComponent i = 0; i < newCoords.GetNumberOfValues(); i++)
{
VTKM_TEST_ASSERT(test_equal(newCoordsP.Get(i)[0], expectedCoords[i][0]),
"result value does not match expected value");
VTKM_TEST_ASSERT(test_equal(newCoordsP.Get(i)[1], expectedCoords[i][1]),
"result value does not match expected value");
VTKM_TEST_ASSERT(test_equal(newCoordsP.Get(i)[2], expectedCoords[i][2]),
"result value does not match expected value");
}
auto newPointvarFieldPortal = newPointvarField.GetPortalConstControl();
for (vtkm::IdComponent i = 0; i < newPointvarField.GetNumberOfValues(); i++)
{
VTKM_TEST_ASSERT(test_equal(newPointvarFieldPortal.Get(static_cast<vtkm::Id>(i)),
expectedPointvar[static_cast<unsigned long>(i)]),
"point field array result does not match expected value");
}
}
void TestSplitSharpEdgesFilterNoSplit(vtkm::cont::DataSet& simpleCubeWithSN,
vtkm::filter::SplitSharpEdges& splitSharpEdgesFilter)
{
// Do nothing
vtkm::FloatDefault featureAngle = 91.0;
splitSharpEdgesFilter.SetFeatureAngle(featureAngle);
splitSharpEdgesFilter.SetActiveField("Normals", vtkm::cont::Field::Association::CELL_SET);
vtkm::cont::DataSet result = splitSharpEdgesFilter.Execute(simpleCubeWithSN);
auto newCoords = result.GetCoordinateSystem().GetData();
vtkm::cont::CellSetExplicit<>& newCellset =
result.GetCellSet().Cast<vtkm::cont::CellSetExplicit<>>();
auto newCoordsP = newCoords.GetPortalConstControl();
vtkm::cont::ArrayHandle<vtkm::FloatDefault> newPointvarField;
result.GetField("pointvar").GetData().CopyTo(newPointvarField);
for (vtkm::IdComponent i = 0; i < newCoords.GetNumberOfValues(); i++)
{
VTKM_TEST_ASSERT(test_equal(newCoordsP.Get(i)[0], expectedCoords[i][0]),
"result value does not match expected value");
VTKM_TEST_ASSERT(test_equal(newCoordsP.Get(i)[1], expectedCoords[i][1]),
"result value does not match expected value");
VTKM_TEST_ASSERT(test_equal(newCoordsP.Get(i)[2], expectedCoords[i][2]),
"result value does not match expected value");
}
const auto& connectivityArray = newCellset.GetConnectivityArray(vtkm::TopologyElementTagPoint(),
vtkm::TopologyElementTagCell());
auto connectivityArrayPortal = connectivityArray.GetPortalConstControl();
for (vtkm::IdComponent i = 0; i < connectivityArray.GetNumberOfValues(); i++)
{
VTKM_TEST_ASSERT(connectivityArrayPortal.Get(static_cast<vtkm::Id>(i)) ==
expectedConnectivityArray91[static_cast<unsigned long>(i)],
"connectivity array result does not match expected value");
}
auto newPointvarFieldPortal = newPointvarField.GetPortalConstControl();
for (vtkm::IdComponent i = 0; i < newPointvarField.GetNumberOfValues(); i++)
{
VTKM_TEST_ASSERT(test_equal(newPointvarFieldPortal.Get(static_cast<vtkm::Id>(i)),
expectedPointvar[static_cast<unsigned long>(i)]),
"point field array result does not match expected value");
}
}
void TestSplitSharpEdgesFilter()
{
vtkm::cont::DataSet simpleCube = Make3DExplicitSimpleCube();
// Generate surface normal field
vtkm::filter::SurfaceNormals surfaceNormalsFilter;
surfaceNormalsFilter.SetGenerateCellNormals(true);
vtkm::cont::DataSet simpleCubeWithSN = surfaceNormalsFilter.Execute(simpleCube);
VTKM_TEST_ASSERT(simpleCubeWithSN.HasField("Normals", vtkm::cont::Field::Association::CELL_SET),
"Cell normals missing.");
VTKM_TEST_ASSERT(simpleCubeWithSN.HasField("pointvar", vtkm::cont::Field::Association::POINTS),
"point field pointvar missing.");
vtkm::filter::SplitSharpEdges splitSharpEdgesFilter;
TestSplitSharpEdgesFilterSplitEveryEdge(simpleCubeWithSN, splitSharpEdgesFilter);
TestSplitSharpEdgesFilterNoSplit(simpleCubeWithSN, splitSharpEdgesFilter);
}
} // anonymous namespace
int UnitTestSplitSharpEdgesFilter(int, char* [])
{
return vtkm::cont::testing::Testing::Run(TestSplitSharpEdgesFilter);
}

1
vtkm/worklet/CMakeLists.txt
View File

@ -67,6 +67,7 @@ set(headers
ScatterIdentity.h
ScatterPermutation.h
ScatterUniform.h
SplitSharpEdges.h
StableSortIndices.h
StreamLineUniformGrid.h
SurfaceNormals.h

569
vtkm/worklet/SplitSharpEdges.h Normal file
View File

@ -0,0 +1,569 @@
//============================================================================
// 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_worklet_SplitSharpEdges_h
#define vtk_m_worklet_SplitSharpEdges_h
#include <map>
#include <vtkm/worklet/CellDeepCopy.h>
#include <vtkm/worklet/DispatcherMapTopology.h>
#include <vtkm/worklet/WorkletMapTopology.h>
#include <vtkm/Bitset.h>
#include <vtkm/CellTraits.h>
#include <vtkm/TypeTraits.h>
#include <vtkm/VectorAnalysis.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/exec/CellEdge.h>
namespace vtkm
{
namespace worklet
{
namespace internal
{
// Given a cell and a point on the cell, find the two edges that are
// associated with this point in canonical index
template <typename PointFromCellSetType>
VTKM_EXEC void FindRelatedEdges(const vtkm::Id& pointIndex,
const vtkm::Id& cellIndexG,
const PointFromCellSetType& pFromCellSet,
vtkm::Id2& edge0G,
vtkm::Id2& edge1G,
const vtkm::exec::FunctorBase& worklet)
{
typename PointFromCellSetType::CellShapeTag cellShape = pFromCellSet.GetCellShape(cellIndexG);
typename PointFromCellSetType::IndicesType cellConnections = pFromCellSet.GetIndices(cellIndexG);
vtkm::IdComponent numPointsInCell = pFromCellSet.GetNumberOfIndices(cellIndexG);
vtkm::IdComponent numEdges =
vtkm::exec::CellEdgeNumberOfEdges(numPointsInCell, cellShape, worklet);
vtkm::IdComponent edgeIndex = -1;
// Find the two edges with the pointIndex
while (true)
{
++edgeIndex;
if (edgeIndex >= numEdges)
{
worklet.RaiseError("Bad cell. Could not find two incident edges.");
return;
}
vtkm::Id2 canonicalEdgeId(cellConnections[vtkm::exec::CellEdgeLocalIndex(
numPointsInCell, 0, edgeIndex, cellShape, worklet)],
cellConnections[vtkm::exec::CellEdgeLocalIndex(
numPointsInCell, 1, edgeIndex, cellShape, worklet)]);
if (canonicalEdgeId[0] == pointIndex || canonicalEdgeId[1] == pointIndex)
{ // Assign value to edge0 first
if ((edge0G[0] == -1) && (edge0G[1] == -1))
{
edge0G = canonicalEdgeId;
}
else
{
edge1G = canonicalEdgeId;
break;
}
}
}
}
// TODO: We should replace this expensive lookup with a WholeCellSetIn<Edge, Cell> map.
// Given an edge on a cell, it would find the neighboring
// cell of this edge in local index. If it's a non manifold edge, -1 would be returned.
template <typename PointFromCellSetType, typename IncidentCellVecType>
VTKM_EXEC int FindNeighborCellInLocalIndex(const vtkm::Id2& eOI,
const PointFromCellSetType& pFromCellSet,
const IncidentCellVecType& incidentCells,
const vtkm::Id currentCellLocalIndex,
const vtkm::exec::FunctorBase& worklet)
{
int neighboringCellIndex = -1;
vtkm::IdComponent numberOfIncidentCells = incidentCells.GetNumberOfComponents();
size_t neighboringCellsCount = 0;
for (vtkm::IdComponent incidentCellIndex = 0; incidentCellIndex < numberOfIncidentCells;
incidentCellIndex++)
{
if (currentCellLocalIndex == incidentCellIndex)
{
continue; // No need to check the current interested cell
}
vtkm::Id cellIndexG = incidentCells[incidentCellIndex]; // Global cell index
typename PointFromCellSetType::CellShapeTag cellShape = pFromCellSet.GetCellShape(cellIndexG);
typename PointFromCellSetType::IndicesType cellConnections =
pFromCellSet.GetIndices(cellIndexG);
vtkm::IdComponent numPointsInCell = pFromCellSet.GetNumberOfIndices(cellIndexG);
vtkm::IdComponent numEdges =
vtkm::exec::CellEdgeNumberOfEdges(numPointsInCell, cellShape, worklet);
vtkm::IdComponent edgeIndex = -1;
// Check if this cell has edge of interest
while (true)
{
++edgeIndex;
if (edgeIndex >= numEdges)
{
break;
}
vtkm::Id2 canonicalEdgeId(cellConnections[vtkm::exec::CellEdgeLocalIndex(
numPointsInCell, 0, edgeIndex, cellShape, worklet)],
cellConnections[vtkm::exec::CellEdgeLocalIndex(
numPointsInCell, 1, edgeIndex, cellShape, worklet)]);
if ((canonicalEdgeId[0] == eOI[0] && canonicalEdgeId[1] == eOI[1]) ||
(canonicalEdgeId[0] == eOI[1] && canonicalEdgeId[1] == eOI[0]))
{
neighboringCellIndex = incidentCellIndex;
neighboringCellsCount++;
break;
}
}
}
return neighboringCellIndex;
}
} // internal namespace
// Split sharp manifold edges where the feature angle between the
// adjacent surfaces are larger than the threshold value
class SplitSharpEdges
{
public:
// This worklet would calculate the needed space for splitting sharp edges.
// For each point, it would have two values as numberOfNewPoint(how many
// times this point needs to be duplicated) and numberOfCellsNeedsUpdate
// (how many neighboring cells need to update connectivity).
// For example, Given a unit cube and feature angle
// as 89 degree, each point would be duplicated twice and there are two cells
// need connectivity update. There is no guarantee on which cell would get which
// new point.
class ClassifyPoint : public vtkm::worklet::WorkletMapCellToPoint
{
public:
ClassifyPoint(vtkm::FloatDefault cosfeatureAngle)
: CosFeatureAngle(cosfeatureAngle)
{
}
using ControlSignature = void(CellSetIn intputCells,
WholeCellSetIn<Point, Cell>, // Query points from cell
FieldInCell<Vec3> faceNormals,
FieldOutPoint<IdType> newPointNum,
FieldOutPoint<IdType> cellNum);
using ExecutionSignature = void(CellIndices incidentCells,
InputIndex pointIndex,
_2 pFromCellSet,
_3 faceNormals,
_4 newPointNum,
_5 cellNum);
using InputDomain = _1;
template <typename IncidentCellVecType,
typename PointFromCellSetType,
typename FaceNormalVecType>
VTKM_EXEC void operator()(const IncidentCellVecType& incidentCells,
vtkm::Id pointIndex,
const PointFromCellSetType& pFromCellSet,
const FaceNormalVecType& faceNormals,
vtkm::Id& newPointNum,
vtkm::Id& cellNum) const
{
vtkm::IdComponent numberOfIncidentCells = incidentCells.GetNumberOfComponents();
VTKM_ASSERT(numberOfIncidentCells < 64);
if (numberOfIncidentCells <= 1)
{
return; // Not enought cells to compare
}
// Initialize a global cell mask to avoid confusion. globalCellIndex->status
// 0 means not visited yet 1 means visited.
vtkm::Bitset<vtkm::UInt64> visitedCells;
vtkm::Id visitedCellsRegionIndex[64] = { 0 };
// Reallocate memory for visitedCellsGroup if needed
int regionIndex = 0;
// Loop through each cell
for (vtkm::IdComponent incidentCellIndex = 0; incidentCellIndex < numberOfIncidentCells;
incidentCellIndex++)
{
vtkm::Id cellIndexG = incidentCells[incidentCellIndex]; // cell index in global order
// If not visited
if (!visitedCells.test(incidentCellIndex))
{
// Mark the cell and track the region
visitedCells.set(incidentCellIndex);
visitedCellsRegionIndex[incidentCellIndex] = regionIndex;
// Find two edges containing the current point in canonial index
vtkm::Id2 edge0G(-1, -1), edge1G(-1, -1);
internal::FindRelatedEdges(pointIndex, cellIndexG, pFromCellSet, edge0G, edge1G, *this);
// Grow the area along each edge
for (size_t i = 0; i < 2; i++)
{ // Reset these two values for each grow operation
vtkm::Id2 currentEdgeG = i == 0 ? edge0G : edge1G;
vtkm::IdComponent currentTestingCellIndex = incidentCellIndex;
while (currentTestingCellIndex >= 0)
{
// Find the neighbor cell of the current cell edge in local index
int neighboringCellIndexQuery = internal::FindNeighborCellInLocalIndex(
currentEdgeG, pFromCellSet, incidentCells, currentTestingCellIndex, *this);
// The edge should be manifold and the neighboring cell should
// have not been visited
if (neighboringCellIndexQuery != -1 && !visitedCells.test(neighboringCellIndexQuery))
{
vtkm::IdComponent neighborCellIndex =
static_cast<vtkm::IdComponent>(neighboringCellIndexQuery);
auto thisNormal = faceNormals[currentTestingCellIndex];
//neighborNormal
auto neighborNormal = faceNormals[neighborCellIndex];
// Try to grow the area
if (vtkm::dot(thisNormal, neighborNormal) > this->CosFeatureAngle)
{ // No need to split.
visitedCells.set(neighborCellIndex);
visitedCellsRegionIndex[neighborCellIndex] = regionIndex;
// Move to examine next cell
currentTestingCellIndex = neighborCellIndex;
vtkm::Id2 neighborCellEdge0G(-1, -1), neighborCellEdge1G(-1, -1);
internal::FindRelatedEdges(pointIndex,
incidentCells[currentTestingCellIndex],
pFromCellSet,
neighborCellEdge0G,
neighborCellEdge1G,
*this);
// Update currentEdgeG
if ((currentEdgeG == neighborCellEdge0G) ||
currentEdgeG == vtkm::Id2(neighborCellEdge0G[1], neighborCellEdge0G[0]))
{
currentEdgeG = neighborCellEdge1G;
}
else
{
currentEdgeG = neighborCellEdge0G;
}
}
else
{
currentTestingCellIndex = -1;
}
}
else
{
currentTestingCellIndex =
-1; // Either seperated by previous visit, boundary or non-manifold
}
// cells is smaller than the thresold and the nighboring cell has not been visited
}
}
regionIndex++;
}
}
// For each new region you need a new point
newPointNum = regionIndex - 1;
vtkm::Id numberOfCellsNeedUpdate = 0;
for (vtkm::IdComponent i = 0; i < numberOfIncidentCells; i++)
{
if (visitedCellsRegionIndex[i] > 0)
{
numberOfCellsNeedUpdate++;
}
}
cellNum = numberOfCellsNeedUpdate;
}
private:
vtkm::FloatDefault CosFeatureAngle; // Cos value of the feature angle
};
// This worklet split the sharp edges and populate the
// cellTopologyUpdateTuples as (cellGlobalId, oldPointId, newPointId).
class SplitSharpEdge : public vtkm::worklet::WorkletMapCellToPoint
{
public:
SplitSharpEdge(vtkm::FloatDefault cosfeatureAngle, vtkm::Id numberOfOldPoints)
: CosFeatureAngle(cosfeatureAngle)
, NumberOfOldPoints(numberOfOldPoints)
{
}
using ControlSignature = void(CellSetIn intputCells,
WholeCellSetIn<Point, Cell>, // Query points from cell
FieldInCell<Vec3> faceNormals,
FieldInPoint<IdType> newPointStartingIndex,
FieldInPoint<IdType> pointCellsStartingIndex,
WholeArrayOut<Id3Type> cellTopologyUpdateTuples);
using ExecutionSignature = void(CellIndices incidentCells,
InputIndex pointIndex,
_2 pFromCellSet,
_3 faceNormals,
_4 newPointStartingIndex,
_5 pointCellsStartingIndex,
_6 cellTopologyUpdateTuples);
using InputDomain = _1;
template <typename IncidentCellVecType,
typename PointFromCellSetType,
typename FaceNormalVecType,
typename CellTopologyUpdateTuples>
VTKM_EXEC void operator()(const IncidentCellVecType& incidentCells,
vtkm::Id pointIndex,
const PointFromCellSetType& pFromCellSet,
const FaceNormalVecType& faceNormals,
const vtkm::Id& newPointStartingIndex,
const vtkm::Id& pointCellsStartingIndex,
CellTopologyUpdateTuples& cellTopologyUpdateTuples) const
{
vtkm::IdComponent numberOfIncidentCells = incidentCells.GetNumberOfComponents();
if (numberOfIncidentCells <= 1)
{
return; // Not enought cells to compare
}
// Initialize a global cell mask to avoid confusion. globalCellIndex->status
// 0 means not visited yet 1 means visited.
vtkm::Bitset<vtkm::UInt64> visitedCells;
vtkm::Id visitedCellsRegionIndex[64] = { 0 };
// Reallocate memory for visitedCellsGroup if needed
int regionIndex = 0;
// Loop through each cell
for (vtkm::IdComponent incidentCellIndex = 0; incidentCellIndex < numberOfIncidentCells;
incidentCellIndex++)
{
vtkm::Id cellIndexG = incidentCells[incidentCellIndex]; // cell index in global order
// If not visited
if (!visitedCells.test(incidentCellIndex))
{
// Mark the cell and track the region
visitedCells.set(incidentCellIndex);
visitedCellsRegionIndex[incidentCellIndex] = regionIndex;
// Find two edges containing the current point in canonial index
vtkm::Id2 edge0G(-1, -1), edge1G(-1, -1);
internal::FindRelatedEdges(pointIndex, cellIndexG, pFromCellSet, edge0G, edge1G, *this);
// Grow the area along each edge
for (size_t i = 0; i < 2; i++)
{ // Reset these two values for each grow operation
vtkm::Id2 currentEdgeG = i == 0 ? edge0G : edge1G;
vtkm::IdComponent currentTestingCellIndex = incidentCellIndex;
while (currentTestingCellIndex >= 0)
{
// Find the neighbor of the current cell edge in local index
int neighboringCellIndexQuery = internal::FindNeighborCellInLocalIndex(
currentEdgeG, pFromCellSet, incidentCells, currentTestingCellIndex, *this);
// The edge should be manifold and the neighboring cell should
// have not been visited
if (neighboringCellIndexQuery != -1 && !visitedCells.test(neighboringCellIndexQuery))
{
vtkm::IdComponent neighborCellIndex =
static_cast<vtkm::IdComponent>(neighboringCellIndexQuery);
// Try to grow the area if the feature angle between current neighbor
auto thisNormal = faceNormals[currentTestingCellIndex];
//neighborNormal
auto neighborNormal = faceNormals[neighborCellIndex];
if (vtkm::dot(thisNormal, neighborNormal) > this->CosFeatureAngle)
{ // No need to split.
visitedCells.set(neighborCellIndex);
// Move to examine next cell
currentTestingCellIndex = neighborCellIndex;
vtkm::Id2 neighborCellEdge0G(-1, -1), neighborCellEdge1G(-1, -1);
internal::FindRelatedEdges(pointIndex,
incidentCells[currentTestingCellIndex],
pFromCellSet,
neighborCellEdge0G,
neighborCellEdge1G,
*this);
// Update currentEdgeG
if ((currentEdgeG == neighborCellEdge0G) ||
currentEdgeG == vtkm::Id2(neighborCellEdge0G[1], neighborCellEdge0G[0]))
{
currentEdgeG = neighborCellEdge1G;
}
else
{
currentEdgeG = neighborCellEdge0G;
}
}
else
{
currentTestingCellIndex = -1;
}
}
else
{
currentTestingCellIndex =
-1; // Either seperated by previous visit, boundary or non-manifold
}
// cells is smaller than the thresold and the nighboring cell has not been visited
}
}
regionIndex++;
}
}
// For each new region you need a new point
// Initialize the offset in the global cellTopologyUpdateTuples;
vtkm::Id cellTopologyUpdateTuplesIndex = pointCellsStartingIndex;
for (vtkm::Id i = 0; i < numberOfIncidentCells; i++)
{
if (visitedCellsRegionIndex[i])
{ // New region generated. Need to update the topology
vtkm::Id replacementPointId =
NumberOfOldPoints + newPointStartingIndex + visitedCellsRegionIndex[i] - 1;
vtkm::Id globalCellId = incidentCells[static_cast<vtkm::IdComponent>(i)];
// (cellGlobalIndex, oldPointId, replacementPointId)
vtkm::Vec<vtkm::Id, 3> tuple =
vtkm::make_Vec(globalCellId, pointIndex, replacementPointId);
cellTopologyUpdateTuples.Set(cellTopologyUpdateTuplesIndex, tuple);
cellTopologyUpdateTuplesIndex++;
}
}
}
private:
vtkm::FloatDefault CosFeatureAngle; // Cos value of the feature angle
vtkm::Id NumberOfOldPoints;
};
template <typename CellSetType,
typename FaceNormalsType,
typename CoordsComType,
typename CoordsInStorageType,
typename CoordsOutStorageType,
typename NewCellSetType,
typename DeviceAdapter>
void Run(
const CellSetType& oldCellset,
const vtkm::FloatDefault featureAngle,
const FaceNormalsType& faceNormals,
const vtkm::cont::ArrayHandle<vtkm::Vec<CoordsComType, 3>, CoordsInStorageType>& oldCoords,
vtkm::cont::ArrayHandle<vtkm::Vec<CoordsComType, 3>, CoordsOutStorageType>& newCoords,
NewCellSetType& newCellset,
DeviceAdapter)
{
using Algorithm = vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapter>;
vtkm::FloatDefault featureAngleR =
featureAngle / static_cast<vtkm::FloatDefault>(180.0) * vtkm::Pi<vtkm::FloatDefault>();
ClassifyPoint classifyPoint(vtkm::Cos(featureAngleR));
vtkm::worklet::DispatcherMapTopology<ClassifyPoint> cpDispatcher(classifyPoint);
cpDispatcher.SetDevice(DeviceAdapter());
// Array of newPointNums and cellNeedUpdateNums
vtkm::cont::ArrayHandle<vtkm::Id> newPointNums, cellNeedUpdateNums;
cpDispatcher.Invoke(oldCellset, oldCellset, faceNormals, newPointNums, cellNeedUpdateNums);
vtkm::Id totalNewPointsNum = Algorithm::Reduce(newPointNums, vtkm::Id(0), vtkm::Add());
// Allocate the size for the updateCellTopologyArray
vtkm::Id cellsNeedUpdateNum = Algorithm::Reduce(cellNeedUpdateNums, vtkm::Id(0), vtkm::Add());
vtkm::cont::ArrayHandle<vtkm::Id3> cellTopologyUpdateTuples;
cellTopologyUpdateTuples.Allocate(cellsNeedUpdateNum);
vtkm::cont::ArrayHandle<vtkm::Id> newpointStartingIndexs, pointCellsStartingIndexs;
Algorithm::ScanExclusive(newPointNums, newpointStartingIndexs);
Algorithm::ScanExclusive(cellNeedUpdateNums, pointCellsStartingIndexs);
SplitSharpEdge splitSharpEdge(vtkm::Cos(featureAngleR), oldCoords.GetNumberOfValues());
vtkm::worklet::DispatcherMapTopology<SplitSharpEdge> sseDispatcher(splitSharpEdge);
sseDispatcher.SetDevice(DeviceAdapter());
sseDispatcher.Invoke(oldCellset,
oldCellset,
faceNormals,
newpointStartingIndexs,
pointCellsStartingIndexs,
cellTopologyUpdateTuples);
auto ctutPortal = cellTopologyUpdateTuples.GetPortalConstControl();
// Create the new point coordinate system and update NewPointsIdArray to
// process point field
this->NewPointsIdArray.Allocate(oldCoords.GetNumberOfValues() + totalNewPointsNum);
auto newPointsIdArrayPortal = this->NewPointsIdArray.GetPortalControl();
for (vtkm::Id i = 0; i < newPointNums.GetNumberOfValues(); i++)
{
newPointsIdArrayPortal.Set(i, i);
}
newCoords.Allocate(oldCoords.GetNumberOfValues() + totalNewPointsNum);
Algorithm::CopySubRange(oldCoords, 0, oldCoords.GetNumberOfValues(), newCoords);
vtkm::Id newCoordsIndex = oldCoords.GetNumberOfValues();
auto oldCoordsPortal = oldCoords.GetPortalConstControl();
auto newCoordsPortal = newCoords.GetPortalControl();
auto newPointNumsPortal = newPointNums.GetPortalControl();
for (vtkm::Id i = 0; i < newPointNums.GetNumberOfValues(); i++)
{ // Find out for each new point, how many times it should be added
for (vtkm::Id j = 0; j < newPointNumsPortal.Get(i); j++)
{
newPointsIdArrayPortal.Set(newCoordsIndex, i);
newCoordsPortal.Set(newCoordsIndex++, oldCoordsPortal.Get(i));
}
}
// Create the new cellset
CellDeepCopy::Run(oldCellset, newCellset, DeviceAdapter());
// FIXME: Since the non const get array function is not in CellSetExplict.h,
// here I just get a non-const copy of the array handle.
auto connectivityArrayHandle = newCellset.GetConnectivityArray(vtkm::TopologyElementTagPoint(),
vtkm::TopologyElementTagCell());
auto connectivityArrayHandleP = connectivityArrayHandle.GetPortalControl();
auto offsetArrayHandle = newCellset.GetIndexOffsetArray(vtkm::TopologyElementTagPoint(),
vtkm::TopologyElementTagCell());
auto offsetArrayHandleP = offsetArrayHandle.GetPortalControl();
for (vtkm::Id i = 0; i < cellTopologyUpdateTuples.GetNumberOfValues(); i++)
{
vtkm::Id cellId(ctutPortal.Get(i)[0]), oldPointId(ctutPortal.Get(i)[1]),
newPointId(ctutPortal.Get(i)[2]);
vtkm::Id bound = (cellId + 1 == offsetArrayHandle.GetNumberOfValues())
? connectivityArrayHandle.GetNumberOfValues()
: offsetArrayHandleP.Get(cellId + 1);
vtkm::Id k = 0;
for (vtkm::Id j = offsetArrayHandleP.Get(cellId); j < bound; j++, k++)
{
if (connectivityArrayHandleP.Get(j) == oldPointId)
{
connectivityArrayHandleP.Set(j, newPointId);
}
}
}
}
template <typename ValueType, typename StorageTag, typename DeviceTag>
vtkm::cont::ArrayHandle<ValueType> ProcessPointField(
const vtkm::cont::ArrayHandle<ValueType, StorageTag> in,
DeviceTag) const
{
using Algo = vtkm::cont::DeviceAdapterAlgorithm<DeviceTag>;
// Use a temporary permutation array to simplify the mapping:
auto tmp = vtkm::cont::make_ArrayHandlePermutation(this->NewPointsIdArray, in);
// Copy into an array with default storage:
vtkm::cont::ArrayHandle<ValueType> result;
Algo::Copy(tmp, result);
return result;
}
private:
vtkm::cont::ArrayHandle<vtkm::Id> NewPointsIdArray;
};
}
} // vtkm::worklet
#endif // vtk_m_worklet_SplitSharpEdges_h

1
vtkm/worklet/testing/CMakeLists.txt
View File

@ -62,6 +62,7 @@ set(unit_tests
UnitTestScatterCounting.cxx
UnitTestScatterPermutation.cxx
UnitTestSplatKernels.cxx
UnitTestSplitSharpEdges.cxx
UnitTestStreamingSine.cxx
UnitTestStreamLineUniformGrid.cxx
UnitTestSurfaceNormals.cxx

246
vtkm/worklet/testing/UnitTestSplitSharpEdges.cxx Normal file
View File

@ -0,0 +1,246 @@
//============================================================================
// 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/worklet/SplitSharpEdges.h>
#include <vtkm/worklet/SurfaceNormals.h>
#include <vtkm/cont/testing/MakeTestDataSet.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/worklet/DispatcherMapTopology.h>
namespace
{
using NormalsArrayHandle = vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::FloatDefault, 3>>;
using DeviceAdapter = VTKM_DEFAULT_DEVICE_ADAPTER_TAG;
const vtkm::Vec<vtkm::FloatDefault, 3> expectedCoords[24] = {
{ 0.0, 0.0, 0.0 }, { 1.0, 0.0, 0.0 }, { 1.0, 0.0, 1.0 }, { 0.0, 0.0, 1.0 }, { 0.0, 1.0, 0.0 },
{ 1.0, 1.0, 0.0 }, { 1.0, 1.0, 1.0 }, { 0.0, 1.0, 1.0 }, { 0.0, 0.0, 0.0 }, { 0.0, 0.0, 0.0 },
{ 1.0, 0.0, 0.0 }, { 1.0, 0.0, 0.0 }, { 1.0, 0.0, 1.0 }, { 1.0, 0.0, 1.0 }, { 0.0, 0.0, 1.0 },
{ 0.0, 0.0, 1.0 }, { 0.0, 1.0, 0.0 }, { 0.0, 1.0, 0.0 }, { 1.0, 1.0, 0.0 }, { 1.0, 1.0, 0.0 },
{ 1.0, 1.0, 1.0 }, { 1.0, 1.0, 1.0 }, { 0.0, 1.0, 1.0 }, { 0.0, 1.0, 1.0 }
};
const std::vector<vtkm::FloatDefault> expectedPointvar{ 10.1f, 20.1f, 30.2f, 40.2f, 50.3f, 60.3f,
70.3f, 80.3f, 10.1f, 10.1f, 20.1f, 20.1f,
30.2f, 30.2f, 40.2f, 40.2f, 50.3f, 50.3f,
60.3f, 60.3f, 70.3f, 70.3f, 80.3f, 80.3f };
const std::vector<vtkm::Id> expectedConnectivityArray91{ 0, 1, 5, 4, 1, 2, 6, 5, 2, 3, 7, 6,
3, 0, 4, 7, 4, 5, 6, 7, 0, 3, 2, 1 };
vtkm::cont::DataSet Make3DExplicitSimpleCube()
{
vtkm::cont::DataSet dataSet;
vtkm::cont::DataSetBuilderExplicit dsb;
const int nVerts = 8;
const int nCells = 6;
using CoordType = vtkm::Vec<vtkm::FloatDefault, 3>;
std::vector<CoordType> coords = {
CoordType(0, 0, 0), // 0
CoordType(1, 0, 0), // 1
CoordType(1, 0, 1), // 2
CoordType(0, 0, 1), // 3
CoordType(0, 1, 0), // 4
CoordType(1, 1, 0), // 5
CoordType(1, 1, 1), // 6
CoordType(0, 1, 1) // 7
};
//Connectivity
std::vector<vtkm::UInt8> shapes;
std::vector<vtkm::IdComponent> numIndices;
for (size_t i = 0; i < 6; i++)
{
shapes.push_back(vtkm::CELL_SHAPE_QUAD);
numIndices.push_back(4);
}
std::vector<vtkm::Id> conn;
// Down face
conn.push_back(0);
conn.push_back(1);
conn.push_back(5);
conn.push_back(4);
// Right face
conn.push_back(1);
conn.push_back(2);
conn.push_back(6);
conn.push_back(5);
// Top face
conn.push_back(2);
conn.push_back(3);
conn.push_back(7);
conn.push_back(6);
// Left face
conn.push_back(3);
conn.push_back(0);
conn.push_back(4);
conn.push_back(7);
// Front face
conn.push_back(4);
conn.push_back(5);
conn.push_back(6);
conn.push_back(7);
// Back face
conn.push_back(0);
conn.push_back(3);
conn.push_back(2);
conn.push_back(1);
//Create the dataset.
dataSet = dsb.Create(coords, shapes, numIndices, conn, "coordinates", "cells");
vtkm::FloatDefault vars[nVerts] = { 10.1f, 20.1f, 30.2f, 40.2f, 50.3f, 60.3f, 70.3f, 80.3f };
vtkm::FloatDefault cellvar[nCells] = { 100.1f, 200.2f, 300.3f, 400.4f, 500.5f, 600.6f };
vtkm::cont::DataSetFieldAdd dsf;
dsf.AddPointField(dataSet, "pointvar", vars, nVerts);
dsf.AddCellField(dataSet, "cellvar", cellvar, nCells, "cells");
return dataSet;
}
void TestSplitSharpEdgesSplitEveryEdge(vtkm::cont::DataSet& simpleCube,
NormalsArrayHandle& faceNormals,
vtkm::worklet::SplitSharpEdges& splitSharpEdges)
{ // Split every edge
vtkm::FloatDefault featureAngle = 89.0;
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::FloatDefault, 3>> newCoords;
vtkm::cont::CellSetExplicit<> newCellset;
splitSharpEdges.Run(simpleCube.GetCellSet(),
featureAngle,
faceNormals,
simpleCube.GetCoordinateSystem().GetData(),
newCoords,
newCellset,
DeviceAdapter());
vtkm::cont::ArrayHandle<vtkm::FloatDefault> pointvar;
simpleCube.GetPointField("pointvar").GetData().CopyTo(pointvar);
vtkm::cont::ArrayHandle<vtkm::FloatDefault> newPointFields =
splitSharpEdges.ProcessPointField(pointvar, DeviceAdapter());
VTKM_TEST_ASSERT(newCoords.GetNumberOfValues() == 24,
"new coordinates"
" number is wrong");
auto newCoordsP = newCoords.GetPortalConstControl();
for (vtkm::Id i = 0; i < newCoords.GetNumberOfValues(); i++)
{
VTKM_TEST_ASSERT(test_equal(newCoordsP.Get(i)[0], expectedCoords[vtkm::IdComponent(i)][0]),
"result value does not match expected value");
VTKM_TEST_ASSERT(test_equal(newCoordsP.Get(i)[1], expectedCoords[vtkm::IdComponent(i)][1]),
"result value does not match expected value");
VTKM_TEST_ASSERT(test_equal(newCoordsP.Get(i)[2], expectedCoords[vtkm::IdComponent(i)][2]),
"result value does not match expected value");
}
auto newPointFieldsPortal = newPointFields.GetPortalConstControl();
for (int i = 0; i < newPointFields.GetNumberOfValues(); i++)
{
VTKM_TEST_ASSERT(
test_equal(newPointFieldsPortal.Get(i), expectedPointvar[static_cast<unsigned long>(i)]),
"point field array result does not match expected value");
}
}
void TestSplitSharpEdgesNoSplit(vtkm::cont::DataSet& simpleCube,
NormalsArrayHandle& faceNormals,
vtkm::worklet::SplitSharpEdges& splitSharpEdges)
{ // Do nothing
vtkm::FloatDefault featureAngle = 91.0;
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::FloatDefault, 3>> newCoords;
vtkm::cont::CellSetExplicit<> newCellset;
splitSharpEdges.Run(simpleCube.GetCellSet(),
featureAngle,
faceNormals,
simpleCube.GetCoordinateSystem().GetData(),
newCoords,
newCellset,
DeviceAdapter());
vtkm::cont::ArrayHandle<vtkm::FloatDefault> pointvar;
simpleCube.GetPointField("pointvar").GetData().CopyTo(pointvar);
vtkm::cont::ArrayHandle<vtkm::FloatDefault> newPointFields =
splitSharpEdges.ProcessPointField(pointvar, DeviceAdapter());
VTKM_TEST_ASSERT(newCoords.GetNumberOfValues() == 8,
"new coordinates"
" number is wrong");
auto newCoordsP = newCoords.GetPortalConstControl();
for (int i = 0; i < newCoords.GetNumberOfValues(); i++)
{
VTKM_TEST_ASSERT(test_equal(newCoordsP.Get(i)[0], expectedCoords[i][0]),
"result value does not match expected value");
VTKM_TEST_ASSERT(test_equal(newCoordsP.Get(i)[1], expectedCoords[i][1]),
"result value does not match expected value");
VTKM_TEST_ASSERT(test_equal(newCoordsP.Get(i)[2], expectedCoords[i][2]),
"result value does not match expected value");
}
const auto& connectivityArray = newCellset.GetConnectivityArray(vtkm::TopologyElementTagPoint(),
vtkm::TopologyElementTagCell());
auto connectivityArrayPortal = connectivityArray.GetPortalConstControl();
for (int i = 0; i < connectivityArray.GetNumberOfValues(); i++)
{
VTKM_TEST_ASSERT(test_equal(connectivityArrayPortal.Get(i),
expectedConnectivityArray91[static_cast<unsigned long>(i)]),
"connectivity array result does not match expected value");
}
auto newPointFieldsPortal = newPointFields.GetPortalConstControl();
for (int i = 0; i < newPointFields.GetNumberOfValues(); i++)
{
VTKM_TEST_ASSERT(
test_equal(newPointFieldsPortal.Get(i), expectedPointvar[static_cast<unsigned long>(i)]),
"point field array result does not match expected value");
}
}
void TestSplitSharpEdges()
{
vtkm::cont::DataSet simpleCube = Make3DExplicitSimpleCube();
NormalsArrayHandle faceNormals;
vtkm::worklet::FacetedSurfaceNormals faceted;
faceted.Run(simpleCube.GetCellSet(),
simpleCube.GetCoordinateSystem().GetData(),
faceNormals,
DeviceAdapter());
vtkm::worklet::SplitSharpEdges splitSharpEdges;
TestSplitSharpEdgesSplitEveryEdge(simpleCube, faceNormals, splitSharpEdges);
TestSplitSharpEdgesNoSplit(simpleCube, faceNormals, splitSharpEdges);
}
} // anonymous namespace
int UnitTestSplitSharpEdges(int, char* [])
{
return vtkm::cont::testing::Testing::Run(TestSplitSharpEdges);
}