Shuenhoy/vtk-m
1
0
Fork 0
mirror of https://gitlab.kitware.com/vtk/vtk-m synced 2024-09-08 13:23:51 +00:00
Code Issues 2 Actions Packages Projects Releases Wiki Activity

Merge topic 'add_mesh_quality'

Browse Source 
b622c7962 Fixed index out of bounds error for the cell counts array
4d61066e9 Removed all modifications in the internal device adapter algorithm header files.
06ac9f721 Revised version of the original mesh quality merge request
e54001367 Added few lines of code missing from cuda device adapter header
3dd34d251 Added custom CopyIf function
417dbcea7 Removed all modifications in the internal device adapter algorithm header files.
50cb805ce Fixed cuda device adapter alg
8c070caa0 Added few lines of code missing from cuda device adapter header
...

Acked-by: Kitware Robot <kwrobot@kitware.com>
Acked-by: Kenneth Moreland <kmorel@sandia.gov>
Merge-request: !1714
This commit is contained in:
Kenneth Moreland 2019-07-17 19:42:12 +00:00 committed by Kitware Robot
commit 79909ecf12
14 changed files with 1759 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
examples/CMakeLists.txt
View File

@ -23,6 +23,7 @@ if(VTKm_ENABLE_EXAMPLES)
add_subdirectory(histogram)
add_subdirectory(isosurface)
add_subdirectory(lagrangian)
add_subdirectory(mesh_quality)
add_subdirectory(multi_backend)
add_subdirectory(oscillator)
add_subdirectory(particle_advection)

37
examples/mesh_quality/CMakeLists.txt Normal file
View File

@ -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()

286
examples/mesh_quality/MeshQuality.cxx Normal file
View File

