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

overhaul workings of mesh quality filter. Unit test a little rocky right now, but heading in the right direction

Browse Source
This commit is contained in:
Hank Childs 2019-09-16 17:19:18 -07:00
parent b5193b2708
commit ed15f684ea
5 changed files with 350 additions and 251 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

278
vtkm/exec/cellmetrics/CellJacobianMetric.h Normal file
View File

@ -0,0 +1,278 @@
//============================================================================
// 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_Jacobian_h
#define vtk_m_exec_cellmetrics_Jacobian_h
/*
* Mesh quality metric functions that computes the jacobian of mesh cells.
* The jacobian of a cell is defined as the determinant of the Jociabian matrix
*
* 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.
*
* 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 CellJacobianMetricOfQuad(const VecType& edgeCalculations,
const VecType& axes)
{
const vtkm::Id numCalculations = edgeCalculations.GetNumberOfComponents();
//Compare partitions of quad to find min
using axesType = typename VecType::ComponentType;
axesType centerCalculation = vtkm::Cross(axes[0], axes[1]);
vtkm::Normalize(centerCalculation);
OutType currCalculation, minCalculation = vtkm::Infinity<OutType>();
for (vtkm::IdComponent i = 0; i < numCalculations; i++)
{
currCalculation = vtkm::Dot(edgeCalculations[i], centerCalculation);
if (currCalculation < minCalculation)
minCalculation = currCalculation;
}
if (minCalculation > 0)
return vtkm::Min(minCalculation, vtkm::Infinity<OutType>()); //normal case
return vtkm::Max(minCalculation, OutType(-1) * vtkm::Infinity<OutType>());
}
template <typename OutType, typename VecType>
VTKM_EXEC inline OutType CellJacobianMetricOfHex(const VecType& matrices)
{
const vtkm::IdComponent numMatrices = matrices.GetNumberOfComponents();
//Compare determinants to find min
OutType currDeterminant, minDeterminant;
//principle axes matrix computed outside of for loop to avoid un-necessary if statement
minDeterminant =
(OutType)vtkm::Dot(matrices[numMatrices - 1][0],
vtkm::Cross(matrices[numMatrices - 1][1], matrices[numMatrices - 1][2]));
minDeterminant /= 64.0;
for (vtkm::IdComponent i = 0; i < numMatrices - 1; i++)
{
currDeterminant =
(OutType)vtkm::Dot(matrices[i][0], vtkm::Cross(matrices[i][1], matrices[i][2]));
if (currDeterminant < minDeterminant)
minDeterminant = currDeterminant;
}
if (minDeterminant > 0)
return vtkm::Min(minDeterminant, vtkm::Infinity<OutType>()); //normal case
return vtkm::Max(minDeterminant, vtkm::NegativeInfinity<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 CellJacobianMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
CellShapeType shape,
const vtkm::exec::FunctorBase&)
{
UNUSED(numPts);
UNUSED(pts);
UNUSED(shape);
return OutType(0.0);
}
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellJacobianMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagPolygon,
const vtkm::exec::FunctorBase& worklet)
{
switch (numPts)
{
case 4:
return CellJacobianMetric<OutType>(numPts, pts, vtkm::CellShapeTagQuad(), worklet);
default:
break;
}
return OutType(-1.0);
}
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellJacobianMetric(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 CellJacobianMetric(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 CellJacobianMetric(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 CellJacobianMetric(const vtkm::IdComponent&,
const PointCoordVecType&,
vtkm::CellShapeTagPyramid,
const vtkm::exec::FunctorBase& worklet)
{
UNUSED(worklet);
return OutType(-1.0);
}
// ========================= 2D cells ==================================
// Compute the jacobian of a quadrilateral.
// Formula: min{Jacobian at each vertex}
// Equals 1 for a unit square
// Acceptable range: [0,FLOAT_MAX]
// Normal range: [0,FLOAT_MAX]
// Full range: [FLOAT_MIN,FLOAT_MAX]
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellJacobianMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagQuad,
const vtkm::exec::FunctorBase& worklet)
{
if (numPts != 4)
{
worklet.RaiseError("Jacobian metric(quad) requires 4 points.");
return OutType(0.0);
}
//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] };
const Edge QuadAxes[2] = { QuadEdges[0] - (pts[2] - pts[3]), QuadEdges[1] - (pts[3] - pts[0]) };
const Edge QuadEdgesToUse[4] = { vtkm::Cross(QuadEdges[3], QuadEdges[0]),
vtkm::Cross(QuadEdges[0], QuadEdges[1]),
vtkm::Cross(QuadEdges[1], QuadEdges[2]),
vtkm::Cross(QuadEdges[2], QuadEdges[3]) };
return vtkm::exec::cellmetrics::CellJacobianMetricOfQuad<OutType>(
vtkm::make_VecC(QuadEdgesToUse, 4), vtkm::make_VecC(QuadAxes, 2));
}
// ============================= 3D Volume cells ==================================
// Compute the jacobian of a hexahedron.
// Formula: min{ {Alpha_i}, Alpha_8*/64}
// -Alpha_i -> jacobian determinant at respective vertex
// -Alpha_8 -> jacobian at center
// Equals 1 for a unit cube
// Acceptable Range: [0, FLOAT_MAX]
// Normal Range: [0, FLOAT_MAX]
// Full range: [FLOAT_MIN ,FLOAT_MAX]
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellJacobianMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagHexahedron,
const vtkm::exec::FunctorBase& worklet)
{
if (numPts != 8)
{
worklet.RaiseError("Jacobian metric(hexahedron) requires 8 points.");
return OutType(0.0);
}
//The 12 edges of a hexahedron
using Edge = typename PointCoordVecType::ComponentType;
Edge HexEdges[12] = { pts[1] - pts[0], pts[2] - pts[1], pts[3] - pts[2], pts[3] - pts[0],
pts[4] - pts[0], pts[5] - pts[1], pts[6] - pts[2], pts[7] - pts[3],
pts[5] - pts[4], pts[6] - pts[5], pts[7] - pts[6], pts[7] - pts[4] };
Edge principleXAxis = HexEdges[0] + (pts[2] - pts[3]) + HexEdges[8] + (pts[6] - pts[7]);
Edge principleYAxis = HexEdges[3] + HexEdges[1] + HexEdges[11] + HexEdges[9];
Edge principleZAxis = HexEdges[5] + HexEdges[6] + HexEdges[7] + HexEdges[8];
const Edge hexMatrices[9][3] = { { HexEdges[0], HexEdges[3], HexEdges[4] },
{ HexEdges[1], -1 * HexEdges[0], HexEdges[5] },
{ HexEdges[2], -1 * HexEdges[1], HexEdges[6] },
{ -1 * HexEdges[3], -1 * HexEdges[2], HexEdges[7] },
{ HexEdges[11], HexEdges[8], -1 * HexEdges[4] },
{ -1 * HexEdges[8], HexEdges[9], -1 * HexEdges[5] },
{ -1 * HexEdges[9], HexEdges[10], -1 * HexEdges[6] },
{ -1 * HexEdges[10], -1 * HexEdges[11], -1 * HexEdges[7] },
{ principleXAxis, principleYAxis, principleZAxis } };
return vtkm::exec::cellmetrics::CellJacobianMetricOfHex<OutType>(
vtkm::make_VecC(hexMatrices, 12));
}
// Compute the jacobian of a tetrahedron.
// Formula: (L2 * L0) * L3
// Equals Sqrt(2) / 2 for unit equilateral tetrahedron
// Acceptable Range: [0, FLOAT_MAX]
// Normal Range: [0, FLOAT_MAX]
// Full range: [FLOAT_MIN,FLOAT_MAX]
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellJacobianMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagTetra,
const vtkm::exec::FunctorBase& worklet)
{
if (numPts != 4)
{
worklet.RaiseError("Jacobian metric requires 4 points");
return OutType(0.0);
}
//the 3 edges we need
using Edge = typename PointCoordVecType::ComponentType;
const Edge EdgesNeeded[3] = { pts[1] - pts[0], pts[0] - pts[2], pts[3] - pts[0] };
return (OutType)vtkm::Dot(vtkm::Cross(EdgesNeeded[1], EdgesNeeded[0]), EdgesNeeded[2]);
}
} // namespace cellmetrics
} // namespace exec
} // namespace vtkm
#endif // vtk_m_exec_cellmetrics_CellEdgeRatioMetric_h

