Commit 0eb767dc authored by Hank's avatar Hank Committed by Kitware Robot

Merge topic 'jacobian-cellmetric'

f9d02b37 remove file that should not have been checked in
d5991747 rename half variable to avoid shadowing CUDA global
954e2587 fix problem with filename
91f7e763 fix copyright, warnings, cmake for Jacobian-cellmetric branch
e040ea05 formatting suggestions from KMorel
f4edb6a8 remove trailing whitespace
ccedb5aa Jacobian working, unit test working
560aa1f6 add support for outputname
...
Acked-by: Kitware Robot's avatarKitware Robot <kwrobot@kitware.com>
Merge-request: !1855
parents bde77265 f9d02b37
......@@ -21,6 +21,11 @@
set(headers
CellDiagonalRatioMetric.h
CellEdgeRatioMetric.h
CellJacobianMetric.h
TypeOfCellHexahedral.h
TypeOfCellQuadrilateral.h
TypeOfCellTetrahedral.h
TypeOfCellTriangle.h
)
......
//============================================================================
// 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.
*
* 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 "TypeOfCellHexahedral.h"
#include "TypeOfCellQuadrilateral.h"
#include "TypeOfCellTetrahedral.h"
#include "TypeOfCellTriangle.h"
#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
{
// ========================= Unsupported cells ==================================
// By default, cells return the metric 0.0 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);
}
// ========================= 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);
}
using Scalar = OutType;
using CollectionOfPoints = PointCoordVecType;
using Vector = typename PointCoordVecType::ComponentType;
const Scalar alpha0 = GetQuadAlpha0<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar alpha1 = GetQuadAlpha1<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar alpha2 = GetQuadAlpha2<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar alpha3 = GetQuadAlpha3<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar q = vtkm::Min(alpha0, vtkm::Min(alpha1, vtkm::Min(alpha2, alpha3)));
return q;
}
// ============================= 3D Volume cells ==================================
// Compute the Jacobian of a hexahedron.
// Formula: min{ {Alpha_i for i in 1..7}, 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);
}
using Scalar = OutType;
using CollectionOfPoints = PointCoordVecType;
using Vector = typename PointCoordVecType::ComponentType;
const Scalar alpha0 = GetHexAlphai<Scalar, Vector, CollectionOfPoints>(pts, vtkm::Id(0));
const Scalar alpha1 = GetHexAlphai<Scalar, Vector, CollectionOfPoints>(pts, vtkm::Id(1));
const Scalar alpha2 = GetHexAlphai<Scalar, Vector, CollectionOfPoints>(pts, vtkm::Id(2));
const Scalar alpha3 = GetHexAlphai<Scalar, Vector, CollectionOfPoints>(pts, vtkm::Id(3));
const Scalar alpha4 = GetHexAlphai<Scalar, Vector, CollectionOfPoints>(pts, vtkm::Id(4));
const Scalar alpha5 = GetHexAlphai<Scalar, Vector, CollectionOfPoints>(pts, vtkm::Id(5));
const Scalar alpha6 = GetHexAlphai<Scalar, Vector, CollectionOfPoints>(pts, vtkm::Id(6));
const Scalar alpha7 = GetHexAlphai<Scalar, Vector, CollectionOfPoints>(pts, vtkm::Id(7));
const Scalar alpha8 = GetHexAlphai<Scalar, Vector, CollectionOfPoints>(pts, vtkm::Id(8));
const Scalar alpha8Div64 = alpha8 / Scalar(64.0);
const Scalar q = vtkm::Min(
alpha0,
vtkm::Min(
alpha1,
vtkm::Min(
alpha2,
vtkm::Min(
alpha3,
vtkm::Min(alpha4,
vtkm::Min(alpha5, vtkm::Min(alpha6, vtkm::Min(alpha7, alpha8Div64))))))));
return q;
}
// Compute the Jacobian of a tetrahedron.
// Formula: (L2 x 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 (tetra) requires 4 points");
return OutType(0.0);
}
using Scalar = OutType;
using CollectionOfPoints = PointCoordVecType;
using Vector = typename PointCoordVecType::ComponentType;
const Vector L0 = GetTetraL0<Scalar, Vector, CollectionOfPoints>(pts);
const Vector L2 = GetTetraL2<Scalar, Vector, CollectionOfPoints>(pts);
const Vector L3 = GetTetraL3<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar q = vtkm::Dot(vtkm::Cross(L2, L0), L3);
return q;
}
} // namespace cellmetrics
} // namespace exec
} // namespace vtkm
#endif // vtk_m_exec_cellmetrics_CellJacobianMetric_h
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
//============================================================================
// 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_TypeOfCellTriangle
#define vtk_m_exec_cellmetrics_TypeOfCellTriangle
/**
* The Verdict manual defines a set of commonly
* used components of a triangle. For example,
* area, side lengths, and so forth.
*
* These definitions can be found starting on
* page 17 of the Verdict manual.
*
* This file contains a set of functions which
* implement return the values of those commonly
* used components for subsequent use in metrics.
*/
/**
* Returns the L0 vector, as defined by the verdict manual.
*
* \param [in] pts The three points which define the triangle.
* \return Returns the vector.
*/
template <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Vector GetTriangleL0(const CollectionOfPoints& pts)
{
const Vector L0(pts[2] - pts[1]);
return L0;
}
/**
* Returns the L1 vector, as defined by the verdict manual.
*
* \param [in] pts The three points which define the triangle.
* \return Returns the vector.
*/
template <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Vector GetTriangleL1(const CollectionOfPoints& pts)
{
const Vector L1(pts[0] - pts[2]);
return L1;
}
/**
* Returns the L2 vector, as defined by the verdict manual.
*
* \param [in] pts The three points which define the triangle.
* \return Returns the vector.
*/
template <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Vector GetTriangleL2(const CollectionOfPoints& pts)
{
const Vector L2(pts[1] - pts[0]);
return L2;
}
/**
* Returns the L0 vector's magnitude, as defined by the verdict manual.
*
* \param [in] pts The three points which define the triangle.
* \return Returns the magnitude of the vector.
*/
template <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetTriangleL0Magnitude(const CollectionOfPoints& pts)
{
const Scalar l0 =
vtkm::Sqrt(vtkm::MagnitudeSquared(GetTriangleL0<Scalar, Vector, CollectionOfPoints>(pts)));
return l0;
}
/**
* Returns the L1 vector's magnitude, as defined by the verdict manual.
*
* \param [in] pts The three points which define the triangle.
* \return Returns the magnitude of the vector.
*/
template <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetTriangleL1Magnitude(const CollectionOfPoints& pts)
{
const Scalar l1 =
vtkm::Sqrt(vtkm::MagnitudeSquared(GetTriangleL1<Scalar, Vector, CollectionOfPoints>(pts)));
return l1;
}
/**
* Returns the L2 vector's magnitude, as defined by the verdict manual.
*
* \param [in] pts The three points which define the triangle.
* \return Returns the magnitude of the vector.
*/
template <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetTriangleL2Magnitude(const CollectionOfPoints& pts)
{
const Scalar l2 =
vtkm::Sqrt(vtkm::MagnitudeSquared(GetTriangleL2<Scalar, Vector, CollectionOfPoints>(pts)));
return l2;
}
/**
* Returns the Max of the magnitude of each vector which makes up the sides of the triangle.
*
* That is to say, the length of the longest side.
*
* \param [in] pts The three points which define the verticies of the triangle.
*
* \return Returns the max of the triangle side lengths.
*/
template <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetTriangleLMax(const CollectionOfPoints& pts)
{
const Scalar l0 = GetTriangleL0Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l1 = GetTriangleL1Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l2 = GetTriangleL2Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar lmax = vtkm::Max(l0, vtkm::Max(l1, l2));
return lmax;
}
/**
* Returns the Min of the magnitude of each vector which makes up the sides of the triangle.
*
* That is to say, the length of the shortest side.
*
* \param [in] pts The three points which define the verticies of the triangle.
*
* \return Returns the max of the triangle side lengths.
*/
template <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetTriangleLMin(const CollectionOfPoints& pts)
{
const Scalar l0 = GetTriangleL0Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l1 = GetTriangleL1Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l2 = GetTriangleL2Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar lmin = vtkm::Min(l0, vtkm::Min(l1, l2));
return lmin;
}
/**
* Returns the area of the triangle.
*
* \param [in] pts The three points which define the verticies of the triangle.
*
* \return Returns the are of the triangle..
*/
template <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetTriangleArea(const CollectionOfPoints& pts)
{
const Vector L0 = GetTriangleL0<Scalar, Vector, CollectionOfPoints>(pts);
const Vector L1 = GetTriangleL1<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar hhalf(0.5);
const Scalar crossProductMagnitude = vtkm::Sqrt(vtkm::MagnitudeSquared(vtkm::Cross(L0, L1)));
const Scalar area = hhalf * crossProductMagnitude;
return area;
}
/**
* Returns the radius of a circle inscribed within the given triangle. This is commonly denoted as 'r'.
*
* \param [in] pts The three points which define the verticies of the triangle.
*
* \return Returns the inradius.
*/
template <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetTriangleInradius(const CollectionOfPoints& pts)
{
const Scalar two(2.0);
const Scalar area = GetTriangleArea<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l0 = GetTriangleL0Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l1 = GetTriangleL1Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l2 = GetTriangleL2Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar inradius = (two * area) / (l0 + l1 + l2);
return inradius;
}
/**
* Returns the radius of a circle circumscribed around the given triangle. This is commonly denoted as 'R'.
*
* \param [in] pts The three points which define the verticies of the triangle.
*
* \return Returns the circumradius.
*/
template <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetTriangleCircumradius(const CollectionOfPoints& pts)
{
const Scalar four(4.0);
const Scalar area = GetTriangleArea<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l0 = GetTriangleL0Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l1 = GetTriangleL1Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l2 = GetTriangleL2Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar circumradius = (l0 * l1 * l2) / (four * area);
return circumradius;
}
#endif
......@@ -33,26 +33,36 @@ 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?
//This must follow the same order as the MetricNames above
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 +78,8 @@ 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);
void SetOutputName(const std::string& s) { this->OutputName = s; };
template <typename T, typename StorageType, typename DerivedPolicy>
VTKM_CONT vtkm::cont::DataSet DoExecute(
......@@ -80,10 +89,8 @@ public:
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;
CellMetric MyMetric;
std::string OutputName;
};
} // namespace filter
......
......@@ -56,14 +56,16 @@ 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;
this->MyMetric = metric;
if (this->MyMetric < CellMetric::AREA || this->MyMetric >= CellMetric::NUMBER_OF_CELL_METRICS)
{
VTKM_ASSERT(true);
}
this->OutputName = MetricNames[(int)this->MyMetric];
}
template <typename T, typename StorageType, typename DerivedPolicy>
......@@ -75,187 +77,22 @@ 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(this->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";
result.AddField(vtkm::cont::make_FieldCell(s, outArray));
result.AddField(vtkm::cont::make_FieldCell(this->OutputName, outArray));
return result;
}
......
......@@ -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"