@ -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 <string>
#include <vector>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/DataSetBuilderExplicit.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/filter/MeshQuality.h>
#include <vtkm/io/reader/VTKDataSetReader.h>
#include <vtkm/io/writer/VTKDataSetWriter.h>
/**
* 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<vtkm::Float64, 3>;
std::vector<CoordType> 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<vtkm::UInt8> shapes;
std::vector<vtkm::IdComponent> numindices;
std::vector<vtkm::Id> 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 <typename T>
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] << " <input.vtk_file> <output.vtk_fileName>" << 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<vtkm::Pair<vtkm::UInt8, vtkm::filter::CellMetric>> 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;
}

3
vtkm/exec/CMakeLists.txt
View File

@ -57,3 +57,6 @@ add_subdirectory(cuda)
#-----------------------------------------------------------------------------
add_subdirectory(testing)
#-----------------------------------------------------------------------------
add_subdirectory(cellmetrics)

28
vtkm/exec/cellmetrics/CMakeLists.txt Normal file
View File

@ -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})

236
vtkm/exec/cellmetrics/CellDiagonalRatioMetric.h Normal file
View File

@ -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 <typename OutType, typename VecType>
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<FloatType>())
return vtkm::Infinity<OutType>();
//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<OutType>()); //normal case
return vtkm::Max(diagonalRatio, OutType(-1) * vtkm::Infinity<OutType>());
}
// ========================= Unsupported cells ==================================
// By default, cells have zero shape unless the shape type template is specialized below.
template <typename OutType, typename PointCoordVecType, typename CellShapeType>
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 <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellDiagonalRatioMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagPolygon,
const vtkm::exec::FunctorBase& worklet)
{
switch (numPts)
{
case 4:
return CellDiagonalRatioMetric<OutType>(numPts, pts, vtkm::CellShapeTagQuad(), worklet);
default:
break;
}
return OutType(-1.0);
}
*/
/*
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellDiagonalRatioMetric(const vtkm::IdComponent&,
const PointCoordVecType&,
vtkm::CellShapeTagLine,
const vtkm::exec::FunctorBase& worklet)
{
UNUSED(worklet);
return OutType(-1.0);
}
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellDiagonalRatioMetric(const vtkm::IdComponent&,
const PointCoordVecType&,
vtkm::CellShapeTagTriangle,
const vtkm::exec::FunctorBase& worklet)
{
UNUSED(worklet);
return OutType(-1.0);
}
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellDiagonalRatioMetric(const vtkm::IdComponent&,
const PointCoordVecType&,
vtkm::CellShapeTagTetra,
const vtkm::exec::FunctorBase& worklet)
{
UNUSED(worklet);
return OutType(-1.0);
}
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellDiagonalRatioMetric(const vtkm::IdComponent&,
const PointCoordVecType&,
vtkm::CellShapeTagWedge,
const vtkm::exec::FunctorBase& worklet)
{
UNUSED(worklet);
return OutType(-1.0);
}
template <typename OutType, typename PointCoordVecType>
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 <typename OutType, typename PointCoordVecType>
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<OutType>(
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 <typename OutType, typename PointCoordVecType>
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<OutType>(
vtkm::make_VecC(HexDiagonals, numDiagonals));
}
} // namespace cellmetrics
} // namespace exec
} // namespace vtkm
#endif // vtk_m_exec_cellmetrics_CellEdgeRatioMetric_h

299
vtkm/exec/cellmetrics/CellEdgeRatioMetric.h Normal file
View File

@ -0,0 +1,299 @@
//============================================================================
// 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 <typename OutType, typename VecType>
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<FloatType>())
return vtkm::Infinity<OutType>();
//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<OutType>()); //normal case
return vtkm::Max(edgeRatio, OutType(-1) * vtkm::Infinity<OutType>());
}
// ========================= Unsupported cells ==================================
// By default, cells have zero shape unless the shape type template is specialized below.
template <typename OutType, typename PointCoordVecType, typename CellShapeType>
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 <typename OutType, typename PointCoordVecType>
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 <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellEdgeRatioMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagTriangle,
const vtkm::exec::FunctorBase& worklet)
{
if (numPts != 3)
{
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<OutType>(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 <typename OutType, typename PointCoordVecType>
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<OutType>(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 <typename OutType, typename PointCoordVecType>
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<OutType>(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 <typename OutType, typename PointCoordVecType>
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<OutType>(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 <typename OutType, typename PointCoordVecType>
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<OutType>(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 <typename OutType, typename PointCoordVecType>
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<OutType>(
vtkm::make_VecC(PyramidEdges, numEdges));
}
} // namespace cellmetrics
} // namespace exec
} // namespace vtkm
#endif // vtk_m_exec_cellmetrics_CellEdgeRatioMetric_h

2
vtkm/filter/CMakeLists.txt
View File

@ -43,6 +43,7 @@ set(headers
MarchingCubes.h
Mask.h
MaskPoints.h
MeshQuality.h
NDEntropy.h
NDHistogram.h
OscillatorSource.h
@ -105,6 +106,7 @@ set(header_template_sources
MarchingCubes.hxx
Mask.hxx
MaskPoints.hxx
MeshQuality.hxx
NDEntropy.hxx
NDHistogram.hxx
OscillatorSource.hxx

101
vtkm/filter/MeshQuality.h Normal file
View File

@ -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 <vtkm/worklet/FieldStatistics.h>
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<MeshQuality>
{
public:
using ShapeMetricsVecType = std::vector<vtkm::Pair<vtkm::UInt8, CellMetric>>;
VTKM_CONT MeshQuality(const ShapeMetricsVecType& metrics);
template <typename T, typename StorageType, typename DerivedPolicy>
VTKM_CONT vtkm::cont::DataSet DoExecute(
const vtkm::cont::DataSet& input,
const vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>, StorageType>& points,
const vtkm::filter::FieldMetadata& fieldMeta,
const vtkm::filter::PolicyBase<DerivedPolicy>& 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<CellMetric> CellTypeMetrics;
};
template <>
class FilterTraits<vtkm::filter::MeshQuality>
{
public:
using InputFieldTypeList = vtkm::TypeListTagFieldVec3;
};
} // namespace filter
} // namespace vtkm
#include <vtkm/filter/MeshQuality.hxx>
#endif // vtk_m_filter_MeshQuality_h

258
vtkm/filter/MeshQuality.hxx Normal file
View File

@ -0,0 +1,258 @@
//============================================================================
// 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 <typename T, typename S = vtkm::cont::DeviceAdapterId>
void MeshQualityDebug(const vtkm::cont::ArrayHandle<T, S>& outputArray, const char* name)
{
typedef vtkm::cont::internal::Storage<T, S> 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 <typename T, typename S>
void MeshQualityDebug(const vtkm::cont::ArrayHandle<T, S>& vtkmNotUsed(outputArray),
const char* vtkmNotUsed(name))
{
}
#endif
} // namespace debug
inline VTKM_CONT MeshQuality::MeshQuality(
const std::vector<vtkm::Pair<vtkm::UInt8, CellMetric>>& metrics)
: vtkm::filter::FilterCell<MeshQuality>()
{
this->SetUseCoordinateSystemAsField(true);
this->CellTypeMetrics.assign(vtkm::NUMBER_OF_CELL_SHAPES, CellMetric::EMPTY);
for (auto p : metrics)
this->CellTypeMetrics[p.first] = p.second;
}
template <typename T, typename StorageType, typename DerivedPolicy>
inline VTKM_CONT vtkm::cont::DataSet MeshQuality::DoExecute(
const vtkm::cont::DataSet& input,
const vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>, StorageType>& points,
const vtkm::filter::FieldMetadata& fieldMeta,
const vtkm::filter::PolicyBase<DerivedPolicy>& policy)
{
VTKM_ASSERT(fieldMeta.IsPointField());
using Algorithm = vtkm::cont::Algorithm;
using ShapeHandle = vtkm::cont::ArrayHandle<vtkm::UInt8>;
using IdHandle = vtkm::cont::ArrayHandle<vtkm::Id>;
using QualityWorklet = vtkm::worklet::MeshQuality<CellMetric>;
using FieldStatsWorklet = vtkm::worklet::FieldStatistics<T>;
//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 uniqueCellCounts;
ShapeHandle uniqueCellShapes, sortedShapes;
Algorithm::Copy(cellShapes, sortedShapes);
Algorithm::Sort(sortedShapes);
Algorithm::ReduceByKey(
sortedShapes,
vtkm::cont::make_ArrayHandleConstant(vtkm::Id(1), cellShapes.GetNumberOfValues()),
uniqueCellShapes,
uniqueCellCounts,
vtkm::Add());
std::cout << "uniqueCellCounts: " << uniqueCellCounts.GetNumberOfValues() << "\n";
const vtkm::Id numUniqueShapes = uniqueCellShapes.GetNumberOfValues();
auto uniqueCellShapesPortal = uniqueCellShapes.GetPortalConstControl();
auto numCellsPerShapePortal = uniqueCellCounts.GetPortalConstControl();
std::vector<vtkm::Id> tempCounts(vtkm::NUMBER_OF_CELL_SHAPES);
for (vtkm::Id i = 0; i < numUniqueShapes; i++)
tempCounts[uniqueCellShapesPortal.Get(i)] = numCellsPerShapePortal.Get(i);
IdHandle cellShapeCounts = vtkm::cont::make_ArrayHandle(tempCounts);
std::cout << "cellShapeCounts: " << cellShapeCounts.GetNumberOfValues() << "\n";
//Invoke the MeshQuality worklet
vtkm::cont::ArrayHandle<T> outArray;
vtkm::cont::ArrayHandle<CellMetric> cellMetrics = vtkm::cont::make_ArrayHandle(CellTypeMetrics);
std::cout << "cellMetrics: " << cellMetrics.GetNumberOfValues() << "\n";
vtkm::worklet::DispatcherMapTopology<QualityWorklet> 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<std::vector<T>> 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.
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<vtkm::UInt8>(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<T> 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

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

@ -37,6 +37,7 @@ set(unit_tests
UnitTestMarchingCubesFilter.cxx
UnitTestMaskFilter.cxx
UnitTestMaskPointsFilter.cxx
UnitTestMeshQualityFilter.cxx
UnitTestMultiBlockFilters.cxx
UnitTestMultiBlockHistogramFilter.cxx
UnitTestNDEntropyFilter.cxx

376
vtkm/filter/testing/UnitTestMeshQualityFilter.cxx Normal file
View File

@ -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 <string>
#include <typeinfo>
#include <vector>
#include <vtkm/cont/ArrayRangeCompute.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/DataSetBuilderExplicit.h>
#include <vtkm/cont/DeviceAdapterAlgorithm.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/filter/MeshQuality.h>
#include <vtkm/io/reader/VTKDataSetReader.h>
#include <vtkm/io/writer/VTKDataSetWriter.h>
//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<vtkm::Float64, 3>;
std::vector<CoordType> 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<vtkm::UInt8> shapes;
std::vector<vtkm::IdComponent> numindices;
std::vector<vtkm::Id> 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 <typename T>
std::vector<std::string> TestMeshQualityFilter(const vtkm::cont::DataSet& input,
const std::vector<vtkm::Float64>& expectedVals,
T filter)
{
std::vector<std::string> 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<std::string>& 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<vtkm::Float64>;
using PairVec = std::vector<vtkm::Pair<vtkm::UInt8, vtkm::filter::CellMetric>>;
using StringVec = std::vector<std::string>;
using CharVec = std::vector<vtkm::UInt8>;
using QualityFilter = vtkm::filter::MeshQuality;
//Test variables
vtkm::cont::DataSet input = Make3DExplicitDataSet();
std::unique_ptr<QualityFilter> 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<std::string> errors =
TestMeshQualityFilter<QualityFilter>(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<QualityFilter>(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<QualityFilter>(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<QualityFilter>(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<QualityFilter>(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);
}

1
vtkm/worklet/CMakeLists.txt
View File

@ -44,6 +44,7 @@ set(headers
MaskNone.h
MaskPoints.h
MaskSelect.h
MeshQuality.h
NDimsEntropy.h
NDimsHistMarginalization.h
NDimsHistogram.h

130
vtkm/worklet/MeshQuality.h Normal file
View File

@ -0,0 +1,130 @@
//============================================================================
// 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 <typename MetricTagType>
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 <typename CellShapeType,
typename PointCoordVecType,
typename CountsArrayType,
typename MetricsArrayType,
typename OutType>
VTKM_EXEC void operator()(CellShapeType shape,
const vtkm::IdComponent& numPoints,
const CountsArrayType& counts,
const MetricsArrayType& metrics,
const PointCoordVecType& pts,
OutType& metricValue) const
{
printf("shape.Id: %u\n", shape.Id);
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<OutType>(
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 <typename OutType,
typename PointCoordVecType,
typename CellShapeType,
typename CellMetricType>
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<CellShapeType>::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<OutType>(numPts, pts, tag, *this);
break;
case MetricTagType::EDGE_RATIO:
metricValue =
vtkm::exec::cellmetrics::CellEdgeRatioMetric<OutType>(numPts, pts, tag, *this);
break;
case MetricTagType::VOLUME:
metricValue = vtkm::exec::CellMeasure<OutType>(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