From 06ac9f72198f50ca983a3fc5b0b25c2bafc9f55a Mon Sep 17 00:00:00 2001 From: Brent Lessley Date: Mon, 14 Jan 2019 15:27:19 -0800 Subject: [PATCH] Revised version of the original mesh quality merge request --- examples/CMakeLists.txt | 41 +- examples/mesh_quality/CMakeLists.txt | 37 + examples/mesh_quality/MeshQuality.cxx | 286 ++++++ .../internal/DeviceAdapterAlgorithmCuda.h | 965 +++++++++--------- vtkm/exec/CMakeLists.txt | 3 + vtkm/exec/cellmetrics/CMakeLists.txt | 28 + .../cellmetrics/CellDiagonalRatioMetric.h | 236 +++++ vtkm/exec/cellmetrics/CellEdgeRatioMetric.h | 300 ++++++ vtkm/filter/CMakeLists.txt | 2 + vtkm/filter/MeshQuality.h | 101 ++ vtkm/filter/MeshQuality.hxx | 250 +++++ vtkm/filter/testing/CMakeLists.txt | 1 + .../testing/UnitTestMeshQualityFilter.cxx | 376 +++++++ vtkm/worklet/CMakeLists.txt | 1 + vtkm/worklet/MeshQuality.h | 129 +++ 15 files changed, 2255 insertions(+), 501 deletions(-) create mode 100644 examples/mesh_quality/CMakeLists.txt create mode 100644 examples/mesh_quality/MeshQuality.cxx create mode 100644 vtkm/exec/cellmetrics/CMakeLists.txt create mode 100644 vtkm/exec/cellmetrics/CellDiagonalRatioMetric.h create mode 100644 vtkm/exec/cellmetrics/CellEdgeRatioMetric.h create mode 100644 vtkm/filter/MeshQuality.h create mode 100644 vtkm/filter/MeshQuality.hxx create mode 100644 vtkm/filter/testing/UnitTestMeshQualityFilter.cxx create mode 100644 vtkm/worklet/MeshQuality.h diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 501cb6ea8..4f222afdf 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -13,29 +13,30 @@ if(VTKm_ENABLE_EXAMPLES) set(CMAKE_PREFIX_PATH ${VTKm_BINARY_DIR}/${VTKm_INSTALL_CONFIG_DIR}) - add_subdirectory(clipping) - add_subdirectory(contour_tree) - add_subdirectory(contour_tree_augmented) - add_subdirectory(cosmotools) - add_subdirectory(demo) - add_subdirectory(game_of_life) - add_subdirectory(hello_world) - add_subdirectory(histogram) - add_subdirectory(isosurface) - add_subdirectory(lagrangian) - add_subdirectory(multi_backend) - add_subdirectory(oscillator) - add_subdirectory(particle_advection) - add_subdirectory(redistribute_points) - add_subdirectory(rendering) - add_subdirectory(streamline) - add_subdirectory(temporal_advection) - add_subdirectory(tetrahedra) + #add_subdirectory(clipping) + #add_subdirectory(contour_tree) + #add_subdirectory(contour_tree_augmented) + #add_subdirectory(cosmotools) + #add_subdirectory(demo) + #add_subdirectory(game_of_life) + #add_subdirectory(hello_world) + #add_subdirectory(histogram) + #add_subdirectory(isosurface) + #add_subdirectory(lagrangian) + add_subdirectory(mesh_quality) + #add_subdirectory(multi_backend) + #add_subdirectory(oscillator) + #add_subdirectory(particle_advection) + #add_subdirectory(redistribute_points) + #add_subdirectory(rendering) + #add_subdirectory(streamline) + #add_subdirectory(temporal_advection) + #add_subdirectory(tetrahedra) endif() if (VTKm_ENABLE_TESTING) # These need to be fast to build as they will # be built each time we run the test - vtkm_test_against_install(rendering) - vtkm_test_against_install(histogram) + #vtkm_test_against_install(rendering) + #vtkm_test_against_install(histogram) endif() diff --git a/examples/mesh_quality/CMakeLists.txt b/examples/mesh_quality/CMakeLists.txt new file mode 100644 index 000000000..7bd7a3f16 --- /dev/null +++ b/examples/mesh_quality/CMakeLists.txt @@ -0,0 +1,37 @@ +##============================================================================= +## +## 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 2015 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +## Copyright 2015 UT-Battelle, LLC. +## Copyright 2015 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. +## +##============================================================================= +cmake_minimum_required(VERSION 3.8...3.14 FATAL_ERROR) +project(MeshQuality CXX) + +#Find the VTK-m package +find_package(VTKm REQUIRED QUIET) + +add_executable(MeshQuality MeshQuality.cxx) +target_link_libraries(MeshQuality PRIVATE vtkm_filter) + +if(TARGET vtkm::tbb) + target_compile_definitions(MeshQuality PRIVATE BUILDING_TBB_VERSION) +endif() + +if(TARGET vtkm::cuda) + set_source_files_properties(MeshQuality.cxx PROPERTIES LANGUAGE "CUDA") +endif() diff --git a/examples/mesh_quality/MeshQuality.cxx b/examples/mesh_quality/MeshQuality.cxx new file mode 100644 index 000000000..f67178a76 --- /dev/null +++ b/examples/mesh_quality/MeshQuality.cxx @@ -0,0 +1,286 @@ +//============================================================================ +// 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 2018 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +// Copyright 2018 UT-Battelle, LLC. +// Copyright 2018 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 +#include +#include +#include +#include +#include +#include +#include + + +/** + * This example program takes an input VTK unstructured grid file and computes + * the mesh quality of its cells. The mesh quality of a cell type (e.g., tetrahedron, + * hexahedron, triangle) is defined by a metric, which is user-specified. + * The summary statistics of this metric, computed over all cells of its assigned type, + * are written as a new field in the output dataset. There is an output mesh quality + * field for each different cell type that has been assigned a non-empty quality metric + * by the user. Additionally, an ending field is written with the output metric value for + * each cell in the input dataset. The output dataset is named according to the + * user-provided input parameter to this file, or "output.vtk" by default if no parameter + * is provided. +*/ + +//Adapted from vtkm/cont/testing/MakeTestDataSet.h +//Modified the content of the Make3DExplicitDataSetZoo() function +inline vtkm::cont::DataSet Make3DExplicitDataSet() +{ + vtkm::cont::DataSet dataSet; + vtkm::cont::DataSetBuilderExplicit dsb; + + using CoordType = vtkm::Vec; + + std::vector coords = { + { 0.00, 0.00, 0.00 }, { 1.00, 0.00, 0.00 }, { 2.00, 0.00, 0.00 }, { 0.00, 0.00, 1.00 }, + { 1.00, 0.00, 1.00 }, { 2.00, 0.00, 1.00 }, { 0.00, 1.00, 0.00 }, { 1.00, 1.00, 0.00 }, + { 2.00, 1.00, 0.00 }, { 0.00, 1.00, 1.00 }, { 1.00, 1.00, 1.00 }, { 2.00, 1.00, 1.00 }, + { 0.00, 2.00, 0.00 }, { 1.00, 2.00, 0.00 }, { 2.00, 2.00, 0.00 }, { 0.00, 2.00, 1.00 }, + { 1.00, 2.00, 1.00 }, { 2.00, 2.00, 1.00 }, { 1.00, 3.00, 1.00 }, { 2.75, 0.00, 1.00 }, + { 3.00, 0.00, 0.75 }, { 3.00, 0.25, 1.00 }, { 3.00, 1.00, 1.00 }, { 3.00, 1.00, 0.00 }, + { 2.57, 2.00, 1.00 }, { 3.00, 1.75, 1.00 }, { 3.00, 1.75, 0.75 }, { 3.00, 0.00, 0.00 }, + { 2.57, 0.42, 0.57 }, { 2.59, 1.43, 0.71 } + }; + + std::vector shapes; + std::vector numindices; + std::vector conn; + + //Construct the shapes/cells of the dataset + //This is a zoo of points, lines, polygons, and polyhedra + shapes.push_back(vtkm::CELL_SHAPE_HEXAHEDRON); + numindices.push_back(8); + conn.push_back(0); + conn.push_back(3); + conn.push_back(4); + conn.push_back(1); + conn.push_back(6); + conn.push_back(9); + conn.push_back(10); + conn.push_back(7); + + shapes.push_back(vtkm::CELL_SHAPE_HEXAHEDRON); + numindices.push_back(8); + conn.push_back(1); + conn.push_back(4); + conn.push_back(5); + conn.push_back(2); + conn.push_back(7); + conn.push_back(10); + conn.push_back(11); + conn.push_back(8); + + shapes.push_back(vtkm::CELL_SHAPE_TETRA); + numindices.push_back(4); + conn.push_back(24); + conn.push_back(26); + conn.push_back(25); + conn.push_back(29); + + shapes.push_back(vtkm::CELL_SHAPE_TETRA); + numindices.push_back(4); + conn.push_back(8); + conn.push_back(17); + conn.push_back(11); + conn.push_back(29); + + shapes.push_back(vtkm::CELL_SHAPE_PYRAMID); + numindices.push_back(5); + conn.push_back(24); + conn.push_back(17); + conn.push_back(8); + conn.push_back(23); + conn.push_back(29); + + shapes.push_back(vtkm::CELL_SHAPE_PYRAMID); + numindices.push_back(5); + conn.push_back(25); + conn.push_back(22); + conn.push_back(11); + conn.push_back(17); + conn.push_back(29); + + shapes.push_back(vtkm::CELL_SHAPE_WEDGE); + numindices.push_back(6); + conn.push_back(8); + conn.push_back(14); + conn.push_back(17); + conn.push_back(7); + conn.push_back(13); + conn.push_back(16); + + shapes.push_back(vtkm::CELL_SHAPE_WEDGE); + numindices.push_back(6); + conn.push_back(11); + conn.push_back(8); + conn.push_back(17); + conn.push_back(10); + conn.push_back(7); + conn.push_back(16); + + shapes.push_back(vtkm::CELL_SHAPE_VERTEX); + numindices.push_back(1); + conn.push_back(0); + + shapes.push_back(vtkm::CELL_SHAPE_VERTEX); + numindices.push_back(1); + conn.push_back(29); + + shapes.push_back(vtkm::CELL_SHAPE_LINE); + numindices.push_back(2); + conn.push_back(0); + conn.push_back(1); + + shapes.push_back(vtkm::CELL_SHAPE_LINE); + numindices.push_back(2); + conn.push_back(11); + conn.push_back(15); + + shapes.push_back(vtkm::CELL_SHAPE_TRIANGLE); + numindices.push_back(3); + conn.push_back(2); + conn.push_back(4); + conn.push_back(15); + + shapes.push_back(vtkm::CELL_SHAPE_TRIANGLE); + numindices.push_back(3); + conn.push_back(5); + conn.push_back(6); + conn.push_back(7); + + shapes.push_back(vtkm::CELL_SHAPE_QUAD); + numindices.push_back(4); + conn.push_back(0); + conn.push_back(3); + conn.push_back(5); + conn.push_back(2); + + shapes.push_back(vtkm::CELL_SHAPE_QUAD); + numindices.push_back(4); + conn.push_back(5); + conn.push_back(4); + conn.push_back(10); + conn.push_back(11); + + shapes.push_back(vtkm::CELL_SHAPE_POLYGON); + numindices.push_back(3); + conn.push_back(4); + conn.push_back(7); + conn.push_back(1); + + shapes.push_back(vtkm::CELL_SHAPE_POLYGON); + numindices.push_back(4); + conn.push_back(1); + conn.push_back(6); + conn.push_back(7); + conn.push_back(2); + + dataSet = dsb.Create(coords, shapes, numindices, conn, "coordinates", "cells"); + + return dataSet; +} + +template +int TestMetrics(const char* outFileName, vtkm::cont::DataSet data, T filter) +{ + vtkm::cont::DataSet outputData; + try + { + vtkm::io::writer::VTKDataSetWriter writer("testZoo_withPolygons.vtk"); + writer.WriteDataSet(data); + + outputData = filter.Execute(data); + std::cout << "filter finished\n"; + } + catch (vtkm::cont::ErrorExecution&) + { + //TODO: need to add something else here... + std::cerr << "Error occured while executing the filter. Exiting" << std::endl; + return 1; + } + try + { + vtkm::io::writer::VTKDataSetWriter writer(outFileName); + writer.WriteDataSet(outputData); + std::cout << "finished writing data\n"; + } + catch (vtkm::io::ErrorIO&) + { + //TODO: need to add something else here... + std::cerr << "Error occured while writing the output data set. Exiting" << std::endl; + return 1; + } + return 0; +} + +int main(int argc, char* argv[]) +{ + // + // Check usage, set filename and read input + // + const char* outFileName = nullptr; + switch (argc) + { + case 1: + std::cerr << "Usage: " << std::endl + << "$ " << argv[0] << " " << std::endl; + return 1; + case 2: + outFileName = "output.vtk"; + break; + default: + outFileName = argv[argc - 1]; + break; + } + vtkm::cont::DataSet input; + vtkm::io::reader::VTKDataSetReader reader(argv[1]); + + + //Assign a cell metric to compute for each different + //shape type that may exist in the input dataset. If no metric + //is specified for a shape type, then it is assumed to be EMPTY + //and no metric is computed. + std::vector> shapeMetrics{ + vtkm::make_Pair(vtkm::CELL_SHAPE_LINE, vtkm::filter::CellMetric::VOLUME), + vtkm::make_Pair(vtkm::CELL_SHAPE_TRIANGLE, vtkm::filter::CellMetric::VOLUME), + vtkm::make_Pair(vtkm::CELL_SHAPE_POLYGON, vtkm::filter::CellMetric::VOLUME), + vtkm::make_Pair(vtkm::CELL_SHAPE_TETRA, vtkm::filter::CellMetric::VOLUME), + vtkm::make_Pair(vtkm::CELL_SHAPE_HEXAHEDRON, vtkm::filter::CellMetric::VOLUME), + vtkm::make_Pair(vtkm::CELL_SHAPE_WEDGE, vtkm::filter::CellMetric::VOLUME), + vtkm::make_Pair(vtkm::CELL_SHAPE_QUAD, vtkm::filter::CellMetric::VOLUME), + vtkm::make_Pair(vtkm::CELL_SHAPE_PYRAMID, vtkm::filter::CellMetric::VOLUME) + }; + + vtkm::filter::MeshQuality filter(shapeMetrics); + try + { + input = reader.ReadDataSet(); //FIELD not supported errors here, but doesnt affect data + //input = Make3DExplicitDataSet(); + } + catch (vtkm::io::ErrorIO&) + { + std::cerr << "Error occured while reading input. Exiting" << std::endl; + return 1; + } + TestMetrics(outFileName, input, filter); + return 0; +} diff --git a/vtkm/cont/cuda/internal/DeviceAdapterAlgorithmCuda.h b/vtkm/cont/cuda/internal/DeviceAdapterAlgorithmCuda.h index c0f5edc09..fea9926e9 100644 --- a/vtkm/cont/cuda/internal/DeviceAdapterAlgorithmCuda.h +++ b/vtkm/cont/cuda/internal/DeviceAdapterAlgorithmCuda.h @@ -1223,549 +1223,552 @@ public: } template - //resized array - vtkm::cont::ArrayHandle temp; - temp.Allocate(copyOutEnd); - CopySubRange(output, 0, outSize, temp); - output = temp; -} -} -CopySubRangePortal(input.PrepareForInput(DeviceAdapterTagCuda()), - inputStartIndex, - numberOfElementsToCopy, - output.PrepareForInPlace(DeviceAdapterTagCuda()), - outputIndex); -return true; -} - -template -VTKM_CONT static void LowerBounds(const vtkm::cont::ArrayHandle& input, - const vtkm::cont::ArrayHandle& values, - vtkm::cont::ArrayHandle& output) -{ - VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); - - vtkm::Id numberOfValues = values.GetNumberOfValues(); - LowerBoundsPortal(input.PrepareForInput(DeviceAdapterTagCuda()), - values.PrepareForInput(DeviceAdapterTagCuda()), - output.PrepareForOutput(numberOfValues, DeviceAdapterTagCuda())); -} - -template -VTKM_CONT static void LowerBounds(const vtkm::cont::ArrayHandle& input, - const vtkm::cont::ArrayHandle& values, - vtkm::cont::ArrayHandle& output, - BinaryCompare binary_compare) -{ - VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); - - vtkm::Id numberOfValues = values.GetNumberOfValues(); - LowerBoundsPortal(input.PrepareForInput(DeviceAdapterTagCuda()), - values.PrepareForInput(DeviceAdapterTagCuda()), - output.PrepareForOutput(numberOfValues, DeviceAdapterTagCuda()), - binary_compare); -} - -template -VTKM_CONT static void LowerBounds(const vtkm::cont::ArrayHandle& input, - vtkm::cont::ArrayHandle& values_output) -{ - VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); - - LowerBoundsPortal(input.PrepareForInput(DeviceAdapterTagCuda()), - values_output.PrepareForInPlace(DeviceAdapterTagCuda())); -} - -template -VTKM_CONT static U Reduce(const vtkm::cont::ArrayHandle& input, U initialValue) -{ - VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); - - const vtkm::Id numberOfValues = input.GetNumberOfValues(); - if (numberOfValues <= 0) + VTKM_CONT static void LowerBounds(const vtkm::cont::ArrayHandle& input, + const vtkm::cont::ArrayHandle& values, + vtkm::cont::ArrayHandle& output) { - return initialValue; - } - return ReducePortal(input.PrepareForInput(DeviceAdapterTagCuda()), initialValue); -} + VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); -template -VTKM_CONT static U Reduce(const vtkm::cont::ArrayHandle& input, - U initialValue, - BinaryFunctor binary_functor) -{ - VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); - - const vtkm::Id numberOfValues = input.GetNumberOfValues(); - if (numberOfValues <= 0) - { - return initialValue; - } - return ReducePortal(input.PrepareForInput(DeviceAdapterTagCuda()), initialValue, binary_functor); -} - -template -VTKM_CONT static void ReduceByKey(const vtkm::cont::ArrayHandle& keys, - const vtkm::cont::ArrayHandle& values, - vtkm::cont::ArrayHandle& keys_output, - vtkm::cont::ArrayHandle& values_output, - BinaryFunctor binary_functor) -{ - VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); - - //there is a concern that by default we will allocate too much - //space for the keys/values output. 1 option is to - const vtkm::Id numberOfValues = keys.GetNumberOfValues(); - if (numberOfValues <= 0) - { - return; - } - vtkm::Id reduced_size = - ReduceByKeyPortal(keys.PrepareForInput(DeviceAdapterTagCuda()), + vtkm::Id numberOfValues = values.GetNumberOfValues(); + LowerBoundsPortal(input.PrepareForInput(DeviceAdapterTagCuda()), values.PrepareForInput(DeviceAdapterTagCuda()), - keys_output.PrepareForOutput(numberOfValues, DeviceAdapterTagCuda()), - values_output.PrepareForOutput(numberOfValues, DeviceAdapterTagCuda()), - binary_functor); - - keys_output.Shrink(reduced_size); - values_output.Shrink(reduced_size); -} - -template -VTKM_CONT static T ScanExclusive(const vtkm::cont::ArrayHandle& input, - vtkm::cont::ArrayHandle& output) -{ - VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); - - const vtkm::Id numberOfValues = input.GetNumberOfValues(); - if (numberOfValues <= 0) - { - output.PrepareForOutput(0, DeviceAdapterTagCuda()); - return vtkm::TypeTraits::ZeroInitialization(); + output.PrepareForOutput(numberOfValues, DeviceAdapterTagCuda())); } - //We need call PrepareForInput on the input argument before invoking a - //function. The order of execution of parameters of a function is undefined - //so we need to make sure input is called before output, or else in-place - //use case breaks. - auto inputPortal = input.PrepareForInput(DeviceAdapterTagCuda()); - return ScanExclusivePortal(inputPortal, - output.PrepareForOutput(numberOfValues, DeviceAdapterTagCuda())); -} - -template -VTKM_CONT static T ScanExclusive(const vtkm::cont::ArrayHandle& input, - vtkm::cont::ArrayHandle& output, - BinaryFunctor binary_functor, - const T& initialValue) -{ - VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); - - const vtkm::Id numberOfValues = input.GetNumberOfValues(); - if (numberOfValues <= 0) + template + VTKM_CONT static void LowerBounds(const vtkm::cont::ArrayHandle& input, + const vtkm::cont::ArrayHandle& values, + vtkm::cont::ArrayHandle& output, + BinaryCompare binary_compare) { - output.PrepareForOutput(0, DeviceAdapterTagCuda()); - return vtkm::TypeTraits::ZeroInitialization(); + VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); + + vtkm::Id numberOfValues = values.GetNumberOfValues(); + LowerBoundsPortal(input.PrepareForInput(DeviceAdapterTagCuda()), + values.PrepareForInput(DeviceAdapterTagCuda()), + output.PrepareForOutput(numberOfValues, DeviceAdapterTagCuda()), + binary_compare); } - //We need call PrepareForInput on the input argument before invoking a - //function. The order of execution of parameters of a function is undefined - //so we need to make sure input is called before output, or else in-place - //use case breaks. - auto inputPortal = input.PrepareForInput(DeviceAdapterTagCuda()); - return ScanExclusivePortal(inputPortal, + template + VTKM_CONT static void LowerBounds(const vtkm::cont::ArrayHandle& input, + vtkm::cont::ArrayHandle& values_output) + { + VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); + + LowerBoundsPortal(input.PrepareForInput(DeviceAdapterTagCuda()), + values_output.PrepareForInPlace(DeviceAdapterTagCuda())); + } + + template + VTKM_CONT static U Reduce(const vtkm::cont::ArrayHandle& input, U initialValue) + { + VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); + + const vtkm::Id numberOfValues = input.GetNumberOfValues(); + if (numberOfValues <= 0) + { + return initialValue; + } + return ReducePortal(input.PrepareForInput(DeviceAdapterTagCuda()), initialValue); + } + + template + VTKM_CONT static U Reduce(const vtkm::cont::ArrayHandle& input, + U initialValue, + BinaryFunctor binary_functor) + { + VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); + + const vtkm::Id numberOfValues = input.GetNumberOfValues(); + if (numberOfValues <= 0) + { + return initialValue; + } + return ReducePortal( + input.PrepareForInput(DeviceAdapterTagCuda()), initialValue, binary_functor); + } + + template + VTKM_CONT static void ReduceByKey(const vtkm::cont::ArrayHandle& keys, + const vtkm::cont::ArrayHandle& values, + vtkm::cont::ArrayHandle& keys_output, + vtkm::cont::ArrayHandle& values_output, + BinaryFunctor binary_functor) + { + VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); + + //there is a concern that by default we will allocate too much + //space for the keys/values output. 1 option is to + const vtkm::Id numberOfValues = keys.GetNumberOfValues(); + if (numberOfValues <= 0) + { + return; + } + vtkm::Id reduced_size = + ReduceByKeyPortal(keys.PrepareForInput(DeviceAdapterTagCuda()), + values.PrepareForInput(DeviceAdapterTagCuda()), + keys_output.PrepareForOutput(numberOfValues, DeviceAdapterTagCuda()), + values_output.PrepareForOutput(numberOfValues, DeviceAdapterTagCuda()), + binary_functor); + + keys_output.Shrink(reduced_size); + values_output.Shrink(reduced_size); + } + + template + VTKM_CONT static T ScanExclusive(const vtkm::cont::ArrayHandle& input, + vtkm::cont::ArrayHandle& output) + { + VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); + + const vtkm::Id numberOfValues = input.GetNumberOfValues(); + if (numberOfValues <= 0) + { + output.PrepareForOutput(0, DeviceAdapterTagCuda()); + return vtkm::TypeTraits::ZeroInitialization(); + } + + //We need call PrepareForInput on the input argument before invoking a + //function. The order of execution of parameters of a function is undefined + //so we need to make sure input is called before output, or else in-place + //use case breaks. + auto inputPortal = input.PrepareForInput(DeviceAdapterTagCuda()); + return ScanExclusivePortal(inputPortal, + output.PrepareForOutput(numberOfValues, DeviceAdapterTagCuda())); + } + + template + VTKM_CONT static T ScanExclusive(const vtkm::cont::ArrayHandle& input, + vtkm::cont::ArrayHandle& output, + BinaryFunctor binary_functor, + const T& initialValue) + { + VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); + + const vtkm::Id numberOfValues = input.GetNumberOfValues(); + if (numberOfValues <= 0) + { + output.PrepareForOutput(0, DeviceAdapterTagCuda()); + return vtkm::TypeTraits::ZeroInitialization(); + } + + //We need call PrepareForInput on the input argument before invoking a + //function. The order of execution of parameters of a function is undefined + //so we need to make sure input is called before output, or else in-place + //use case breaks. + auto inputPortal = input.PrepareForInput(DeviceAdapterTagCuda()); + return ScanExclusivePortal(inputPortal, + output.PrepareForOutput(numberOfValues, DeviceAdapterTagCuda()), + binary_functor, + initialValue); + } + + template + VTKM_CONT static T ScanInclusive(const vtkm::cont::ArrayHandle& input, + vtkm::cont::ArrayHandle& output) + { + VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); + + const vtkm::Id numberOfValues = input.GetNumberOfValues(); + if (numberOfValues <= 0) + { + output.PrepareForOutput(0, DeviceAdapterTagCuda()); + return vtkm::TypeTraits::ZeroInitialization(); + } + + //We need call PrepareForInput on the input argument before invoking a + //function. The order of execution of parameters of a function is undefined + //so we need to make sure input is called before output, or else in-place + //use case breaks. + auto inputPortal = input.PrepareForInput(DeviceAdapterTagCuda()); + return ScanInclusivePortal(inputPortal, + output.PrepareForOutput(numberOfValues, DeviceAdapterTagCuda())); + } + + template + VTKM_CONT static T ScanInclusive(const vtkm::cont::ArrayHandle& input, + vtkm::cont::ArrayHandle& output, + BinaryFunctor binary_functor) + { + VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); + + const vtkm::Id numberOfValues = input.GetNumberOfValues(); + if (numberOfValues <= 0) + { + output.PrepareForOutput(0, DeviceAdapterTagCuda()); + return vtkm::TypeTraits::ZeroInitialization(); + } + + //We need call PrepareForInput on the input argument before invoking a + //function. The order of execution of parameters of a function is undefined + //so we need to make sure input is called before output, or else in-place + //use case breaks. + auto inputPortal = input.PrepareForInput(DeviceAdapterTagCuda()); + return ScanInclusivePortal( + inputPortal, output.PrepareForOutput(numberOfValues, DeviceAdapterTagCuda()), binary_functor); + } + + template + VTKM_CONT static void ScanInclusiveByKey(const vtkm::cont::ArrayHandle& keys, + const vtkm::cont::ArrayHandle& values, + vtkm::cont::ArrayHandle& output) + { + VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); + + const vtkm::Id numberOfValues = keys.GetNumberOfValues(); + if (numberOfValues <= 0) + { + output.PrepareForOutput(0, DeviceAdapterTagCuda()); + } + + //We need call PrepareForInput on the input argument before invoking a + //function. The order of execution of parameters of a function is undefined + //so we need to make sure input is called before output, or else in-place + //use case breaks. + auto keysPortal = keys.PrepareForInput(DeviceAdapterTagCuda()); + auto valuesPortal = values.PrepareForInput(DeviceAdapterTagCuda()); + ScanInclusiveByKeyPortal( + keysPortal, valuesPortal, output.PrepareForOutput(numberOfValues, DeviceAdapterTagCuda())); + } + + template + VTKM_CONT static void ScanInclusiveByKey(const vtkm::cont::ArrayHandle& keys, + const vtkm::cont::ArrayHandle& values, + vtkm::cont::ArrayHandle& output, + BinaryFunctor binary_functor) + { + VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); + + const vtkm::Id numberOfValues = keys.GetNumberOfValues(); + if (numberOfValues <= 0) + { + output.PrepareForOutput(0, DeviceAdapterTagCuda()); + } + + //We need call PrepareForInput on the input argument before invoking a + //function. The order of execution of parameters of a function is undefined + //so we need to make sure input is called before output, or else in-place + //use case breaks. + auto keysPortal = keys.PrepareForInput(DeviceAdapterTagCuda()); + auto valuesPortal = values.PrepareForInput(DeviceAdapterTagCuda()); + ScanInclusiveByKeyPortal(keysPortal, + valuesPortal, output.PrepareForOutput(numberOfValues, DeviceAdapterTagCuda()), - binary_functor, - initialValue); -} - -template -VTKM_CONT static T ScanInclusive(const vtkm::cont::ArrayHandle& input, - vtkm::cont::ArrayHandle& output) -{ - VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); - - const vtkm::Id numberOfValues = input.GetNumberOfValues(); - if (numberOfValues <= 0) - { - output.PrepareForOutput(0, DeviceAdapterTagCuda()); - return vtkm::TypeTraits::ZeroInitialization(); + ::thrust::equal_to(), + binary_functor); } - //We need call PrepareForInput on the input argument before invoking a - //function. The order of execution of parameters of a function is undefined - //so we need to make sure input is called before output, or else in-place - //use case breaks. - auto inputPortal = input.PrepareForInput(DeviceAdapterTagCuda()); - return ScanInclusivePortal(inputPortal, - output.PrepareForOutput(numberOfValues, DeviceAdapterTagCuda())); -} - -template -VTKM_CONT static T ScanInclusive(const vtkm::cont::ArrayHandle& input, - vtkm::cont::ArrayHandle& output, - BinaryFunctor binary_functor) -{ - VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); - - const vtkm::Id numberOfValues = input.GetNumberOfValues(); - if (numberOfValues <= 0) + template + VTKM_CONT static void ScanExclusiveByKey(const vtkm::cont::ArrayHandle& keys, + const vtkm::cont::ArrayHandle& values, + vtkm::cont::ArrayHandle& output) { - output.PrepareForOutput(0, DeviceAdapterTagCuda()); - return vtkm::TypeTraits::ZeroInitialization(); + VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); + + const vtkm::Id numberOfValues = keys.GetNumberOfValues(); + if (numberOfValues <= 0) + { + output.PrepareForOutput(0, DeviceAdapterTagCuda()); + return; + } + + //We need call PrepareForInput on the input argument before invoking a + //function. The order of execution of parameters of a function is undefined + //so we need to make sure input is called before output, or else in-place + //use case breaks. + auto keysPortal = keys.PrepareForInput(DeviceAdapterTagCuda()); + auto valuesPortal = values.PrepareForInput(DeviceAdapterTagCuda()); + ScanExclusiveByKeyPortal(keysPortal, + valuesPortal, + output.PrepareForOutput(numberOfValues, DeviceAdapterTagCuda()), + vtkm::TypeTraits::ZeroInitialization(), + ::thrust::equal_to(), + vtkm::Add()); } - //We need call PrepareForInput on the input argument before invoking a - //function. The order of execution of parameters of a function is undefined - //so we need to make sure input is called before output, or else in-place - //use case breaks. - auto inputPortal = input.PrepareForInput(DeviceAdapterTagCuda()); - return ScanInclusivePortal( - inputPortal, output.PrepareForOutput(numberOfValues, DeviceAdapterTagCuda()), binary_functor); -} - -template -VTKM_CONT static void ScanInclusiveByKey(const vtkm::cont::ArrayHandle& keys, - const vtkm::cont::ArrayHandle& values, - vtkm::cont::ArrayHandle& output) -{ - VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); - - const vtkm::Id numberOfValues = keys.GetNumberOfValues(); - if (numberOfValues <= 0) + template + VTKM_CONT static void ScanExclusiveByKey(const vtkm::cont::ArrayHandle& keys, + const vtkm::cont::ArrayHandle& values, + vtkm::cont::ArrayHandle& output, + const U& initialValue, + BinaryFunctor binary_functor) { - output.PrepareForOutput(0, DeviceAdapterTagCuda()); + VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); + + const vtkm::Id numberOfValues = keys.GetNumberOfValues(); + if (numberOfValues <= 0) + { + output.PrepareForOutput(0, DeviceAdapterTagCuda()); + return; + } + + //We need call PrepareForInput on the input argument before invoking a + //function. The order of execution of parameters of a function is undefined + //so we need to make sure input is called before output, or else in-place + //use case breaks. + auto keysPortal = keys.PrepareForInput(DeviceAdapterTagCuda()); + auto valuesPortal = values.PrepareForInput(DeviceAdapterTagCuda()); + ScanExclusiveByKeyPortal(keysPortal, + valuesPortal, + output.PrepareForOutput(numberOfValues, DeviceAdapterTagCuda()), + initialValue, + ::thrust::equal_to(), + binary_functor); } - //We need call PrepareForInput on the input argument before invoking a - //function. The order of execution of parameters of a function is undefined - //so we need to make sure input is called before output, or else in-place - //use case breaks. - auto keysPortal = keys.PrepareForInput(DeviceAdapterTagCuda()); - auto valuesPortal = values.PrepareForInput(DeviceAdapterTagCuda()); - ScanInclusiveByKeyPortal( - keysPortal, valuesPortal, output.PrepareForOutput(numberOfValues, DeviceAdapterTagCuda())); -} - -template -VTKM_CONT static void ScanInclusiveByKey(const vtkm::cont::ArrayHandle& keys, - const vtkm::cont::ArrayHandle& values, - vtkm::cont::ArrayHandle& output, - BinaryFunctor binary_functor) -{ - VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); - - const vtkm::Id numberOfValues = keys.GetNumberOfValues(); - if (numberOfValues <= 0) + // we use cuda pinned memory to reduce the amount of synchronization + // and mem copies between the host and device. + struct VTKM_CONT_EXPORT PinnedErrorArray { - output.PrepareForOutput(0, DeviceAdapterTagCuda()); - } + char* HostPtr = nullptr; + char* DevicePtr = nullptr; + vtkm::Id Size = 0; + }; - //We need call PrepareForInput on the input argument before invoking a - //function. The order of execution of parameters of a function is undefined - //so we need to make sure input is called before output, or else in-place - //use case breaks. - auto keysPortal = keys.PrepareForInput(DeviceAdapterTagCuda()); - auto valuesPortal = values.PrepareForInput(DeviceAdapterTagCuda()); - ScanInclusiveByKeyPortal(keysPortal, - valuesPortal, - output.PrepareForOutput(numberOfValues, DeviceAdapterTagCuda()), - ::thrust::equal_to(), - binary_functor); -} + VTKM_CONT_EXPORT + static const PinnedErrorArray& GetPinnedErrorArray(); -template -VTKM_CONT static void ScanExclusiveByKey(const vtkm::cont::ArrayHandle& keys, - const vtkm::cont::ArrayHandle& values, - vtkm::cont::ArrayHandle& output) -{ - VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); + VTKM_CONT_EXPORT + static void CheckForErrors(); // throws vtkm::cont::ErrorExecution - const vtkm::Id numberOfValues = keys.GetNumberOfValues(); - if (numberOfValues <= 0) - { - output.PrepareForOutput(0, DeviceAdapterTagCuda()); - return; - } + VTKM_CONT_EXPORT + static void SetupErrorBuffer(vtkm::exec::cuda::internal::TaskStrided& functor); - //We need call PrepareForInput on the input argument before invoking a - //function. The order of execution of parameters of a function is undefined - //so we need to make sure input is called before output, or else in-place - //use case breaks. - auto keysPortal = keys.PrepareForInput(DeviceAdapterTagCuda()); - auto valuesPortal = values.PrepareForInput(DeviceAdapterTagCuda()); - ScanExclusiveByKeyPortal(keysPortal, - valuesPortal, - output.PrepareForOutput(numberOfValues, DeviceAdapterTagCuda()), - vtkm::TypeTraits::ZeroInitialization(), - ::thrust::equal_to(), - vtkm::Add()); -} + VTKM_CONT_EXPORT + static void GetBlocksAndThreads(vtkm::UInt32& blocks, + vtkm::UInt32& threadsPerBlock, + vtkm::Id size); -template -VTKM_CONT static void ScanExclusiveByKey(const vtkm::cont::ArrayHandle& keys, - const vtkm::cont::ArrayHandle& values, - vtkm::cont::ArrayHandle& output, - const U& initialValue, - BinaryFunctor binary_functor) -{ - VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); + VTKM_CONT_EXPORT + static void GetBlocksAndThreads(vtkm::UInt32& blocks, dim3& threadsPerBlock, const dim3& size); - const vtkm::Id numberOfValues = keys.GetNumberOfValues(); - if (numberOfValues <= 0) - { - output.PrepareForOutput(0, DeviceAdapterTagCuda()); - return; - } + VTKM_CONT_EXPORT + static void LogKernelLaunch(const cudaFuncAttributes& func_attrs, + const std::type_info& worklet_info, + vtkm::UInt32 blocks, + vtkm::UInt32 threadsPerBlock, + vtkm::Id size); - //We need call PrepareForInput on the input argument before invoking a - //function. The order of execution of parameters of a function is undefined - //so we need to make sure input is called before output, or else in-place - //use case breaks. - auto keysPortal = keys.PrepareForInput(DeviceAdapterTagCuda()); - auto valuesPortal = values.PrepareForInput(DeviceAdapterTagCuda()); - ScanExclusiveByKeyPortal(keysPortal, - valuesPortal, - output.PrepareForOutput(numberOfValues, DeviceAdapterTagCuda()), - initialValue, - ::thrust::equal_to(), - binary_functor); -} - -// we use cuda pinned memory to reduce the amount of synchronization -// and mem copies between the host and device. -struct VTKM_CONT_EXPORT PinnedErrorArray -{ - char* HostPtr = nullptr; - char* DevicePtr = nullptr; - vtkm::Id Size = 0; -}; - -VTKM_CONT_EXPORT -static const PinnedErrorArray& GetPinnedErrorArray(); - -VTKM_CONT_EXPORT -static void CheckForErrors(); // throws vtkm::cont::ErrorExecution - -VTKM_CONT_EXPORT -static void SetupErrorBuffer(vtkm::exec::cuda::internal::TaskStrided& functor); - -VTKM_CONT_EXPORT -static void GetBlocksAndThreads(vtkm::UInt32& blocks, vtkm::UInt32& threadsPerBlock, vtkm::Id size); - -VTKM_CONT_EXPORT -static void GetBlocksAndThreads(vtkm::UInt32& blocks, dim3& threadsPerBlock, const dim3& size); - -VTKM_CONT_EXPORT -static void LogKernelLaunch(const cudaFuncAttributes& func_attrs, - const std::type_info& worklet_info, - vtkm::UInt32 blocks, - vtkm::UInt32 threadsPerBlock, - vtkm::Id size); - -VTKM_CONT_EXPORT -static void LogKernelLaunch(const cudaFuncAttributes& func_attrs, - const std::type_info& worklet_info, - vtkm::UInt32 blocks, - dim3 threadsPerBlock, - const dim3& size); + VTKM_CONT_EXPORT + static void LogKernelLaunch(const cudaFuncAttributes& func_attrs, + const std::type_info& worklet_info, + vtkm::UInt32 blocks, + dim3 threadsPerBlock, + const dim3& size); public: -template -static void ScheduleTask(vtkm::exec::cuda::internal::TaskStrided1D& functor, - vtkm::Id numInstances) -{ - VTKM_ASSERT(numInstances >= 0); - if (numInstances < 1) + template + static void ScheduleTask(vtkm::exec::cuda::internal::TaskStrided1D& functor, + vtkm::Id numInstances) { - // No instances means nothing to run. Just return. - return; - } + VTKM_ASSERT(numInstances >= 0); + if (numInstances < 1) + { + // No instances means nothing to run. Just return. + return; + } - CheckForErrors(); - SetupErrorBuffer(functor); + CheckForErrors(); + SetupErrorBuffer(functor); - vtkm::UInt32 blocks, threadsPerBlock; - GetBlocksAndThreads(blocks, threadsPerBlock, numInstances); + vtkm::UInt32 blocks, threadsPerBlock; + GetBlocksAndThreads(blocks, threadsPerBlock, numInstances); #ifdef VTKM_ENABLE_LOGGING - if (GetStderrLogLevel() >= vtkm::cont::LogLevel::KernelLaunches) - { - using FunctorType = vtkm::exec::cuda::internal::TaskStrided1D; - cudaFuncAttributes empty_kernel_attrs; - VTKM_CUDA_CALL( - cudaFuncGetAttributes(&empty_kernel_attrs, cuda::internal::TaskStrided1DLaunch)); - LogKernelLaunch(empty_kernel_attrs, typeid(WType), blocks, threadsPerBlock, numInstances); - } + if (GetStderrLogLevel() >= vtkm::cont::LogLevel::KernelLaunches) + { + using FunctorType = vtkm::exec::cuda::internal::TaskStrided1D; + cudaFuncAttributes empty_kernel_attrs; + VTKM_CUDA_CALL(cudaFuncGetAttributes(&empty_kernel_attrs, + cuda::internal::TaskStrided1DLaunch)); + LogKernelLaunch(empty_kernel_attrs, typeid(WType), blocks, threadsPerBlock, numInstances); + } #endif - cuda::internal::TaskStrided1DLaunch<<>>( - functor, numInstances); -} - -template -static void ScheduleTask(vtkm::exec::cuda::internal::TaskStrided3D& functor, - vtkm::Id3 rangeMax) -{ - VTKM_ASSERT((rangeMax[0] >= 0) && (rangeMax[1] >= 0) && (rangeMax[2] >= 0)); - if ((rangeMax[0] < 1) || (rangeMax[1] < 1) || (rangeMax[2] < 1)) - { - // No instances means nothing to run. Just return. - return; + cuda::internal::TaskStrided1DLaunch<<>>( + functor, numInstances); } - CheckForErrors(); - SetupErrorBuffer(functor); + template + static void ScheduleTask(vtkm::exec::cuda::internal::TaskStrided3D& functor, + vtkm::Id3 rangeMax) + { + VTKM_ASSERT((rangeMax[0] >= 0) && (rangeMax[1] >= 0) && (rangeMax[2] >= 0)); + if ((rangeMax[0] < 1) || (rangeMax[1] < 1) || (rangeMax[2] < 1)) + { + // No instances means nothing to run. Just return. + return; + } - const dim3 ranges(static_cast(rangeMax[0]), - static_cast(rangeMax[1]), - static_cast(rangeMax[2])); + CheckForErrors(); + SetupErrorBuffer(functor); - vtkm::UInt32 blocks; - dim3 threadsPerBlock; - GetBlocksAndThreads(blocks, threadsPerBlock, ranges); + const dim3 ranges(static_cast(rangeMax[0]), + static_cast(rangeMax[1]), + static_cast(rangeMax[2])); + + vtkm::UInt32 blocks; + dim3 threadsPerBlock; + GetBlocksAndThreads(blocks, threadsPerBlock, ranges); #ifdef VTKM_ENABLE_LOGGING - if (GetStderrLogLevel() >= vtkm::cont::LogLevel::KernelLaunches) - { - using FunctorType = vtkm::exec::cuda::internal::TaskStrided3D; - cudaFuncAttributes empty_kernel_attrs; - VTKM_CUDA_CALL( - cudaFuncGetAttributes(&empty_kernel_attrs, cuda::internal::TaskStrided3DLaunch)); - LogKernelLaunch(empty_kernel_attrs, typeid(WType), blocks, threadsPerBlock, ranges); - } + if (GetStderrLogLevel() >= vtkm::cont::LogLevel::KernelLaunches) + { + using FunctorType = vtkm::exec::cuda::internal::TaskStrided3D; + cudaFuncAttributes empty_kernel_attrs; + VTKM_CUDA_CALL(cudaFuncGetAttributes(&empty_kernel_attrs, + cuda::internal::TaskStrided3DLaunch)); + LogKernelLaunch(empty_kernel_attrs, typeid(WType), blocks, threadsPerBlock, ranges); + } #endif - cuda::internal::TaskStrided3DLaunch<<>>(functor, - ranges); -} + cuda::internal::TaskStrided3DLaunch<<>>( + functor, ranges); + } -template -VTKM_CONT static void Schedule(Functor functor, vtkm::Id numInstances) -{ - VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); + template + VTKM_CONT static void Schedule(Functor functor, vtkm::Id numInstances) + { + VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); - vtkm::exec::cuda::internal::TaskStrided1D kernel(functor); + vtkm::exec::cuda::internal::TaskStrided1D kernel(functor); - ScheduleTask(kernel, numInstances); -} + ScheduleTask(kernel, numInstances); + } -template -VTKM_CONT static void Schedule(Functor functor, const vtkm::Id3& rangeMax) -{ - VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); + template + VTKM_CONT static void Schedule(Functor functor, const vtkm::Id3& rangeMax) + { + VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); - vtkm::exec::cuda::internal::TaskStrided3D kernel(functor); - ScheduleTask(kernel, rangeMax); -} + vtkm::exec::cuda::internal::TaskStrided3D kernel(functor); + ScheduleTask(kernel, rangeMax); + } -template -VTKM_CONT static void Sort(vtkm::cont::ArrayHandle& values) -{ - VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); + template + VTKM_CONT static void Sort(vtkm::cont::ArrayHandle& values) + { + VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); - SortPortal(values.PrepareForInPlace(DeviceAdapterTagCuda())); -} + SortPortal(values.PrepareForInPlace(DeviceAdapterTagCuda())); + } -template -VTKM_CONT static void Sort(vtkm::cont::ArrayHandle& values, - BinaryCompare binary_compare) -{ - VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); - - SortPortal(values.PrepareForInPlace(DeviceAdapterTagCuda()), binary_compare); -} - -template -VTKM_CONT static void SortByKey(vtkm::cont::ArrayHandle& keys, - vtkm::cont::ArrayHandle& values) -{ - VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); - - SortByKeyPortal(keys.PrepareForInPlace(DeviceAdapterTagCuda()), - values.PrepareForInPlace(DeviceAdapterTagCuda())); -} - -template -VTKM_CONT static void SortByKey(vtkm::cont::ArrayHandle& keys, - vtkm::cont::ArrayHandle& values, - BinaryCompare binary_compare) -{ - VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); - - SortByKeyPortal(keys.PrepareForInPlace(DeviceAdapterTagCuda()), - values.PrepareForInPlace(DeviceAdapterTagCuda()), - binary_compare); -} - -template -VTKM_CONT static void Unique(vtkm::cont::ArrayHandle& values) -{ - VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); - - vtkm::Id newSize = UniquePortal(values.PrepareForInPlace(DeviceAdapterTagCuda())); - - values.Shrink(newSize); -} - -template -VTKM_CONT static void Unique(vtkm::cont::ArrayHandle& values, + template + VTKM_CONT static void Sort(vtkm::cont::ArrayHandle& values, BinaryCompare binary_compare) -{ - VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); + { + VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); - vtkm::Id newSize = UniquePortal(values.PrepareForInPlace(DeviceAdapterTagCuda()), binary_compare); + SortPortal(values.PrepareForInPlace(DeviceAdapterTagCuda()), binary_compare); + } - values.Shrink(newSize); -} + template + VTKM_CONT static void SortByKey(vtkm::cont::ArrayHandle& keys, + vtkm::cont::ArrayHandle& values) + { + VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); -template -VTKM_CONT static void UpperBounds(const vtkm::cont::ArrayHandle& input, - const vtkm::cont::ArrayHandle& values, - vtkm::cont::ArrayHandle& output) -{ - VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); + SortByKeyPortal(keys.PrepareForInPlace(DeviceAdapterTagCuda()), + values.PrepareForInPlace(DeviceAdapterTagCuda())); + } - vtkm::Id numberOfValues = values.GetNumberOfValues(); - UpperBoundsPortal(input.PrepareForInput(DeviceAdapterTagCuda()), - values.PrepareForInput(DeviceAdapterTagCuda()), - output.PrepareForOutput(numberOfValues, DeviceAdapterTagCuda())); -} - -template -VTKM_CONT static void UpperBounds(const vtkm::cont::ArrayHandle& input, - const vtkm::cont::ArrayHandle& values, - vtkm::cont::ArrayHandle& output, + template + VTKM_CONT static void SortByKey(vtkm::cont::ArrayHandle& keys, + vtkm::cont::ArrayHandle& values, BinaryCompare binary_compare) -{ - VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); + { + VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); - vtkm::Id numberOfValues = values.GetNumberOfValues(); - UpperBoundsPortal(input.PrepareForInput(DeviceAdapterTagCuda()), - values.PrepareForInput(DeviceAdapterTagCuda()), - output.PrepareForOutput(numberOfValues, DeviceAdapterTagCuda()), + SortByKeyPortal(keys.PrepareForInPlace(DeviceAdapterTagCuda()), + values.PrepareForInPlace(DeviceAdapterTagCuda()), binary_compare); -} + } -template -VTKM_CONT static void UpperBounds(const vtkm::cont::ArrayHandle& input, - vtkm::cont::ArrayHandle& values_output) -{ - VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); + template + VTKM_CONT static void Unique(vtkm::cont::ArrayHandle& values) + { + VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); - UpperBoundsPortal(input.PrepareForInput(DeviceAdapterTagCuda()), - values_output.PrepareForInPlace(DeviceAdapterTagCuda())); -} + vtkm::Id newSize = UniquePortal(values.PrepareForInPlace(DeviceAdapterTagCuda())); -VTKM_CONT static void Synchronize() -{ - VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); + values.Shrink(newSize); + } - VTKM_CUDA_CALL(cudaStreamSynchronize(cudaStreamPerThread)); - CheckForErrors(); -} -} -; + template + VTKM_CONT static void Unique(vtkm::cont::ArrayHandle& values, + BinaryCompare binary_compare) + { + VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); + + vtkm::Id newSize = + UniquePortal(values.PrepareForInPlace(DeviceAdapterTagCuda()), binary_compare); + + values.Shrink(newSize); + } + + template + VTKM_CONT static void UpperBounds(const vtkm::cont::ArrayHandle& input, + const vtkm::cont::ArrayHandle& values, + vtkm::cont::ArrayHandle& output) + { + VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); + + vtkm::Id numberOfValues = values.GetNumberOfValues(); + UpperBoundsPortal(input.PrepareForInput(DeviceAdapterTagCuda()), + values.PrepareForInput(DeviceAdapterTagCuda()), + output.PrepareForOutput(numberOfValues, DeviceAdapterTagCuda())); + } + + template + VTKM_CONT static void UpperBounds(const vtkm::cont::ArrayHandle& input, + const vtkm::cont::ArrayHandle& values, + vtkm::cont::ArrayHandle& output, + BinaryCompare binary_compare) + { + VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); + + vtkm::Id numberOfValues = values.GetNumberOfValues(); + UpperBoundsPortal(input.PrepareForInput(DeviceAdapterTagCuda()), + values.PrepareForInput(DeviceAdapterTagCuda()), + output.PrepareForOutput(numberOfValues, DeviceAdapterTagCuda()), + binary_compare); + } + + template + VTKM_CONT static void UpperBounds(const vtkm::cont::ArrayHandle& input, + vtkm::cont::ArrayHandle& values_output) + { + VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); + + UpperBoundsPortal(input.PrepareForInput(DeviceAdapterTagCuda()), + values_output.PrepareForInPlace(DeviceAdapterTagCuda())); + } + + VTKM_CONT static void Synchronize() + { + VTKM_LOG_SCOPE_FUNCTION(vtkm::cont::LogLevel::Perf); + + VTKM_CUDA_CALL(cudaStreamSynchronize(cudaStreamPerThread)); + CheckForErrors(); + } +}; template <> class DeviceTaskTypes diff --git a/vtkm/exec/CMakeLists.txt b/vtkm/exec/CMakeLists.txt index 210f9ed0a..9f310994d 100644 --- a/vtkm/exec/CMakeLists.txt +++ b/vtkm/exec/CMakeLists.txt @@ -56,3 +56,6 @@ add_subdirectory(cuda) #----------------------------------------------------------------------------- add_subdirectory(testing) + +#----------------------------------------------------------------------------- +add_subdirectory(cellmetrics) diff --git a/vtkm/exec/cellmetrics/CMakeLists.txt b/vtkm/exec/cellmetrics/CMakeLists.txt new file mode 100644 index 000000000..c878006c2 --- /dev/null +++ b/vtkm/exec/cellmetrics/CMakeLists.txt @@ -0,0 +1,28 @@ +##============================================================================ +## 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. +##============================================================================ + +set(headers + CellDiagonalRatioMetric.h + CellEdgeRatioMetric.h + ) + + +#----------------------------------------------------------------------------- +vtkm_declare_headers(${headers}) diff --git a/vtkm/exec/cellmetrics/CellDiagonalRatioMetric.h b/vtkm/exec/cellmetrics/CellDiagonalRatioMetric.h new file mode 100644 index 000000000..00925fb64 --- /dev/null +++ b/vtkm/exec/cellmetrics/CellDiagonalRatioMetric.h @@ -0,0 +1,236 @@ +//============================================================================ +// 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 2018 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +// Copyright 2018 UT-Battelle, LLC. +// Copyright 2018 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_cellmetrics_CellDiagonalRatioMetric_h +#define vtk_m_exec_cellmetrics_CellDiagonalRatioMetric_h + +/* + * Mesh quality metric functions that compute the diagonal ratio of mesh cells. + * The diagonal ratio of a cell is defined as the length (magnitude) of the longest + * cell diagonal length divided by the length of the shortest cell diagonal length. + * + * These metric computations are adapted from the VTK implementation of the Verdict library, + * which provides a set of mesh/cell metrics for evaluating the geometric qualities of regions + * of mesh spaces. + * + * The edge ratio computations for a pyramid cell types is not defined in the + * VTK implementation, but is provided here. + * + * See: The Verdict Library Reference Manual (for per-cell-type metric formulae) + * See: vtk/ThirdParty/verdict/vtkverdict (for VTK code implementation of this metric) + */ + +#include "vtkm/CellShape.h" +#include "vtkm/CellTraits.h" +#include "vtkm/VecTraits.h" +#include "vtkm/VectorAnalysis.h" +#include "vtkm/exec/FunctorBase.h" + +#define UNUSED(expr) (void)(expr); + +namespace vtkm +{ +namespace exec +{ +namespace cellmetrics +{ + +using FloatType = vtkm::FloatDefault; + + +template +VTKM_EXEC inline OutType ComputeDiagonalRatio(const VecType& diagonals) +{ + const vtkm::Id numDiagonals = diagonals.GetNumberOfComponents(); + + //Compare diagonal lengths to determine the longest and shortest + //TODO: Could we use lambda expression here? + + FloatType d0Len = (FloatType)vtkm::MagnitudeSquared(diagonals[0]); + FloatType currLen, minLen = d0Len, maxLen = d0Len; + for (int i = 1; i < numDiagonals; i++) + { + currLen = (FloatType)vtkm::MagnitudeSquared(diagonals[i]); + if (currLen < minLen) + minLen = currLen; + if (currLen > maxLen) + maxLen = currLen; + } + + if (minLen < vtkm::NegativeInfinity()) + return vtkm::Infinity(); + + //Take square root because we only did magnitude squared before + OutType diagonalRatio = (OutType)vtkm::Sqrt(maxLen / minLen); + if (diagonalRatio > 0) + return vtkm::Min(diagonalRatio, vtkm::Infinity()); //normal case + + return vtkm::Max(diagonalRatio, OutType(-1) * vtkm::Infinity()); +} + + +// ========================= Unsupported cells ================================== + +// By default, cells have zero shape unless the shape type template is specialized below. +template +VTKM_EXEC OutType CellDiagonalRatioMetric(const vtkm::IdComponent& numPts, + const PointCoordVecType& pts, + CellShapeType shape, + const vtkm::exec::FunctorBase&) +{ + UNUSED(numPts); + UNUSED(pts); + UNUSED(shape); + return OutType(0.0); +} + +/* +//TODO: Should polygons be supported? Maybe call Quad or Triangle function... +template +VTKM_EXEC OutType CellDiagonalRatioMetric(const vtkm::IdComponent& numPts, + const PointCoordVecType& pts, + vtkm::CellShapeTagPolygon, + const vtkm::exec::FunctorBase& worklet) +{ + switch (numPts) + { + case 4: + return CellDiagonalRatioMetric(numPts, pts, vtkm::CellShapeTagQuad(), worklet); + default: + break; + } + return OutType(-1.0); +} +*/ + +/* +template +VTKM_EXEC OutType CellDiagonalRatioMetric(const vtkm::IdComponent&, + const PointCoordVecType&, + vtkm::CellShapeTagLine, + const vtkm::exec::FunctorBase& worklet) +{ + UNUSED(worklet); + return OutType(-1.0); +} + +template +VTKM_EXEC OutType CellDiagonalRatioMetric(const vtkm::IdComponent&, + const PointCoordVecType&, + vtkm::CellShapeTagTriangle, + const vtkm::exec::FunctorBase& worklet) +{ + UNUSED(worklet); + return OutType(-1.0); +} + +template +VTKM_EXEC OutType CellDiagonalRatioMetric(const vtkm::IdComponent&, + const PointCoordVecType&, + vtkm::CellShapeTagTetra, + const vtkm::exec::FunctorBase& worklet) +{ + UNUSED(worklet); + return OutType(-1.0); +} + +template +VTKM_EXEC OutType CellDiagonalRatioMetric(const vtkm::IdComponent&, + const PointCoordVecType&, + vtkm::CellShapeTagWedge, + const vtkm::exec::FunctorBase& worklet) +{ + UNUSED(worklet); + return OutType(-1.0); +} + +template +VTKM_EXEC OutType CellDiagonalRatioMetric(const vtkm::IdComponent&, + const PointCoordVecType&, + vtkm::CellShapeTagPyramid, + const vtkm::exec::FunctorBase& worklet) +{ + UNUSED(worklet); + return OutType(-1.0); +} +*/ + +// ========================= 2D cells ================================== +// Compute the diagonal ratio of a quadrilateral. +// Formula: Maximum diagonal length divided by minimum diagonal length +// Equals 1 for a unit square +// Full range: [1,FLOAT_MAX] +template +VTKM_EXEC OutType CellDiagonalRatioMetric(const vtkm::IdComponent& numPts, + const PointCoordVecType& pts, + vtkm::CellShapeTagQuad, + const vtkm::exec::FunctorBase& worklet) +{ + if (numPts != 4) + { + worklet.RaiseError("Diagonal ratio metric(quad) requires 4 points."); + return OutType(0.0); + } + + vtkm::IdComponent numDiagonals = 2; //pts.GetNumberOfComponents(); + + //The 2 diagonals of a quadrilateral + using Diagonal = typename PointCoordVecType::ComponentType; + const Diagonal QuadDiagonals[2] = { pts[2] - pts[0], pts[3] - pts[1] }; + + return vtkm::exec::cellmetrics::ComputeDiagonalRatio( + vtkm::make_VecC(QuadDiagonals, numDiagonals)); +} + +// ============================= 3D Volume cells ================================== +// Compute the diagonal ratio of a hexahedron. +// Formula: Maximum diagonal length divided by minimum diagonal length +// Equals 1 for a unit cube +// Acceptable Range: [0.65, 1] +// Normal Range: [0, 1] +// Full range: [1,FLOAT_MAX] +template +VTKM_EXEC OutType CellDiagonalRatioMetric(const vtkm::IdComponent& numPts, + const PointCoordVecType& pts, + vtkm::CellShapeTagHexahedron, + const vtkm::exec::FunctorBase& worklet) +{ + if (numPts != 8) + { + worklet.RaiseError("Diagonal ratio metric(hexahedron) requires 8 points."); + return OutType(0.0); + } + + vtkm::IdComponent numDiagonals = 4; //pts.GetNumberOfComponents(); + + //The 4 diagonals of a hexahedron + using Diagonal = typename PointCoordVecType::ComponentType; + const Diagonal HexDiagonals[4] = { + pts[6] - pts[0], pts[7] - pts[1], pts[4] - pts[2], pts[5] - pts[3] + }; + + return vtkm::exec::cellmetrics::ComputeDiagonalRatio( + vtkm::make_VecC(HexDiagonals, numDiagonals)); +} + +} // namespace cellmetrics +} // namespace exec +} // namespace vtkm + +#endif // vtk_m_exec_cellmetrics_CellEdgeRatioMetric_h diff --git a/vtkm/exec/cellmetrics/CellEdgeRatioMetric.h b/vtkm/exec/cellmetrics/CellEdgeRatioMetric.h new file mode 100644 index 000000000..20d7d1300 --- /dev/null +++ b/vtkm/exec/cellmetrics/CellEdgeRatioMetric.h @@ -0,0 +1,300 @@ +//============================================================================ +// 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_exec_cellmetrics_CellEdgeRatioMetric_h +#define vtk_m_exec_cellmetrics_CellEdgeRatioMetric_h + +/* + * Mesh quality metric functions that compute the edge ratio of mesh cells. + * The edge ratio of a cell is defined as the length (magnitude) of the longest + * cell edge divided by the length of the shortest cell edge. + * + * These metric computations are adapted from the VTK implementation of the Verdict library, + * which provides a set of mesh/cell metrics for evaluating the geometric qualities of regions + * of mesh spaces. + * + * The edge ratio computations for a pyramid cell types is not defined in the + * VTK implementation, but is provided here. + * + * See: The Verdict Library Reference Manual (for per-cell-type metric formulae) + * See: vtk/ThirdParty/verdict/vtkverdict (for VTK code implementation of this metric) + */ + +#include "vtkm/CellShape.h" +#include "vtkm/CellTraits.h" +#include "vtkm/VecTraits.h" +#include "vtkm/VectorAnalysis.h" +#include "vtkm/exec/FunctorBase.h" + +#define UNUSED(expr) (void)(expr); + +namespace vtkm +{ +namespace exec +{ +namespace cellmetrics +{ + +using FloatType = vtkm::FloatDefault; + + +template +VTKM_EXEC inline OutType ComputeEdgeRatio(const VecType& edges) +{ + const vtkm::Id numEdges = edges.GetNumberOfComponents(); + + //Compare edge lengths to determine the longest and shortest + //TODO: Could we use lambda expression here? + + FloatType e0Len = (FloatType)vtkm::MagnitudeSquared(edges[0]); + FloatType currLen, minLen = e0Len, maxLen = e0Len; + for (int i = 1; i < numEdges; i++) + { + currLen = (FloatType)vtkm::MagnitudeSquared(edges[i]); + if (currLen < minLen) + minLen = currLen; + if (currLen > maxLen) + maxLen = currLen; + } + + if (minLen < vtkm::NegativeInfinity()) + return vtkm::Infinity(); + + //Take square root because we only did magnitude squared before + OutType edgeRatio = (OutType)vtkm::Sqrt(maxLen / minLen); + if (edgeRatio > 0) + return vtkm::Min(edgeRatio, vtkm::Infinity()); //normal case + + return vtkm::Max(edgeRatio, OutType(-1) * vtkm::Infinity()); +} + + +// ========================= Unsupported cells ================================== + +// By default, cells have zero shape unless the shape type template is specialized below. +template +VTKM_EXEC OutType CellEdgeRatioMetric(const vtkm::IdComponent& numPts, + const PointCoordVecType& pts, + CellShapeType shape, + const vtkm::exec::FunctorBase&) +{ + UNUSED(numPts); + UNUSED(pts); + UNUSED(shape); + return OutType(0.0); +} + +// ========================= 2D cells ================================== + + +// Compute the edge ratio of a line. +// Formula: Maximum edge length divided by minimum edge length +// Trivially equals 1, since only a single edge +template +VTKM_EXEC OutType CellEdgeRatioMetric(const vtkm::IdComponent& numPts, + const PointCoordVecType& pts, + vtkm::CellShapeTagLine, + const vtkm::exec::FunctorBase& worklet) +{ + UNUSED(pts); + if (numPts < 2) + { + worklet.RaiseError("Degenerate line has no edge ratio."); + return OutType(0.0); + } + return OutType(1.0); +} + +// Compute the edge ratio of a triangle. +// Formula: Maximum edge length divided by minimum edge length +// Equals 1 for an equilateral unit triangle +// Acceptable range: [1,1.3] +// Full range: [1,FLOAT_MAX] +template +VTKM_EXEC OutType CellEdgeRatioMetric(const vtkm::IdComponent& numPts, + const PointCoordVecType& pts, + vtkm::CellShapeTagTriangle, + const vtkm::exec::FunctorBase& worklet) +{ + if (numPts != 3) + { + std::cout << "numPts triangle: " << numPts << "\n"; + worklet.RaiseError("Edge ratio metric(triangle) requires 3 points."); + return OutType(0.0); + } + + const vtkm::IdComponent numEdges = 3; //pts.GetNumberOfComponents(); + + //The 3 edges of a triangle + using Edge = typename PointCoordVecType::ComponentType; + const Edge TriEdges[3] = { pts[1] - pts[0], pts[2] - pts[1], pts[0] - pts[2] }; + return vtkm::exec::cellmetrics::ComputeEdgeRatio(vtkm::make_VecC(TriEdges, numEdges)); +} + + +// Compute the edge ratio of a quadrilateral. +// Formula: Maximum edge length divided by minimum edge length +// Equals 1 for a unit square +// Acceptable range: [1,1.3] +// Full range: [1,FLOAT_MAX] +template +VTKM_EXEC OutType CellEdgeRatioMetric(const vtkm::IdComponent& numPts, + const PointCoordVecType& pts, + vtkm::CellShapeTagQuad, + const vtkm::exec::FunctorBase& worklet) +{ + if (numPts != 4) + { + worklet.RaiseError("Edge ratio metric(quad) requires 4 points."); + return OutType(0.0); + } + + vtkm::IdComponent numEdges = 4; //pts.GetNumberOfComponents(); + + //The 4 edges of a quadrilateral + using Edge = typename PointCoordVecType::ComponentType; + const Edge QuadEdges[4] = { pts[1] - pts[0], pts[2] - pts[1], pts[3] - pts[2], pts[0] - pts[3] }; + + return vtkm::exec::cellmetrics::ComputeEdgeRatio(vtkm::make_VecC(QuadEdges, numEdges)); +} + + + +// ============================= 3D Volume cells ==================================i + +// Compute the edge ratio of a tetrahedron. +// Formula: Maximum edge length divided by minimum edge length +// Equals 1 for a unit equilateral tetrahedron +// Acceptable range: [1,3] +// Full range: [1,FLOAT_MAX] +template +VTKM_EXEC OutType CellEdgeRatioMetric(const vtkm::IdComponent& numPts, + const PointCoordVecType& pts, + vtkm::CellShapeTagTetra, + const vtkm::exec::FunctorBase& worklet) +{ + if (numPts != 4) + { + worklet.RaiseError("Edge ratio metric(tetrahedron) requires 4 points."); + return OutType(0.0); + } + + vtkm::IdComponent numEdges = 6; //pts.GetNumberOfComponents(); + + //The 6 edges of a tetrahedron + using Edge = typename PointCoordVecType::ComponentType; + const Edge TetEdges[6] = { pts[1] - pts[0], pts[2] - pts[1], pts[0] - pts[2], + pts[3] - pts[0], pts[3] - pts[1], pts[3] - pts[2] }; + + return vtkm::exec::cellmetrics::ComputeEdgeRatio(vtkm::make_VecC(TetEdges, numEdges)); +} + + +// Compute the edge ratio of a hexahedron. +// Formula: Maximum edge length divided by minimum edge length +// Equals 1 for a unit cube +// Full range: [1,FLOAT_MAX] +template +VTKM_EXEC OutType CellEdgeRatioMetric(const vtkm::IdComponent& numPts, + const PointCoordVecType& pts, + vtkm::CellShapeTagHexahedron, + const vtkm::exec::FunctorBase& worklet) +{ + if (numPts != 8) + { + worklet.RaiseError("Edge ratio metric(hexahedron) requires 8 points."); + return OutType(0.0); + } + + vtkm::IdComponent numEdges = 12; //pts.GetNumberOfComponents(); + + //The 12 edges of a hexahedron + using Edge = typename PointCoordVecType::ComponentType; + const Edge HexEdges[12] = { pts[1] - pts[0], pts[2] - pts[1], pts[3] - pts[2], pts[0] - pts[3], + pts[5] - pts[4], pts[6] - pts[5], pts[7] - pts[6], pts[4] - pts[7], + pts[4] - pts[0], pts[5] - pts[1], pts[6] - pts[2], pts[7] - pts[3] }; + + return vtkm::exec::cellmetrics::ComputeEdgeRatio(vtkm::make_VecC(HexEdges, numEdges)); +} + + +// Compute the edge ratio of a wedge/prism. +// Formula: Maximum edge length divided by minimum edge length +// Equals 1 for a right unit wedge +// Full range: [1,FLOAT_MAX] +template +VTKM_EXEC OutType CellEdgeRatioMetric(const vtkm::IdComponent& numPts, + const PointCoordVecType& pts, + vtkm::CellShapeTagWedge, + const vtkm::exec::FunctorBase& worklet) +{ + if (numPts != 6) + { + worklet.RaiseError("Edge ratio metric(wedge) requires 6 points."); + return OutType(0.0); + } + + vtkm::IdComponent numEdges = 9; //pts.GetNumberOfComponents(); + + //The 9 edges of a wedge/prism + using Edge = typename PointCoordVecType::ComponentType; + const Edge WedgeEdges[9] = { pts[1] - pts[0], pts[2] - pts[1], pts[0] - pts[2], + pts[4] - pts[3], pts[5] - pts[4], pts[3] - pts[5], + pts[3] - pts[0], pts[4] - pts[1], pts[5] - pts[2] }; + + return vtkm::exec::cellmetrics::ComputeEdgeRatio(vtkm::make_VecC(WedgeEdges, numEdges)); +} + +// Compute the edge ratio of a pyramid. +// Formula: Maximum edge length divided by minimum edge length +// TODO: Equals 1 for a right unit (square base?) pyramid (?) +// Full range: [1,FLOAT_MAX] +// TODO: Verdict/VTK don't define this metric for a pyramid. What does VisIt output? +template +VTKM_EXEC OutType CellEdgeRatioMetric(const vtkm::IdComponent& numPts, + const PointCoordVecType& pts, + vtkm::CellShapeTagPyramid, + const vtkm::exec::FunctorBase& worklet) +{ + if (numPts != 5) + { + worklet.RaiseError("Edge ratio metric(pyramid) requires 5 points."); + return OutType(0.0); + } + + vtkm::IdComponent numEdges = 8; // pts.GetNumberOfComponents(); + + //The 8 edges of a pyramid (4 quadrilateral base edges + 4 edges to apex) + using Edge = typename PointCoordVecType::ComponentType; + const Edge PyramidEdges[8] = { + pts[1] - pts[0], pts[2] - pts[1], pts[2] - pts[3], pts[3] - pts[0], + pts[4] - pts[0], pts[4] - pts[1], pts[4] - pts[2], pts[4] - pts[3] + }; + + return vtkm::exec::cellmetrics::ComputeEdgeRatio( + vtkm::make_VecC(PyramidEdges, numEdges)); +} + + + +} // namespace cellmetrics +} // namespace exec +} // namespace vtkm + +#endif // vtk_m_exec_cellmetrics_CellEdgeRatioMetric_h diff --git a/vtkm/filter/CMakeLists.txt b/vtkm/filter/CMakeLists.txt index 39a3ce794..df7eba3d0 100644 --- a/vtkm/filter/CMakeLists.txt +++ b/vtkm/filter/CMakeLists.txt @@ -43,6 +43,7 @@ set(headers MarchingCubes.h Mask.h MaskPoints.h + MeshQuality.h NDEntropy.h NDHistogram.h OscillatorSource.h @@ -104,6 +105,7 @@ set(header_template_sources MarchingCubes.hxx Mask.hxx MaskPoints.hxx + MeshQuality.hxx NDEntropy.hxx NDHistogram.hxx OscillatorSource.hxx diff --git a/vtkm/filter/MeshQuality.h b/vtkm/filter/MeshQuality.h new file mode 100644 index 000000000..e5edfccee --- /dev/null +++ b/vtkm/filter/MeshQuality.h @@ -0,0 +1,101 @@ +//============================================================================ +// 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_MeshQuality_h +#define vtk_m_filter_MeshQuality_h + +#include "vtkm/CellShape.h" +#include "vtkm/filter/FilterCell.h" +#include "vtkm/worklet/MeshQuality.h" +#include + +namespace vtkm +{ +namespace filter +{ + +//Names of the available cell metrics, for use in +//the output dataset fields +//TODO: static const? +static const std::string MetricNames[] = { "empty", + "diagonalRatio", + "edgeRatio", + //"skew", + "oddy", + "relativeSizeSquared", + "volume" }; + +//Different cell metrics available to use +//TODO: static const? +enum class CellMetric +{ + EMPTY, //0 + DIAGONAL_RATIO, + EDGE_RATIO, + ODDY, + RELATIVE_SIZE, + VOLUME, + NUMBER_OF_CELL_METRICS //(num metrics = NUMBER_OF_CELL_METRICS - 2) +}; + +/** \brief Computes the quality of an unstructured cell-based mesh. The quality is defined in terms of the + * summary statistics (frequency, mean, variance, min, max) of metrics computed over the mesh + * cells. One of several different metrics can be specified for a given cell type, and the mesh + * can consist of one or more different cell types. The resulting mesh quality is stored as one + * or more new fields in the output dataset of this filter, with a separate field for each cell type. + * Each field contains the metric summary statistics for the cell type. + * Summary statists with all 0 values imply that the specified metric does not support the cell type. + */ +class MeshQuality : public vtkm::filter::FilterCell +{ +public: + using ShapeMetricsVecType = std::vector>; + + VTKM_CONT MeshQuality(const ShapeMetricsVecType& metrics); + + template + VTKM_CONT vtkm::cont::DataSet DoExecute( + const vtkm::cont::DataSet& input, + const vtkm::cont::ArrayHandle, StorageType>& points, + const vtkm::filter::FieldMetadata& fieldMeta, + const vtkm::filter::PolicyBase& policy); + +private: + //A user-assigned cell metric per shape/cell type + //Empty metric if not provided by user + //Length of vector is the number of different VTK-m cell types + std::vector CellTypeMetrics; +}; + + +template <> +class FilterTraits +{ +public: + using InputFieldTypeList = vtkm::TypeListTagFieldVec3; +}; + + +} // namespace filter +} // namespace vtkm + +#include + +#endif // vtk_m_filter_MeshQuality_h diff --git a/vtkm/filter/MeshQuality.hxx b/vtkm/filter/MeshQuality.hxx new file mode 100644 index 000000000..99ff29411 --- /dev/null +++ b/vtkm/filter/MeshQuality.hxx @@ -0,0 +1,250 @@ +//============================================================================ +// 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/DynamicCellSet.h" +#include "vtkm/cont/ErrorFilterExecution.h" +#include "vtkm/cont/Field.h" +#include "vtkm/filter/internal/CreateResult.h" +#include "vtkm/worklet/DispatcherMapTopology.h" + +#define DEBUG_PRINT + +namespace vtkm +{ +namespace filter +{ + +namespace debug +{ +#ifdef DEBUG_PRINT +//---------------------------------------------------------------------------- +template +void MeshQualityDebug(const vtkm::cont::ArrayHandle& outputArray, const char* name) +{ + typedef vtkm::cont::internal::Storage StorageType; + typedef typename StorageType::PortalConstType PortalConstType; + PortalConstType readPortal = outputArray.GetPortalConstControl(); + vtkm::Id numElements = readPortal.GetNumberOfValues(); + std::cout << name << "= " << numElements << " ["; + for (vtkm::Id i = 0; i < numElements; i++) + std::cout << (int)readPortal.Get(i) << " "; + std::cout << "]\n"; +} +#else +template +void MeshQualityDebug(const vtkm::cont::ArrayHandle& vtkmNotUsed(outputArray), + const char* vtkmNotUsed(name)) +{ +} +#endif +} // namespace debug + + +inline VTKM_CONT MeshQuality::MeshQuality( + const std::vector>& metrics) + : vtkm::filter::FilterCell() +{ + this->SetUseCoordinateSystemAsField(true); + this->CellTypeMetrics.assign(vtkm::NUMBER_OF_CELL_SHAPES, CellMetric::EMPTY); + for (auto p : metrics) + this->CellTypeMetrics[p.first] = p.second; +} + +template +inline VTKM_CONT vtkm::cont::DataSet MeshQuality::DoExecute( + const vtkm::cont::DataSet& input, + const vtkm::cont::ArrayHandle, StorageType>& points, + const vtkm::filter::FieldMetadata& fieldMeta, + const vtkm::filter::PolicyBase& policy) +{ + VTKM_ASSERT(fieldMeta.IsPointField()); + + using Algorithm = vtkm::cont::Algorithm; + using ShapeHandle = vtkm::cont::ArrayHandle; + using IdHandle = vtkm::cont::ArrayHandle; + using QualityWorklet = vtkm::worklet::MeshQuality; + using FieldStatsWorklet = vtkm::worklet::FieldStatistics; + + //TODO: Should other cellset types be supported? + vtkm::cont::CellSetExplicit<> cellSet; + input.GetCellSet(this->GetActiveCellSetIndex()).CopyTo(cellSet); + + ShapeHandle cellShapes = + cellSet.GetShapesArray(vtkm::TopologyElementTagPoint(), vtkm::TopologyElementTagCell()); + + //Obtain the frequency counts of each cell type in the input dataset + IdHandle cellShapeCounts; + ShapeHandle uniqueCellShapes; + ShapeHandle sortedShapes; + Algorithm::Copy(cellShapes, sortedShapes); + Algorithm::Sort(sortedShapes); + Algorithm::ReduceByKey( + sortedShapes, + vtkm::cont::make_ArrayHandleConstant(vtkm::Id(1), cellShapes.GetNumberOfValues()), + uniqueCellShapes, + cellShapeCounts, + vtkm::Add()); + + //Invoke the MeshQuality worklet + vtkm::cont::ArrayHandle outArray; + vtkm::cont::ArrayHandle cellMetrics = vtkm::cont::make_ArrayHandle(CellTypeMetrics); + vtkm::worklet::DispatcherMapTopology dispatcher; + dispatcher.Invoke( + vtkm::filter::ApplyPolicy(cellSet, policy), cellShapeCounts, cellMetrics, points, outArray); + + //Build the output dataset: a separate field for each cell type that has a specified metric + vtkm::cont::DataSet result; + result.CopyStructure(input); //clone of the input dataset + + auto cellShapePortal = cellShapes.GetPortalConstControl(); + auto metricValuesPortal = outArray.GetPortalConstControl(); + + const vtkm::Id numCells = outArray.GetNumberOfValues(); + T currMetric = 0; + vtkm::UInt8 currShape = 0; + + //Output metric values stored in separate containers + //based on shape type. Unsupported shape types in VTK-m + //are represented with an empty "placeholder" container. + std::vector> metricValsPerShape = { + { /*placeholder*/ }, { /*vertices*/ }, { /*placeholder*/ }, { /*lines*/ }, + { /*placeholder*/ }, { /*triangles*/ }, { /*placeholder*/ }, { /*polygons*/ }, + { /*placeholder*/ }, { /*quads*/ }, { /*tetrahedrons*/ }, { /*placeholder*/ }, + { /*hexahedrons*/ }, { /*wedges*/ }, { /*pyramids*/ } + }; + + for (vtkm::Id metricArrayIndex = 0; metricArrayIndex < numCells; metricArrayIndex++) + { + currShape = cellShapePortal.Get(metricArrayIndex); + currMetric = metricValuesPortal.Get(metricArrayIndex); + metricValsPerShape[currShape].emplace_back(currMetric); + } + + //Compute the mesh quality for each shape type. This consists + //of computing the summary statistics of the metric values for + //each cell of the given shape type. + const vtkm::Id numUniqueShapes = uniqueCellShapes.GetNumberOfValues(); + auto uniqueCellShapesPortal = uniqueCellShapes.GetPortalConstControl(); + auto numCellsPerShapePortal = cellShapeCounts.GetPortalConstControl(); + std::string fieldName = "", metricName = ""; + vtkm::UInt8 cellShape = 0; + vtkm::Id cellCount = 0; + bool skipShape = false; + for (vtkm::Id shapeIndex = 0; shapeIndex < numUniqueShapes; shapeIndex++) + { + cellShape = uniqueCellShapesPortal.Get(shapeIndex); + cellCount = numCellsPerShapePortal.Get(shapeIndex); + metricName = MetricNames[static_cast(CellTypeMetrics[cellShape])]; + + //Skip over shapes with an empty/unspecified metric; + //don't include a field for them + if (CellTypeMetrics[cellShape] == CellMetric::EMPTY) + continue; + + switch (cellShape) + { + case vtkm::CELL_SHAPE_EMPTY: + skipShape = true; + break; + case vtkm::CELL_SHAPE_VERTEX: + fieldName = "vertices"; + break; + case vtkm::CELL_SHAPE_LINE: + fieldName = "lines"; + break; + case vtkm::CELL_SHAPE_TRIANGLE: + fieldName = "triangles"; + break; + case vtkm::CELL_SHAPE_POLYGON: + fieldName = "polygons"; + break; + case vtkm::CELL_SHAPE_QUAD: + fieldName = "quads"; + break; + case vtkm::CELL_SHAPE_TETRA: + fieldName = "tetrahedrons"; + break; + case vtkm::CELL_SHAPE_HEXAHEDRON: + fieldName = "hexahedrons"; + break; + case vtkm::CELL_SHAPE_WEDGE: + fieldName = "wedges"; + break; + case vtkm::CELL_SHAPE_PYRAMID: + fieldName = "pyramids"; + break; + default: + skipShape = true; + break; + } + + //Skip over shapes of empty cell type; don't include a field for them + if (skipShape) + continue; + + fieldName += "-" + metricName; + auto shapeMetricVals = metricValsPerShape[cellShape]; + auto shapeMetricValsHandle = vtkm::cont::make_ArrayHandle(std::move(shapeMetricVals)); + + //Invoke the field stats worklet on the array of metric values for this shape type + typename FieldStatsWorklet::StatInfo statinfo; + FieldStatsWorklet().Run(shapeMetricValsHandle, statinfo); + + //Retrieve summary stats from the output stats struct. + //These stats define the mesh quality with respect to this shape type. + std::vector shapeMeshQuality = { + T(cellCount), statinfo.mean, statinfo.variance, statinfo.minimum, statinfo.maximum + }; + + //Append the summary stats into the output dataset as a new field + result.AddField(vtkm::cont::make_Field(fieldName, + vtkm::cont::Field::Association::CELL_SET, + "cells", + shapeMeshQuality, + vtkm::CopyFlag::On)); + + std::cout << "-----------------------------------------------------\n" + << "Mesh quality of " << fieldName << ":\n" + << "Number of cells: " << cellCount << "\n" + << "Mean: " << statinfo.mean << "\n" + << "Variance: " << statinfo.variance << "\n" + << "Minimum: " << statinfo.minimum << "\n" + << "Maximum: " << statinfo.maximum << "\n" + << "-----------------------------------------------------\n"; + } + + auto metricValsPortal = outArray.GetPortalConstControl(); + std::cout << "-----------------------------------------------------\n" + << "Metric values - all cells:\n"; + for (vtkm::Id v = 0; v < outArray.GetNumberOfValues(); v++) + std::cout << metricValsPortal.Get(v) << "\n"; + std::cout << "-----------------------------------------------------\n"; + + //Append the metric values of all cells into the output + //dataset as a new field + std::string s = "allCells-metricValues"; + result.AddField( + vtkm::cont::Field(s, vtkm::cont::Field::Association::CELL_SET, "cells", outArray)); + + return result; +} + +} // namespace filter +} // namespace vtkm diff --git a/vtkm/filter/testing/CMakeLists.txt b/vtkm/filter/testing/CMakeLists.txt index 31008f5a1..a0823d1cf 100644 --- a/vtkm/filter/testing/CMakeLists.txt +++ b/vtkm/filter/testing/CMakeLists.txt @@ -37,6 +37,7 @@ set(unit_tests UnitTestMarchingCubesFilter.cxx UnitTestMaskFilter.cxx UnitTestMaskPointsFilter.cxx + UnitTestMeshQualityFilter.cxx UnitTestMultiBlockFilters.cxx UnitTestMultiBlockHistogramFilter.cxx UnitTestNDEntropyFilter.cxx diff --git a/vtkm/filter/testing/UnitTestMeshQualityFilter.cxx b/vtkm/filter/testing/UnitTestMeshQualityFilter.cxx new file mode 100644 index 000000000..373c55b7a --- /dev/null +++ b/vtkm/filter/testing/UnitTestMeshQualityFilter.cxx @@ -0,0 +1,376 @@ +//============================================================================ +// 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 2018 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +// Copyright 2018 UT-Battelle, LLC. +// Copyright 2018 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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +//Adapted from vtkm/cont/testing/MakeTestDataSet.h +//Modified the content of the Make3DExplicitDataSetZoo() function +inline vtkm::cont::DataSet Make3DExplicitDataSet() +{ + vtkm::cont::DataSet dataSet; + vtkm::cont::DataSetBuilderExplicit dsb; + + using CoordType = vtkm::Vec; + + std::vector coords = { + { 0.00, 0.00, 0.00 }, { 1.00, 0.00, 0.00 }, { 2.00, 0.00, 0.00 }, { 0.00, 0.00, 1.00 }, + { 1.00, 0.00, 1.00 }, { 2.00, 0.00, 1.00 }, { 0.00, 1.00, 0.00 }, { 1.00, 1.00, 0.00 }, + { 2.00, 1.00, 0.00 }, { 0.00, 1.00, 1.00 }, { 1.00, 1.00, 1.00 }, { 2.00, 1.00, 1.00 }, + { 0.00, 2.00, 0.00 }, { 1.00, 2.00, 0.00 }, { 2.00, 2.00, 0.00 }, { 0.00, 2.00, 1.00 }, + { 1.00, 2.00, 1.00 }, { 2.00, 2.00, 1.00 }, { 1.00, 3.00, 1.00 }, { 2.75, 0.00, 1.00 }, + { 3.00, 0.00, 0.75 }, { 3.00, 0.25, 1.00 }, { 3.00, 1.00, 1.00 }, { 3.00, 1.00, 0.00 }, + { 2.57, 2.00, 1.00 }, { 3.00, 1.75, 1.00 }, { 3.00, 1.75, 0.75 }, { 3.00, 0.00, 0.00 }, + { 2.57, 0.42, 0.57 }, { 2.59, 1.43, 0.71 } + }; + + std::vector shapes; + std::vector numindices; + std::vector conn; + + //Construct the shapes/cells of the dataset + //This is a zoo of points, lines, polygons, and polyhedra + shapes.push_back(vtkm::CELL_SHAPE_HEXAHEDRON); + numindices.push_back(8); + conn.push_back(0); + conn.push_back(3); + conn.push_back(4); + conn.push_back(1); + conn.push_back(6); + conn.push_back(9); + conn.push_back(10); + conn.push_back(7); + + shapes.push_back(vtkm::CELL_SHAPE_HEXAHEDRON); + numindices.push_back(8); + conn.push_back(1); + conn.push_back(4); + conn.push_back(5); + conn.push_back(2); + conn.push_back(7); + conn.push_back(10); + conn.push_back(11); + conn.push_back(8); + + shapes.push_back(vtkm::CELL_SHAPE_TETRA); + numindices.push_back(4); + conn.push_back(24); + conn.push_back(26); + conn.push_back(25); + conn.push_back(29); + + shapes.push_back(vtkm::CELL_SHAPE_TETRA); + numindices.push_back(4); + conn.push_back(8); + conn.push_back(17); + conn.push_back(11); + conn.push_back(29); + + shapes.push_back(vtkm::CELL_SHAPE_PYRAMID); + numindices.push_back(5); + conn.push_back(24); + conn.push_back(17); + conn.push_back(8); + conn.push_back(23); + conn.push_back(29); + + shapes.push_back(vtkm::CELL_SHAPE_PYRAMID); + numindices.push_back(5); + conn.push_back(25); + conn.push_back(22); + conn.push_back(11); + conn.push_back(17); + conn.push_back(29); + + shapes.push_back(vtkm::CELL_SHAPE_WEDGE); + numindices.push_back(6); + conn.push_back(8); + conn.push_back(14); + conn.push_back(17); + conn.push_back(7); + conn.push_back(13); + conn.push_back(16); + + shapes.push_back(vtkm::CELL_SHAPE_WEDGE); + numindices.push_back(6); + conn.push_back(11); + conn.push_back(8); + conn.push_back(17); + conn.push_back(10); + conn.push_back(7); + conn.push_back(16); + + shapes.push_back(vtkm::CELL_SHAPE_VERTEX); + numindices.push_back(1); + conn.push_back(0); + + shapes.push_back(vtkm::CELL_SHAPE_VERTEX); + numindices.push_back(1); + conn.push_back(29); + + shapes.push_back(vtkm::CELL_SHAPE_LINE); + numindices.push_back(2); + conn.push_back(0); + conn.push_back(1); + + shapes.push_back(vtkm::CELL_SHAPE_LINE); + numindices.push_back(2); + conn.push_back(15); + conn.push_back(16); + + shapes.push_back(vtkm::CELL_SHAPE_TRIANGLE); + numindices.push_back(3); + conn.push_back(2); + conn.push_back(4); + conn.push_back(15); + + shapes.push_back(vtkm::CELL_SHAPE_TRIANGLE); + numindices.push_back(3); + conn.push_back(5); + conn.push_back(6); + conn.push_back(7); + + shapes.push_back(vtkm::CELL_SHAPE_QUAD); + numindices.push_back(4); + conn.push_back(0); + conn.push_back(3); + conn.push_back(5); + conn.push_back(2); + + shapes.push_back(vtkm::CELL_SHAPE_QUAD); + numindices.push_back(4); + conn.push_back(5); + conn.push_back(4); + conn.push_back(10); + conn.push_back(11); + + shapes.push_back(vtkm::CELL_SHAPE_POLYGON); + numindices.push_back(3); + conn.push_back(4); + conn.push_back(7); + conn.push_back(1); + + shapes.push_back(vtkm::CELL_SHAPE_POLYGON); + numindices.push_back(4); + conn.push_back(1); + conn.push_back(6); + conn.push_back(7); + conn.push_back(2); + + dataSet = dsb.Create(coords, shapes, numindices, conn, "coordinates", "cells"); + + return dataSet; +} + +template +std::vector TestMeshQualityFilter(const vtkm::cont::DataSet& input, + const std::vector& expectedVals, + T filter) +{ + std::vector errors; + vtkm::cont::DataSet output; + try + { + output = filter.Execute(input); + } + catch (vtkm::cont::ErrorExecution&) + { + errors.push_back("Error occured while executing filter. Exiting..."); + return errors; + } + + //Test the computed metric values (for all cells) and expected metric + //values for equality. + const vtkm::Id numFields = output.GetNumberOfFields(); + vtkm::cont::testing::TestEqualResult result = vtkm::cont::testing::test_equal_ArrayHandles( + vtkm::cont::make_ArrayHandle(expectedVals), output.GetField(numFields - 1).GetData()); + if (!result) + result.PushMessage(std::string("Data doesn't match")); + + return result.GetMessages(); +} + +//Either an error occurred during execution of the +//filter, or mismatches exist between the expected output +//and the computed filter output. +void CheckForErrors(const std::vector& messages) +{ + if (!messages.empty()) + { + std::cout << "FAIL\n"; + for (std::string m : messages) + std::cout << m << "\n"; + } + else + std::cout << "SUCCESS\n"; +} + +int TestMeshQuality() +{ + using FloatVec = std::vector; + using PairVec = std::vector>; + using StringVec = std::vector; + using CharVec = std::vector; + using QualityFilter = vtkm::filter::MeshQuality; + + //Test variables + vtkm::cont::DataSet input = Make3DExplicitDataSet(); + std::unique_ptr filter; + std::string metricSuffix; + StringVec fieldNames; + CharVec testedShapes; + PairVec shapeMetricPairs; + FloatVec expectedValues; + + /*************************************************** + * Test 1: Volume metric + ***************************************************/ + + std::cout << "Testing MeshQuality filter: Volume metric" + << "\n++++++++++++++++++++++++++++++++++++++++++++++++++\n"; + + //Assign a cell metric to compute for each different + //shape type that may exist in the input dataset. If no metric + //is specified for a shape type, then it is assumed to be EMPTY + //and no metric is computed. + testedShapes = { vtkm::CELL_SHAPE_TETRA, vtkm::CELL_SHAPE_HEXAHEDRON, vtkm::CELL_SHAPE_WEDGE, + vtkm::CELL_SHAPE_PYRAMID, vtkm::CELL_SHAPE_POLYGON, vtkm::CELL_SHAPE_LINE, + vtkm::CELL_SHAPE_QUAD, vtkm::CELL_SHAPE_TRIANGLE }; + shapeMetricPairs.clear(); + for (auto s : testedShapes) + shapeMetricPairs.push_back(vtkm::make_Pair(s, vtkm::filter::CellMetric::VOLUME)); + + //The ground truth metric value for each cell in the input dataset. + //These values are generated from VisIt using the equivalent pseudocolor + //mesh quality metric. + expectedValues = { 1, 1, 0.0100042, 0.0983333, 0.0732667, 0.0845833, -0.5, -0.5, 0, + 0, 1, 1, 1.5, 0.7071068, 2, 1, 0.5, 1 }; + + filter.reset(new QualityFilter(shapeMetricPairs)); + std::vector errors = + TestMeshQualityFilter(input, expectedValues, *filter); + std::cout << "Volume metric test: "; + CheckForErrors(errors); + + /*************************************************** + * Test 2: Edge Ratio metric + ***************************************************/ + + std::cout << "\nTesting MeshQuality filter: Edge Ratio metric" + << "\n++++++++++++++++++++++++++++++++++++++++++++++++++\n"; + + testedShapes = { vtkm::CELL_SHAPE_TETRA, vtkm::CELL_SHAPE_HEXAHEDRON, vtkm::CELL_SHAPE_WEDGE, + vtkm::CELL_SHAPE_PYRAMID, vtkm::CELL_SHAPE_POLYGON, vtkm::CELL_SHAPE_LINE, + vtkm::CELL_SHAPE_QUAD, vtkm::CELL_SHAPE_TRIANGLE }; + + shapeMetricPairs.clear(); + for (auto s : testedShapes) + shapeMetricPairs.push_back(vtkm::make_Pair(s, vtkm::filter::CellMetric::EDGE_RATIO)); + + expectedValues = { 1, 1, 2.55938, 1.80027, 2.59323, 1.73099, 1.41421, 1.41421, 0, + 0, 1, 1, 2.12132, 2.44949, 2, 1, 1.41421, 1.41421 }; + + filter.reset(new QualityFilter(shapeMetricPairs)); + errors = TestMeshQualityFilter(input, expectedValues, *filter); + std::cout << "Edge ratio metric test: "; + CheckForErrors(errors); + + /*************************************************** + * Test 3: Diagonal Ratio metric + ***************************************************/ + + std::cout << "Testing MeshQuality filter: Diagonal Ratio metric" + << "\n++++++++++++++++++++++++++++++++++++++++++++++++++\n"; + + testedShapes = { vtkm::CELL_SHAPE_HEXAHEDRON, vtkm::CELL_SHAPE_POLYGON, vtkm::CELL_SHAPE_QUAD }; + + shapeMetricPairs.clear(); + for (auto s : testedShapes) + shapeMetricPairs.push_back(vtkm::make_Pair(s, vtkm::filter::CellMetric::DIAGONAL_RATIO)); + + expectedValues = { 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 2.23607 }; + + filter.reset(new QualityFilter(shapeMetricPairs)); + errors = TestMeshQualityFilter(input, expectedValues, *filter); + std::cout << "Diagonal ratio metric test: "; + CheckForErrors(errors); + +#if 0 + /*************************************************** + * Test 4: Oddy metric + ***************************************************/ + + std::cout << "Testing MeshQuality filter: Oddy metric" + << "\n++++++++++++++++++++++++++++++++++++++++++++++++++\n"; + + testedShapes = {vtkm::CELL_SHAPE_HEXAHEDRON, vtkm::CELL_SHAPE_POLYGON, + vtkm::CELL_SHAPE_QUAD}; + + shapeMetricPairs.clear(); + for (auto s : testedShapes) + shapeMetricPairs.push_back(vtkm::make_Pair(s, vtkm::filter::CellMetric::ODDY)); + expectedValues = {0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 1.125, 0, + 0, 2.5}; /*(all cells, volume value)*/ + + filter.reset(new QualityFilter(shapeMetricPairs)); + errors = TestMeshQualityFilter(input, expectedValues, *filter); + std::cout << "Oddy metric test: "; + CheckForErrors(errors); + + /*************************************************** + * Test 4: Relative Size Squared metric + ***************************************************/ + + std::cout << "Testing MeshQuality filter: Relative Size Squared metric" + << "\n++++++++++++++++++++++++++++++++++++++++++++++++++\n"; + + testedShapes = {vtkm::CELL_SHAPE_HEXAHEDRON, vtkm::CELL_SHAPE_POLYGON, + vtkm::CELL_SHAPE_QUAD, vtkm::CELL_SHAPE_TRIANGLE, + vtkm::CELL_SHAPE_TETRA}; + + shapeMetricPairs.clear(); + for (auto s : testedShapes) + shapeMetricPairs.push_back(vtkm::make_Pair(s, vtkm::filter::CellMetric::RELATIVE_SIZE)); + expectedValues = {1, 1, 0.0341086, 0.303456, 0, 0, 0, + 0, 0, 0, 0, 0, 0.361898, 0.614047, 1, 1, + 0.307024, 1}; + + filter.reset(new QualityFilter(shapeMetricPairs)); + errors = TestMeshQualityFilter(input, expectedValues, *filter); + std::cout << "Relative Size Square metric test: "; + CheckForErrors(errors); +#endif + + return 0; +} + +int UnitTestMeshQualityFilter(int argc, char* argv[]) +{ + return vtkm::cont::testing::Testing::Run(TestMeshQuality, argc, argv); +} diff --git a/vtkm/worklet/CMakeLists.txt b/vtkm/worklet/CMakeLists.txt index 13fdf0d95..26a20c1e4 100644 --- a/vtkm/worklet/CMakeLists.txt +++ b/vtkm/worklet/CMakeLists.txt @@ -44,6 +44,7 @@ set(headers MaskNone.h MaskPoints.h MaskSelect.h + MeshQuality.h NDimsEntropy.h NDimsHistMarginalization.h NDimsHistogram.h diff --git a/vtkm/worklet/MeshQuality.h b/vtkm/worklet/MeshQuality.h new file mode 100644 index 000000000..b108fb29d --- /dev/null +++ b/vtkm/worklet/MeshQuality.h @@ -0,0 +1,129 @@ +//============================================================================ +// 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 2018 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +// Copyright 2018 UT-Battelle, LLC. +// Copyright 2018 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_MeshQuality_h +#define vtk_m_worklet_MeshQuality_h + +#include "vtkm/exec/CellMeasure.h" +#include "vtkm/exec/cellmetrics/CellDiagonalRatioMetric.h" +#include "vtkm/exec/cellmetrics/CellEdgeRatioMetric.h" +#include "vtkm/worklet/WorkletMapTopology.h" + +namespace vtkm +{ +namespace worklet +{ + +/** + * Worklet that computes mesh quality metric values for each cell in + * the input mesh. A metric is specified per cell type in the calling filter, + * and this metric is invoked over all cells of that cell type. An array of + * the computed metric values (one per cell) is returned as output. + */ +template +class MeshQuality : public vtkm::worklet::WorkletMapPointToCell +{ +public: + using ControlSignature = void(CellSetIn cellset, + WholeArrayIn counts, + WholeArrayIn metrics, + FieldInPoint pointCoords, + FieldOutCell metricOut); + using ExecutionSignature = void(CellShape, PointCount, _2, _3, _4, _5); + using InputDomain = _1; + + template + VTKM_EXEC void operator()(CellShapeType shape, + const vtkm::IdComponent& numPoints, + const CountsArrayType& counts, + const MetricsArrayType& metrics, + const PointCoordVecType& pts, + OutType& metricValue) const + { + vtkm::UInt8 thisId = shape.Id; + if (shape.Id == vtkm::CELL_SHAPE_POLYGON) + { + if (numPoints == 3) + thisId = vtkm::CELL_SHAPE_TRIANGLE; + else if (numPoints == 4) + thisId = vtkm::CELL_SHAPE_QUAD; + } + switch (thisId) + { + vtkmGenericCellShapeMacro( + metricValue = this->ComputeMetric( + numPoints, pts, counts.Get(shape.Id), CellShapeTag(), metrics.Get(CellShapeTag().Id))); + default: + this->RaiseError("Asked for metric of unknown cell type."); + metricValue = OutType(0.0); + } + } + +protected: + template + VTKM_EXEC OutType ComputeMetric(const vtkm::IdComponent& numPts, + const PointCoordVecType& pts, + const vtkm::Id& numShapes, + CellShapeType tag, + CellMetricType metric) const + { + UNUSED(numShapes); + constexpr vtkm::IdComponent dims = vtkm::CellTraits::TOPOLOGICAL_DIMENSIONS; + + //Only compute the metric for 2D and 3D shapes; return 0 otherwise + OutType metricValue = OutType(0.0); + if (dims > 0) + { + switch (metric) + { + case MetricTagType::DIAGONAL_RATIO: + metricValue = + vtkm::exec::cellmetrics::CellDiagonalRatioMetric(numPts, pts, tag, *this); + break; + case MetricTagType::EDGE_RATIO: + metricValue = + vtkm::exec::cellmetrics::CellEdgeRatioMetric(numPts, pts, tag, *this); + break; + case MetricTagType::VOLUME: + metricValue = vtkm::exec::CellMeasure(numPts, pts, tag, *this); + break; + case MetricTagType::EMPTY: + break; + default: + //Only call metric function if a metric is specified for this shape type + this->RaiseError("Asked for unknown metric."); + } + } + + return metricValue; + } +}; + +} // namespace worklet +} // namespace vtkm + +#endif // vtk_m_worklet_MeshQuality_h