41
vtkm/filter/MeshQuality.h
View File

@ -33,26 +33,35 @@ 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" };
static const std::string MetricNames[] = {
"area", "aspectGamma", "aspectRatio", "condition", "diagonalRatio", "jacobian",
"minAngle", "maxAngle", "oddy", "relativeSize", "scaledJacobian", "shape",
"shear", "skew", "stretch", "taper", "volume", "warpage"
};
//Different cell metrics available to use
//TODO: static const?
enum class CellMetric
{
EMPTY, //0
AREA,
ASPECT_GAMMA,
ASPECT_RATIO,
CONDITION,
DIAGONAL_RATIO,
EDGE_RATIO,
JACOBIAN,
MIN_ANGLE,
MAX_ANGLE,
ODDY,
RELATIVE_SIZE,
SCALED_JACOBIAN,
SHAPE,
SHEAR,
SKEW,
STRETCH,
TAPER,
VOLUME,
NUMBER_OF_CELL_METRICS //(num metrics = NUMBER_OF_CELL_METRICS - 2)
WARPAGE,
NUMBER_OF_CELL_METRICS,
EMPTY
};
/** \brief Computes the quality of an unstructured cell-based mesh. The quality is defined in terms of the
@ -68,9 +77,7 @@ class MeshQuality : public vtkm::filter::FilterCell<MeshQuality>
public:
using SupportedTypes = vtkm::TypeListTagFieldVec3;
using ShapeMetricsVecType = std::vector<vtkm::Pair<vtkm::UInt8, CellMetric>>;
VTKM_CONT MeshQuality(const ShapeMetricsVecType& metrics);
VTKM_CONT MeshQuality(CellMetric);
template <typename T, typename StorageType, typename DerivedPolicy>
VTKM_CONT vtkm::cont::DataSet DoExecute(
@ -81,9 +88,7 @@ public:
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;
CellMetric myMetric;
};
} // namespace filter

178
vtkm/filter/MeshQuality.hxx
View File

@ -56,14 +56,11 @@ void MeshQualityDebug(const vtkm::cont::ArrayHandle<T, S>& vtkmNotUsed(outputArr
} // namespace debug
inline VTKM_CONT MeshQuality::MeshQuality(
const std::vector<vtkm::Pair<vtkm::UInt8, CellMetric>>& metrics)
inline VTKM_CONT MeshQuality::MeshQuality(CellMetric metric)
: 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;
myMetric = metric;
}
template <typename T, typename StorageType, typename DerivedPolicy>
@ -75,183 +72,18 @@ inline VTKM_CONT vtkm::cont::DataSet MeshQuality::DoExecute(
{
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().CopyTo(cellSet);
ShapeHandle cellShapes =
cellSet.GetShapesArray(vtkm::TopologyElementTagCell(), vtkm::TopologyElementTagPoint());
//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());
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);
//Invoke the MeshQuality worklet
vtkm::cont::ArrayHandle<T> outArray;
vtkm::cont::ArrayHandle<CellMetric> cellMetrics = vtkm::cont::make_ArrayHandle(CellTypeMetrics);
this->Invoke(QualityWorklet{},
vtkm::filter::ApplyPolicyCellSet(cellSet, policy),
cellShapeCounts,
cellMetrics,
points,
outArray);
vtkm::worklet::MeshQuality<CellMetric> qualityWorklet;
qualityWorklet.SetMetric(myMetric);
this->Invoke(qualityWorklet, vtkm::filter::ApplyPolicyCellSet(cellSet, policy), 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.
vtkm::cont::ArrayHandle<T> shapeMeshQuality;
shapeMeshQuality.Allocate(5);
{
auto portal = shapeMeshQuality.GetPortalControl();
portal.Set(0, T(cellCount));
portal.Set(1, statinfo.mean);
portal.Set(2, statinfo.variance);
portal.Set(3, statinfo.minimum);
portal.Set(4, statinfo.maximum);
}
//Append the summary stats into the output dataset as a new field
result.AddField(vtkm::cont::make_FieldCell(fieldName, shapeMeshQuality));
#ifdef DEBUG_PRINT
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";
#endif
}
#ifdef DEBUG_PRINT
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";
#endif
//Append the metric values of all cells into the output
//dataset as a new field
const std::string s = "allCells-metricValues";

71
vtkm/filter/testing/UnitTestMeshQualityFilter.cxx
View File

@ -17,6 +17,7 @@
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#include <stdio.h>
#include <string>
#include <typeinfo>
#include <vector>
@ -189,17 +190,21 @@ inline vtkm::cont::DataSet Make3DExplicitDataSet()
template <typename T>
std::vector<std::string> TestMeshQualityFilter(const vtkm::cont::DataSet& input,
const std::vector<vtkm::Float64>& expectedVals,
const std::vector<vtkm::FloatDefault>& expectedVals,
T filter)
{
fprintf(stderr, "In TMQF\n");
std::vector<std::string> errors;
vtkm::cont::DataSet output;
try
{
fprintf(stderr, "Try execute\n");
output = filter.Execute(input);
fprintf(stderr, "Complete execute\n");
}
catch (vtkm::cont::ErrorExecution&)
{
fprintf(stderr, "Catch error\n");
errors.push_back("Error occured while executing filter. Exiting...");
return errors;
}
@ -232,7 +237,7 @@ void CheckForErrors(const std::vector<std::string>& messages)
int TestMeshQuality()
{
using FloatVec = std::vector<vtkm::Float64>;
using FloatVec = std::vector<vtkm::FloatDefault>;
using PairVec = std::vector<vtkm::Pair<vtkm::UInt8, vtkm::filter::CellMetric>>;
using StringVec = std::vector<std::string>;
using CharVec = std::vector<vtkm::UInt8>;
@ -254,75 +259,51 @@ int TestMeshQuality()
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));
filter.reset(new QualityFilter(vtkm::filter::CellMetric::VOLUME));
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
* Test 2: 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));
filter.reset(new QualityFilter(vtkm::filter::CellMetric::DIAGONAL_RATIO));
errors = TestMeshQualityFilter<QualityFilter>(input, expectedValues, *filter);
std::cout << "Diagonal ratio metric test: ";
CheckForErrors(errors);
/***************************************************
* Test 3: Jacobian metric
***************************************************/
std::cout << "Testing MeshQuality filter: Jacobian metric"
<< "\n++++++++++++++++++++++++++++++++++++++++++++++++++\n";
expectedValues = { 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 2.23607 };
filter.reset(new QualityFilter(vtkm::filter::CellMetric::JACOBIAN));
errors = TestMeshQualityFilter<QualityFilter>(input, expectedValues, *filter);
std::cout << "Jacobian metric test: ";
CheckForErrors(errors);
#if 0
/***************************************************
* Test 4: Oddy metric
* Test 5: Oddy metric
***************************************************/
std::cout << "Testing MeshQuality filter: Oddy metric"

33
vtkm/worklet/MeshQuality.h
View File

@ -24,6 +24,7 @@
#include "vtkm/exec/CellMeasure.h"
#include "vtkm/exec/cellmetrics/CellDiagonalRatioMetric.h"
#include "vtkm/exec/cellmetrics/CellEdgeRatioMetric.h"
#include "vtkm/exec/cellmetrics/CellJacobianMetric.h"
#include "vtkm/worklet/WorkletMapTopology.h"
namespace vtkm
@ -42,26 +43,22 @@ class MeshQuality : public vtkm::worklet::WorkletVisitCellsWithPoints
{
public:
using ControlSignature = void(CellSetIn cellset,
WholeArrayIn counts,
WholeArrayIn metrics,
FieldInPoint pointCoords,
FieldOutCell metricOut);
using ExecutionSignature = void(CellShape, PointCount, _2, _3, _4, _5);
using ExecutionSignature = void(CellShape, PointCount, _2, _3);
using InputDomain = _1;
template <typename CellShapeType,
typename PointCoordVecType,
typename CountsArrayType,
typename MetricsArrayType,
typename OutType>
void SetMetric(MetricTagType m) { metric = m; }
template <typename CellShapeType, typename PointCoordVecType, typename OutType>
VTKM_EXEC void operator()(CellShapeType shape,
const vtkm::IdComponent& numPoints,
const CountsArrayType& counts,
const MetricsArrayType& metrics,
//const CountsArrayType& counts,
//const MetricsArrayType& metrics,
//MetricTagType metric,
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)
{
@ -73,8 +70,7 @@ public:
switch (thisId)
{
vtkmGenericCellShapeMacro(
metricValue = this->ComputeMetric<OutType>(
numPoints, pts, counts.Get(shape.Id), CellShapeTag(), metrics.Get(CellShapeTag().Id)));
metricValue = this->ComputeMetric<OutType>(numPoints, pts, CellShapeTag(), metric));
default:
this->RaiseError("Asked for metric of unknown cell type.");
metricValue = OutType(0.0);
@ -82,17 +78,18 @@ public:
}
protected:
// data member
MetricTagType metric;
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
@ -105,10 +102,16 @@ protected:
metricValue =
vtkm::exec::cellmetrics::CellDiagonalRatioMetric<OutType>(numPts, pts, tag, *this);
break;
/* DEPRECATING
case MetricTagType::EDGE_RATIO:
metricValue =
vtkm::exec::cellmetrics::CellEdgeRatioMetric<OutType>(numPts, pts, tag, *this);
break;
*/
case MetricTagType::JACOBIAN:
metricValue =
vtkm::exec::cellmetrics::CellJacobianMetric<OutType>(numPts, pts, tag, *this);
break;
case MetricTagType::VOLUME:
metricValue = vtkm::exec::CellMeasure<OutType>(numPts, pts, tag, *this);
break;