From ccedb5aac10394bdfbc84845c3379356d59b068c Mon Sep 17 00:00:00 2001 From: Hank Childs Date: Tue, 17 Sep 2019 13:41:53 -0700 Subject: [PATCH] Jacobian working, unit test working --- vtkm/exec/cellmetrics/CellJacobianMetric.h | 223 ++--- vtkm/exec/cellmetrics/TypeOfCellHexahedral.h | 765 ++++++++++++++++++ .../cellmetrics/TypeOfCellQuadrilateral.h | 468 +++++++++++ vtkm/exec/cellmetrics/TypeOfCellTetrahedral.h | 295 +++++++ vtkm/exec/cellmetrics/TypeOfCellTriangle.h | 192 +++++ .../testing/UnitTestMeshQualityFilter.cxx | 387 +++------ 6 files changed, 1900 insertions(+), 430 deletions(-) create mode 100644 vtkm/exec/cellmetrics/TypeOfCellHexahedral.h create mode 100644 vtkm/exec/cellmetrics/TypeOfCellQuadrilateral.h create mode 100644 vtkm/exec/cellmetrics/TypeOfCellTetrahedral.h create mode 100644 vtkm/exec/cellmetrics/TypeOfCellTriangle.h diff --git a/vtkm/exec/cellmetrics/CellJacobianMetric.h b/vtkm/exec/cellmetrics/CellJacobianMetric.h index 1e9f2e1c8..51f54beb1 100644 --- a/vtkm/exec/cellmetrics/CellJacobianMetric.h +++ b/vtkm/exec/cellmetrics/CellJacobianMetric.h @@ -1,3 +1,4 @@ + //============================================================================ // Copyright (c) Kitware, Inc. // All rights reserved. @@ -21,17 +22,20 @@ #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 + * 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. + * 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" @@ -47,63 +51,9 @@ namespace exec namespace cellmetrics { -using FloatType = vtkm::FloatDefault; - -template -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(); - - 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()); //normal case - - return vtkm::Max(minCalculation, OutType(-1) * vtkm::Infinity()); -} - -template -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()); //normal case - - return vtkm::Max(minDeterminant, vtkm::NegativeInfinity()); -} - - // ========================= Unsupported cells ================================== -// By default, cells have zero shape unless the shape type template is specialized below. +// By default, cells return the metric 0.0 unless the shape type template is specialized below. template VTKM_EXEC OutType CellJacobianMetric(const vtkm::IdComponent& numPts, const PointCoordVecType& pts, @@ -116,63 +66,9 @@ VTKM_EXEC OutType CellJacobianMetric(const vtkm::IdComponent& numPts, return OutType(0.0); } -template -VTKM_EXEC OutType CellJacobianMetric(const vtkm::IdComponent& numPts, - const PointCoordVecType& pts, - vtkm::CellShapeTagPolygon, - const vtkm::exec::FunctorBase& worklet) -{ - switch (numPts) - { - case 4: - return CellJacobianMetric(numPts, pts, vtkm::CellShapeTagQuad(), worklet); - default: - break; - } - return OutType(-1.0); -} - -template -VTKM_EXEC OutType CellJacobianMetric(const vtkm::IdComponent&, - const PointCoordVecType&, - vtkm::CellShapeTagLine, - const vtkm::exec::FunctorBase& worklet) -{ - UNUSED(worklet); - return OutType(-1.0); -} - -template -VTKM_EXEC OutType CellJacobianMetric(const vtkm::IdComponent&, - const PointCoordVecType&, - vtkm::CellShapeTagTriangle, - const vtkm::exec::FunctorBase& worklet) -{ - UNUSED(worklet); - return OutType(-1.0); -} - -template -VTKM_EXEC OutType CellJacobianMetric(const vtkm::IdComponent&, - const PointCoordVecType&, - vtkm::CellShapeTagWedge, - const vtkm::exec::FunctorBase& worklet) -{ - UNUSED(worklet); - return OutType(-1.0); -} - -template -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. + +// Compute the Jacobian of a quadrilateral. // Formula: min{Jacobian at each vertex} // Equals 1 for a unit square // Acceptable range: [0,FLOAT_MAX] @@ -190,24 +86,25 @@ VTKM_EXEC OutType CellJacobianMetric(const vtkm::IdComponent& numPts, return OutType(0.0); } + using Scalar = OutType; + using CollectionOfPoints = PointCoordVecType; + using Vector = typename PointCoordVecType::ComponentType; - //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( - vtkm::make_VecC(QuadEdgesToUse, 4), vtkm::make_VecC(QuadAxes, 2)); + const Scalar alpha0 = GetQuadAlpha0(pts); + const Scalar alpha1 = GetQuadAlpha1(pts); + const Scalar alpha2 = GetQuadAlpha2(pts); + const Scalar alpha3 = GetQuadAlpha3(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}, Alpha_8*/64} -// -Alpha_i -> jacobian determinant at respective vertex -// -Alpha_8 -> jacobian at center +// 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] @@ -224,31 +121,37 @@ VTKM_EXEC OutType CellJacobianMetric(const vtkm::IdComponent& numPts, return OutType(0.0); } + using Scalar = OutType; + using CollectionOfPoints = PointCoordVecType; + using Vector = typename PointCoordVecType::ComponentType; - //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 Scalar alpha0 = GetHexAlphai(pts, vtkm::Id(0)); + const Scalar alpha1 = GetHexAlphai(pts, vtkm::Id(1)); + const Scalar alpha2 = GetHexAlphai(pts, vtkm::Id(2)); + const Scalar alpha3 = GetHexAlphai(pts, vtkm::Id(3)); + const Scalar alpha4 = GetHexAlphai(pts, vtkm::Id(4)); + const Scalar alpha5 = GetHexAlphai(pts, vtkm::Id(5)); + const Scalar alpha6 = GetHexAlphai(pts, vtkm::Id(6)); + const Scalar alpha7 = GetHexAlphai(pts, vtkm::Id(7)); + const Scalar alpha8 = GetHexAlphai(pts, vtkm::Id(8)); + const Scalar alpha8Div64 = alpha8 / Scalar(64.0); - 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( - vtkm::make_VecC(hexMatrices, 12)); + 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 * L0) * L3 +// 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] @@ -261,18 +164,24 @@ VTKM_EXEC OutType CellJacobianMetric(const vtkm::IdComponent& numPts, { if (numPts != 4) { - worklet.RaiseError("Jacobian metric requires 4 points"); + worklet.RaiseError("Jacobian metric (tetra) 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]); -} + using Scalar = OutType; + using CollectionOfPoints = PointCoordVecType; + using Vector = typename PointCoordVecType::ComponentType; + const Vector L0 = GetTetraL0(pts); + const Vector L2 = GetTetraL2(pts); + const Vector L3 = GetTetraL3(pts); + + const Scalar q = vtkm::Dot(vtkm::Cross(L2, L0), L3); + + return q; +} } // namespace cellmetrics } // namespace exec } // namespace vtkm -#endif // vtk_m_exec_cellmetrics_CellEdgeRatioMetric_h +#endif // vtk_m_exec_cellmetrics_CellJacobianMetric_h diff --git a/vtkm/exec/cellmetrics/TypeOfCellHexahedral.h b/vtkm/exec/cellmetrics/TypeOfCellHexahedral.h new file mode 100644 index 000000000..428c4874b --- /dev/null +++ b/vtkm/exec/cellmetrics/TypeOfCellHexahedral.h @@ -0,0 +1,765 @@ +#ifndef vtk_m_exec_cellmetrics_TypeOfCellHexahedral +#define vtk_m_exec_cellmetrics_TypeOfCellHexahedral +/** + * The Verdict manual defines a set of commonly + * used components of a hexahedra (hex). For example, + * area, edge lengths, and so forth. + * + * These definitions can be found starting on + * page 77 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 eight points which define the hexahedra. + * \return Returns the vector. + */ +template +VTKM_EXEC Vector GetHexL0(const CollectionOfPoints& pts) +{ + const Vector L0(pts[1] - pts[0]); + return L0; +} + +/** + * Returns the L1 vector, as defined by the verdict manual. + * + * \param [in] pts The eight points which define the hexahedra. + * \return Returns the vector. + */ +template +VTKM_EXEC Vector GetHexL1(const CollectionOfPoints& pts) +{ + const Vector L1(pts[2] - pts[1]); + return L1; +} + +/** + * Returns the L2 vector, as defined by the verdict manual. + * + * \param [in] pts The eight points which define the hexahedra. + * \return Returns the vector. + */ +template +VTKM_EXEC Vector GetHexL2(const CollectionOfPoints& pts) +{ + const Vector L2(pts[3] - pts[2]); + return L2; +} + +/** + * Returns the L3 vector, as defined by the verdict manual. + * + * \param [in] pts The eight points which define the hexahedra. + * \return Returns the vector. + */ +template +VTKM_EXEC Vector GetHexL3(const CollectionOfPoints& pts) +{ + const Vector L3(pts[3] - pts[0]); + return L3; +} + +/** + * Returns the L4 vector, as defined by the verdict manual. + * + * \param [in] pts The eight points which define the hexahedra. + * \return Returns the vector. + */ +template +VTKM_EXEC Vector GetHexL4(const CollectionOfPoints& pts) +{ + const Vector L4(pts[4] - pts[0]); + return L4; +} + +/** + * Returns the L5 vector, as defined by the verdict manual. + * + * \param [in] pts The eight points which define the hexahedra. + * \return Returns the vector. + */ +template +VTKM_EXEC Vector GetHexL5(const CollectionOfPoints& pts) +{ + const Vector L5(pts[5] - pts[1]); + return L5; +} + +/** + * Returns the L6 vector, as defined by the verdict manual. + * + * \param [in] pts The eight points which define the hexahedra. + * \return Returns the vector. + */ +template +VTKM_EXEC Vector GetHexL6(const CollectionOfPoints& pts) +{ + const Vector L6(pts[6] - pts[2]); + return L6; +} + +/** + * Returns the L7 vector, as defined by the verdict manual. + * + * \param [in] pts The eight points which define the hexahedra. + * \return Returns the vector. + */ +template +VTKM_EXEC Vector GetHexL7(const CollectionOfPoints& pts) +{ + const Vector L7(pts[7] - pts[3]); + return L7; +} + +/** + * Returns the L8 vector, as defined by the verdict manual. + * + * \param [in] pts The eight points which define the hexahedra. + * \return Returns the vector. + */ +template +VTKM_EXEC Vector GetHexL8(const CollectionOfPoints& pts) +{ + const Vector L8(pts[5] - pts[4]); + return L8; +} + +/** + * Returns the L9 vector, as defined by the verdict manual. + * + * \param [in] pts The eight points which define the hexahedra. + * \return Returns the vector. + */ +template +VTKM_EXEC Vector GetHexL9(const CollectionOfPoints& pts) +{ + const Vector L9(pts[6] - pts[5]); + return L9; +} + +/** + * Returns the L10 vector, as defined by the verdict manual. + * + * \param [in] pts The eight points which define the hexahedra. + * \return Returns the vector. + */ +template +VTKM_EXEC Vector GetHexL10(const CollectionOfPoints& pts) +{ + const Vector L10(pts[7] - pts[6]); + return L10; +} +/** + * Returns the L11 vector, as defined by the verdict manual. + * + * \param [in] pts The eight points which define the hexahedra. + * \return Returns the vector. + */ +template +VTKM_EXEC Vector GetHexL11(const CollectionOfPoints& pts) +{ + const Vector L11(pts[7] - pts[4]); + return L11; +} +/** + * Returns the L0 vector's magnitude, as defined by the verdict manual. + * + * \param [in] pts The eight points which define the Hex. + * \return Returns the magnitude of the vector. + */ +template +VTKM_EXEC Scalar GetHexL0Magnitude(const CollectionOfPoints& pts) +{ + const Scalar l0 = + vtkm::Sqrt(vtkm::MagnitudeSquared(GetHexL0(pts))); + return l0; +} +/** + * Returns the L1 vector's magnitude, as defined by the verdict manual. + * + * \param [in] pts The eight points which define the Hex. + * \return Returns the magnitude of the vector. + */ +template +VTKM_EXEC Scalar GetHexL1Magnitude(const CollectionOfPoints& pts) +{ + const Scalar l1 = + vtkm::Sqrt(vtkm::MagnitudeSquared(GetHexL1(pts))); + return l1; +} +/** + * Returns the L2 vector's magnitude, as defined by the verdict manual. + * + * \param [in] pts The eight points which define the Hex. + * \return Returns the magnitude of the vector. + */ +template +VTKM_EXEC Scalar GetHexL2Magnitude(const CollectionOfPoints& pts) +{ + const Scalar l2 = + vtkm::Sqrt(vtkm::MagnitudeSquared(GetHexL2(pts))); + return l2; +} +/** + * Returns the L3 vector's magnitude, as defined by the verdict manual. + * + * \param [in] pts The eight points which define the Hex. + * \return Returns the magnitude of the vector. + */ +template +VTKM_EXEC Scalar GetHexL3Magnitude(const CollectionOfPoints& pts) +{ + const Scalar l3 = + vtkm::Sqrt(vtkm::MagnitudeSquared(GetHexL3(pts))); + return l3; +} +/** + * Returns the L4 vector's magnitude, as defined by the verdict manual. + * + * \param [in] pts The eight points which define the Hex. + * \return Returns the magnitude of the vector. + */ +template +VTKM_EXEC Scalar GetHexL4Magnitude(const CollectionOfPoints& pts) +{ + const Scalar l4 = + vtkm::Sqrt(vtkm::MagnitudeSquared(GetHexL4(pts))); + return l4; +} +/** + * Returns the L5 vector's magnitude, as defined by the verdict manual. + * + * \param [in] pts The eight points which define the Hex. + * \return Returns the magnitude of the vector. + */ +template +VTKM_EXEC Scalar GetHexL5Magnitude(const CollectionOfPoints& pts) +{ + const Scalar l5 = + vtkm::Sqrt(vtkm::MagnitudeSquared(GetHexL5(pts))); + return l5; +} +/** + * Returns the L6 vector's magnitude, as defined by the verdict manual. + * + * \param [in] pts The eight points which define the Hex. + * \return Returns the magnitude of the vector. + */ +template +VTKM_EXEC Scalar GetHexL6Magnitude(const CollectionOfPoints& pts) +{ + const Scalar l6 = + vtkm::Sqrt(vtkm::MagnitudeSquared(GetHexL6(pts))); + return l6; +} +/** + * Returns the L7 vector's magnitude, as defined by the verdict manual. + * + * \param [in] pts The eight points which define the Hex. + * \return Returns the magnitude of the vector. + */ +template +VTKM_EXEC Scalar GetHexL7Magnitude(const CollectionOfPoints& pts) +{ + const Scalar l7 = + vtkm::Sqrt(vtkm::MagnitudeSquared(GetHexL7(pts))); + return l7; +} +/** + * Returns the L8 vector's magnitude, as defined by the verdict manual. + * + * \param [in] pts The eight points which define the Hex. + * \return Returns the magnitude of the vector. + */ +template +VTKM_EXEC Scalar GetHexL8Magnitude(const CollectionOfPoints& pts) +{ + const Scalar l8 = + vtkm::Sqrt(vtkm::MagnitudeSquared(GetHexL8(pts))); + return l8; +} +/** + * Returns the L9 vector's magnitude, as defined by the verdict manual. + * + * \param [in] pts The eight points which define the Hex. + * \return Returns the magnitude of the vector. + */ +template +VTKM_EXEC Scalar GetHexL9Magnitude(const CollectionOfPoints& pts) +{ + const Scalar l9 = + vtkm::Sqrt(vtkm::MagnitudeSquared(GetHexL9(pts))); + return l9; +} +/** + * Returns the L10 vector's magnitude, as defined by the verdict manual. + * + * \param [in] pts The eight points which define the Hex. + * \return Returns the magnitude of the vector. + */ +template +VTKM_EXEC Scalar GetHexL10Magnitude(const CollectionOfPoints& pts) +{ + const Scalar l10 = + vtkm::Sqrt(vtkm::MagnitudeSquared(GetHexL10(pts))); + return l10; +} +/** + * Returns the L11 vector's magnitude, as defined by the verdict manual. + * + * \param [in] pts The eight points which define the Hex. + * \return Returns the magnitude of the vector. + */ +template +VTKM_EXEC Scalar GetHexL11Magnitude(const CollectionOfPoints& pts) +{ + const Scalar l11 = + vtkm::Sqrt(vtkm::MagnitudeSquared(GetHexL11(pts))); + return l11; +} +/** + * Returns the Max of the magnitude of each vector which makes up the sides of the Hex. + * + * \param [in] pts The eight points which define the Hex. + * \return Returns the magnitude of the vector. + */ +template +VTKM_EXEC Scalar GetHexLMax(const CollectionOfPoints& pts) +{ + const Scalar l0 = GetHexL0Magnitude(pts); + const Scalar l1 = GetHexL1Magnitude(pts); + const Scalar l2 = GetHexL2Magnitude(pts); + const Scalar l3 = GetHexL3Magnitude(pts); + const Scalar l4 = GetHexL4Magnitude(pts); + const Scalar l5 = GetHexL5Magnitude(pts); + const Scalar l6 = GetHexL6Magnitude(pts); + const Scalar l7 = GetHexL7Magnitude(pts); + const Scalar l8 = GetHexL8Magnitude(pts); + const Scalar l9 = GetHexL9Magnitude(pts); + const Scalar l10 = GetHexL10Magnitude(pts); + const Scalar l11 = GetHexL11Magnitude(pts); + + const Scalar lmax = vtkm::Max( + l0, + vtkm::Max( + l1, + vtkm::Max( + l2, + vtkm::Max( + l3, + vtkm::Max( + l4, + vtkm::Max( + l5, + vtkm::Max(l6, vtkm::Max(l7, vtkm::Max(l8, vtkm::Max(l9, vtkm::Max(l10, l11))))))))))); + + return lmax; +} + +/** + * Returns the Min of the magnitude of each vector which makes up the sides of the Hex. + * + * \param [in] pts The eight points which define the Hex. + * \return Returns the magnitude of the vector. + */ +template +VTKM_EXEC Scalar GetHexLMin(const CollectionOfPoints& pts) +{ + const Scalar l0 = GetHexL0Magnitude(pts); + const Scalar l1 = GetHexL1Magnitude(pts); + const Scalar l2 = GetHexL2Magnitude(pts); + const Scalar l3 = GetHexL3Magnitude(pts); + const Scalar l4 = GetHexL4Magnitude(pts); + const Scalar l5 = GetHexL5Magnitude(pts); + const Scalar l6 = GetHexL6Magnitude(pts); + const Scalar l7 = GetHexL7Magnitude(pts); + const Scalar l8 = GetHexL8Magnitude(pts); + const Scalar l9 = GetHexL9Magnitude(pts); + const Scalar l10 = GetHexL10Magnitude(pts); + const Scalar l11 = GetHexL11Magnitude(pts); + + const Scalar lmin = vtkm::Min( + l0, + vtkm::Min( + l1, + vtkm::Min( + l2, + vtkm::Min( + l3, + vtkm::Min( + l4, + vtkm::Min( + l5, + vtkm::Min(l6, vtkm::Min(l7, vtkm::Min(l8, vtkm::Min(l9, vtkm::Min(l10, l11))))))))))); + + return lmin; +} + +/** + * Returns the D0 vector, as defined by the verdict manual. + * + * \param [in] pts The eight points which define the hexahedra. + * \return Returns the vector. + */ +template +VTKM_EXEC Vector GetHexD0(const CollectionOfPoints& pts) +{ + const Vector D0((pts[6] - pts[0])); + return D0; +} + +/** + * Returns the D1 vector, as defined by the verdict manual. + * + * \param [in] pts The eight points which define the hexahedra. + * \return Returns the vector. + */ +template +VTKM_EXEC Vector GetHexD1(const CollectionOfPoints& pts) +{ + const Vector D1(pts[7] - pts[1]); + return D1; +} + +/** + * Returns the D2 vector, as defined by the verdict manual. + * + * \param [in] pts The eight points which define the hexahedra. + * \return Returns the vector. + */ +template +VTKM_EXEC Vector GetHexD2(const CollectionOfPoints& pts) +{ + const Vector D2(pts[4] - pts[2]); + return D2; +} + +/** + * Returns the D3 vector, as defined by the verdict manual. + * + * \param [in] pts The eight points which define the hexahedra. + * \return Returns the vector. + */ +template +VTKM_EXEC Vector GetHexD3(const CollectionOfPoints& pts) +{ + const Vector D3(pts[5] - pts[3]); + return D3; +} + +/** + * Returns the D0 vector's magnitude, as defined by the verdict manual. + * + * \param [in] pts The eight points which define the Hex. + * \return Returns the magnitude of the vector. + */ +template +VTKM_EXEC Scalar GetHexD0Magnitude(const CollectionOfPoints& pts) +{ + const Scalar d0 = + vtkm::Sqrt(vtkm::MagnitudeSquared(GetHexD0(pts))); + return d0; +} + +/** + * Returns the D1 vector's magnitude, as defined by the verdict manual. + * + * \param [in] pts The eight points which define the Hex. + * \return Returns the magnitude of the vector. + */ +template +VTKM_EXEC Scalar GetHexD1Magnitude(const CollectionOfPoints& pts) +{ + const Scalar d1 = + vtkm::Sqrt(vtkm::MagnitudeSquared(GetHexD1(pts))); + return d1; +} + +/** + * Returns the D2 vector's magnitude, as defined by the verdict manual. + * + * \param [in] pts The eight points which define the Hex. + * \return Returns the magnitude of the vector. + */ +template +VTKM_EXEC Scalar GetHexD2Magnitude(const CollectionOfPoints& pts) +{ + const Scalar d2 = + vtkm::Sqrt(vtkm::MagnitudeSquared(GetHexD2(pts))); + return d2; +} + +/** + * Returns the D3 vector's magnitude, as defined by the verdict manual. + * + * \param [in] pts The eight points which define the Hex. + * \return Returns the magnitude of the vector. + */ +template +VTKM_EXEC Scalar GetHexD3Magnitude(const CollectionOfPoints& pts) +{ + const Scalar d3 = + vtkm::Sqrt(vtkm::MagnitudeSquared(GetHexD3(pts))); + return d3; +} + +/** + * Returns the Min of the magnitude of each vector which makes up the sides of the Hex. + * + * \param [in] pts The eight points which define the Hex. + * \return Returns the magnitude of the vector. + */ +template +VTKM_EXEC Scalar GetHexDMin(const CollectionOfPoints& pts) +{ + const Scalar d0 = GetHexD0Magnitude(pts); + const Scalar d1 = GetHexD1Magnitude(pts); + const Scalar d2 = GetHexD2Magnitude(pts); + const Scalar d3 = GetHexD3Magnitude(pts); + + const Scalar dmin = vtkm::Min(d0, vtkm::Min(d1, vtkm::Min(d2, d3))); + + return dmin; +} + +/** + * Returns the Min of the magnitude of each vector which makes up the sides of the Hex. + * + * \param [in] pts The eight points which define the Hex. + * \return Returns the magnitude of the vector. + */ +template +VTKM_EXEC Scalar GetHexDMax(const CollectionOfPoints& pts) +{ + const Scalar d0 = GetHexD0Magnitude(pts); + const Scalar d1 = GetHexD1Magnitude(pts); + const Scalar d2 = GetHexD2Magnitude(pts); + const Scalar d3 = GetHexD3Magnitude(pts); + + const Scalar dmax = vtkm::Max(d0, vtkm::Max(d1, vtkm::Max(d2, d3))); + + return dmax; +} + +/** + * Returns the X1 vector defined in the verdict manual. + * + * \param [in] pts The eight points which define the Hex. + * \return Returns the vector. + */ +template +VTKM_EXEC Vector GetHexX1(const CollectionOfPoints& pts) +{ + const Vector X1((pts[1] - pts[0]) + (pts[2] - pts[3]) + (pts[5] - pts[4]) + (pts[6] - pts[7])); + return X1; +} + +/** + * Returns the X2 vector defined in the verdict manual. + * + * \param [in] pts The eight points which define the Hex. + * \return Returns the vector. + */ +template +VTKM_EXEC Vector GetHexX2(const CollectionOfPoints& pts) +{ + const Vector X2((pts[3] - pts[0]) + (pts[2] - pts[1]) + (pts[7] - pts[4]) + (pts[6] - pts[5])); + return X2; +} + +/** + * Returns the X3 vector defined in the verdict manual. + * + * \param [in] pts The eight points which define the Hex. + * \return Returns the vector. + */ +template +VTKM_EXEC Vector GetHexX3(const CollectionOfPoints& pts) +{ + const Vector X3((pts[4] - pts[0]) + (pts[5] - pts[1]) + (pts[6] - pts[2]) + (pts[7] - pts[3])); + return X3; +} + +/** + * Returns the A_i matrix defined in the verdict manual. + * + * \param [in] pts The eight points which define the Hex. + * \param [in] index The index of A to load. + * \return Returns the vector. + */ +template +VTKM_EXEC vtkm::Vec GetHexAi(const CollectionOfPoints& pts, const vtkm::Id& index) +{ + const Scalar neg1(-1.0); + if (index == 0) + { + const Vector v0 = GetHexL0(pts); + const Vector v1 = GetHexL3(pts); + const Vector v2 = GetHexL4(pts); + const vtkm::Vec A = { v0, v1, v2 }; + return A; + } + else if (index == 1) + { + const Vector v0 = GetHexL1(pts); + const Vector v1 = neg1 * GetHexL0(pts); + const Vector v2 = GetHexL5(pts); + const vtkm::Vec A = { v0, v1, v2 }; + return A; + } + else if (index == 2) + { + const Vector v0 = GetHexL2(pts); + const Vector v1 = neg1 * GetHexL1(pts); + const Vector v2 = GetHexL6(pts); + const vtkm::Vec A = { v0, v1, v2 }; + return A; + } + else if (index == 3) + { + const Vector v0 = neg1 * GetHexL3(pts); + const Vector v1 = neg1 * GetHexL2(pts); + const Vector v2 = GetHexL7(pts); + const vtkm::Vec A = { v0, v1, v2 }; + return A; + } + else if (index == 4) + { + const Vector v0 = GetHexL11(pts); + const Vector v1 = GetHexL8(pts); + const Vector v2 = neg1 * GetHexL4(pts); + const vtkm::Vec A = { v0, v1, v2 }; + return A; + } + else if (index == 5) + { + const Vector v0 = neg1 * GetHexL8(pts); + const Vector v1 = GetHexL9(pts); + const Vector v2 = neg1 * GetHexL5(pts); + const vtkm::Vec A = { v0, v1, v2 }; + return A; + } + else if (index == 6) + { + const Vector v0 = neg1 * GetHexL9(pts); + const Vector v1 = GetHexL10(pts); + const Vector v2 = neg1 * GetHexL6(pts); + const vtkm::Vec A = { v0, v1, v2 }; + return A; + } + else if (index == 7) + { + const Vector v0 = neg1 * GetHexL10(pts); + const Vector v1 = neg1 * GetHexL11(pts); + const Vector v2 = neg1 * GetHexL7(pts); + const vtkm::Vec A = { v0, v1, v2 }; + return A; + } + else + { + const Vector v0 = GetHexX1(pts); + const Vector v1 = GetHexX2(pts); + const Vector v2 = GetHexX3(pts); + const vtkm::Vec A = { v0, v1, v2 }; + return A; + } +} + +/** + * Returns ||A_i||^2 as defined in the verdict manual. + * + * \param [in] pts The eight points which define the Hex. + * \param [in] index The index of A to load. + * \return Returns the vector. + */ +template +VTKM_EXEC Scalar GetHexAiNormSquared(const CollectionOfPoints& pts, const vtkm::Id& index) +{ + const vtkm::Vec Ai = GetHexAi(pts, index); + const Scalar magSquared0 = vtkm::MagnitudeSquared(Ai[0]); + const Scalar magSquared1 = vtkm::MagnitudeSquared(Ai[1]); + const Scalar magSquared2 = vtkm::MagnitudeSquared(Ai[2]); + + const Scalar AiNormSquared = magSquared0 + magSquared1 + magSquared2; + + return AiNormSquared; +} + +/** + * Returns ||adj(A_i)||^2 as defined in the verdict manual. + * + * \param [in] pts The eight points which define the Hex. + * \param [in] index The index of A to load. + * \return Returns the vector. + */ +template +VTKM_EXEC Scalar GetHexAiAdjNormSquared(const CollectionOfPoints& pts, const vtkm::Id& index) +{ + const vtkm::Vec Ai = GetHexAi(pts, index); + const Scalar magSquared0 = vtkm::MagnitudeSquared(vtkm::Cross(Ai[0], Ai[1])); + const Scalar magSquared1 = vtkm::MagnitudeSquared(vtkm::Cross(Ai[1], Ai[2])); + const Scalar magSquared2 = vtkm::MagnitudeSquared(vtkm::Cross(Ai[2], Ai[0])); + + const Scalar AiAdjNormSquared = magSquared0 + magSquared1 + magSquared2; + + return AiAdjNormSquared; +} + +/** + * Returns alpha_i, the determinant of A_i, as defined in the verdict manual. + * + * \param [in] pts The eight points which define the Hex. + * \param [in] index The index of A to load. + * \return Returns the vector. + */ +template +VTKM_EXEC Scalar GetHexAlphai(const CollectionOfPoints& pts, const vtkm::Id& index) +{ + const vtkm::Vec Ai = GetHexAi(pts, index); + const Scalar alpha_i = vtkm::Dot(Ai[0], vtkm::Cross(Ai[1], Ai[2])); + + return alpha_i; +} + +/** + * Returns hat{A}_i, the "normalized" version of A_i, as defined in the verdict manual. + * + * \param [in] pts The eight points which define the Hex. + * \param [in] index The index of A to load. + * \return Returns the vector. + */ +template +VTKM_EXEC vtkm::Vec GetHexAiHat(const CollectionOfPoints& pts, const vtkm::Id& index) +{ + const vtkm::Vec Ai = GetHexAi(pts, index); + const Vector v0hat = Ai[0] / vtkm::Sqrt(vtkm::MagnitudeSquared(Ai[0])); + const Vector v1hat = Ai[1] / vtkm::Sqrt(vtkm::MagnitudeSquared(Ai[1])); + const Vector v2hat = Ai[2] / vtkm::Sqrt(vtkm::MagnitudeSquared(Ai[2])); + + const vtkm::Vec Ahat = { v0hat, v1hat, v2hat }; + return Ahat; + + return Ahat; +} +/** + * Returns hat{alpha}_i, the determinant of hat{A}_i, as defined in the verdict manual. + * + * \param [in] pts The eight points which define the Hex. + * \param [in] index The index of A to load. + * \return Returns the vector. + */ +template +VTKM_EXEC Scalar GetHexAlphaiHat(const CollectionOfPoints& pts, const vtkm::Id& index) +{ + const vtkm::Vec Ai = GetHexAiHat(pts, index); + const Scalar hatAlpha_i = vtkm::Dot(Ai[0], vtkm::Cross(Ai[1], Ai[2])); + + return hatAlpha_i; +} +#endif diff --git a/vtkm/exec/cellmetrics/TypeOfCellQuadrilateral.h b/vtkm/exec/cellmetrics/TypeOfCellQuadrilateral.h new file mode 100644 index 000000000..1091615ab --- /dev/null +++ b/vtkm/exec/cellmetrics/TypeOfCellQuadrilateral.h @@ -0,0 +1,468 @@ +#ifndef vtk_m_exec_cellmetrics_TypeOfCellQuadrilateral +#define vtk_m_exec_cellmetrics_TypeOfCellQuadrilateral +/** + * The Verdict manual defines a set of commonly + * used components of a quadrilateral (quad). For example, + * area, side lengths, and so forth. + * + * These definitions can be found starting on + * page 32 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 four points which define the quadrilateral. + * \return Returns the vector. + */ +template +VTKM_EXEC Vector GetQuadL0(const CollectionOfPoints& pts) +{ + const Vector L0(pts[1] - pts[0]); + return L0; +} + +/** + * Returns the L1 vector, as defined by the verdict manual. + * + * \param [in] pts The four points which define the quadrilateral. + * \return Returns the vector. + */ +template +VTKM_EXEC Vector GetQuadL1(const CollectionOfPoints& pts) +{ + const Vector L1(pts[2] - pts[1]); + return L1; +} + +/** + * Returns the L2 vector, as defined by the verdict manual. + * + * \param [in] pts The four points which define the quadrilateral. + * \return Returns the vector. + */ +template +VTKM_EXEC Vector GetQuadL2(const CollectionOfPoints& pts) +{ + const Vector L2(pts[3] - pts[2]); + return L2; +} + +/** + * Returns the L3 vector, as defined by the verdict manual. + * + * \param [in] pts The four points which define the quadrilateral. + * \return Returns the vector. + */ +template +VTKM_EXEC Vector GetQuadL3(const CollectionOfPoints& pts) +{ + const Vector L3(pts[0] - pts[3]); + return L3; +} + +/** + * Returns the L0 vector's magnitude, as defined by the verdict manual. + * + * \param [in] pts The four points which define the Quad. + * \return Returns the magnitude of the vector. + */ +template +VTKM_EXEC Scalar GetQuadL0Magnitude(const CollectionOfPoints& pts) +{ + const Scalar l0 = + vtkm::Sqrt(vtkm::MagnitudeSquared(GetQuadL0(pts))); + return l0; +} + +/** + * Returns the L1 vector's magnitude, as defined by the verdict manual. + * + * \param [in] pts The four points which define the Quad. + * \return Returns the magnitude of the vector. + */ +template +VTKM_EXEC Scalar GetQuadL1Magnitude(const CollectionOfPoints& pts) +{ + const Scalar l1 = + vtkm::Sqrt(vtkm::MagnitudeSquared(GetQuadL1(pts))); + return l1; +} + +/** + * Returns the L2 vector's magnitude, as defined by the verdict manual. + * + * \param [in] pts The four points which define the Quad. + * \return Returns the magnitude of the vector. + */ +template +VTKM_EXEC Scalar GetQuadL2Magnitude(const CollectionOfPoints& pts) +{ + const Scalar l2 = + vtkm::Sqrt(vtkm::MagnitudeSquared(GetQuadL2(pts))); + return l2; +} + +/** + * Returns the L3 vector's magnitude, as defined by the verdict manual. + * + * \param [in] pts The four points which define the Quad. + * \return Returns the magnitude of the vector. + */ +template +VTKM_EXEC Scalar GetQuadL3Magnitude(const CollectionOfPoints& pts) +{ + const Scalar l3 = + vtkm::Sqrt(vtkm::MagnitudeSquared(GetQuadL3(pts))); + return l3; +} + +/** + * Returns the Max of the magnitude of each vector which makes up the sides of the Quad. + * + * \param [in] pts The four points which define the Quad. + * \return Returns the magnitude of the vector. + */ +template +VTKM_EXEC Scalar GetQuadLMax(const CollectionOfPoints& pts) +{ + const Scalar l0 = GetQuadL0Magnitude(pts); + const Scalar l1 = GetQuadL1Magnitude(pts); + const Scalar l2 = GetQuadL2Magnitude(pts); + const Scalar l3 = GetQuadL3Magnitude(pts); + const Scalar lmax = vtkm::Max(l0, vtkm::Max(l1, vtkm::Max(l2, l3))); + return lmax; +} + +/** + * Returns the Min of the magnitude of each vector which makes up the sides of the Quad. + * + * \param [in] pts The four points which define the Quad. + * \return Returns the magnitude of the vector. + */ +template +VTKM_EXEC Scalar GetQuadLMin(const CollectionOfPoints& pts) +{ + const Scalar l0 = GetQuadL0Magnitude(pts); + const Scalar l1 = GetQuadL1Magnitude(pts); + const Scalar l2 = GetQuadL2Magnitude(pts); + const Scalar l3 = GetQuadL3Magnitude(pts); + const Scalar lmin = vtkm::Min(l0, vtkm::Min(l1, vtkm::Min(l2, l3))); + return lmin; +} + +/** + * Returns the D0 vector, as defined by the verdict manual. + * + * \param [in] pts The four points which define the quadrilateral. + * \return Returns the vector. + */ +template +VTKM_EXEC Vector GetQuadD0(const CollectionOfPoints& pts) +{ + const Vector D0(pts[2] - pts[0]); + return D0; +} + +/** + * Returns the D1 vector, as defined by the verdict manual. + * + * \param [in] pts The four points which define the quadrilateral. + * \return Returns the vector. + */ +template +VTKM_EXEC Vector GetQuadD1(const CollectionOfPoints& pts) +{ + const Vector D1(pts[3] - pts[1]); + return D1; +} + +/** + * Returns the D0 vector's magnitude, as defined by the verdict manual. + * + * \param [in] pts The four points which define the Quad. + * \return Returns the magnitude of the vector. + */ +template +VTKM_EXEC Scalar GetQuadD0Magnitude(const CollectionOfPoints& pts) +{ + const Scalar d0 = + vtkm::Sqrt(vtkm::MagnitudeSquared(GetQuadD0(pts))); + return d0; +} + +/** + * Returns the D0 vector's magnitude, as defined by the verdict manual. + * + * \param [in] pts The four points which define the Quad. + * \return Returns the magnitude of the vector. + */ +template +VTKM_EXEC Scalar GetQuadD1Magnitude(const CollectionOfPoints& pts) +{ + const Scalar d1 = + vtkm::Sqrt(vtkm::MagnitudeSquared(GetQuadD1(pts))); + return d1; +} + +/** + * Returns the Max of the magnitude of each vector which makes up the diagonals of the Quad. + * + * \param [in] pts The four points which define the Quad. + * \return Returns the magnitude of the vector. + */ +template +VTKM_EXEC Scalar GetQuadDMax(const CollectionOfPoints& pts) +{ + const Scalar d0 = GetQuadD0Magnitude(pts); + const Scalar d1 = GetQuadD1Magnitude(pts); + + const Scalar dmax = vtkm::Max(d0, d1); + return dmax; +} + +/** + * Returns the X0 vector, as defined by the verdict manual. + * + * \param [in] pts The four points which define the quadrilateral. + * \return Returns the vector. + */ +template +VTKM_EXEC Vector GetQuadX0(const CollectionOfPoints& pts) +{ + const Vector X0((pts[1] - pts[0]) + (pts[2] - pts[3])); + return X0; +} + +/** + * Returns the X1 vector, as defined by the verdict manual. + * + * \param [in] pts The four points which define the quadrilateral. + * \return Returns the vector. + */ +template +VTKM_EXEC Vector GetQuadX1(const CollectionOfPoints& pts) +{ + const Vector X1((pts[2] - pts[1]) + (pts[3] - pts[0])); + return X1; +} + +/** + * Returns the N0 vector, as defined by the verdict manual. + * + * \param [in] pts The four points which define the quadrilateral. + * \return Returns the vector. + */ +template +VTKM_EXEC Vector GetQuadN0(const CollectionOfPoints& pts) +{ + const Vector A = GetQuadL3(pts); + const Vector B = GetQuadL0(pts); + const Vector N0 = vtkm::Cross(A, B); + return N0; +} + +/** + * Returns the N1 vector, as defined by the verdict manual. + * + * \param [in] pts The four points which define the quadrilateral. + * \return Returns the vector. + */ +template +VTKM_EXEC Vector GetQuadN1(const CollectionOfPoints& pts) +{ + const Vector A = GetQuadL0(pts); + const Vector B = GetQuadL1(pts); + const Vector N1 = vtkm::Cross(A, B); + return N1; +} + +/** + * Returns the N2 vector, as defined by the verdict manual. + * + * \param [in] pts The four points which define the quadrilateral. + * \return Returns the vector. + */ +template +VTKM_EXEC Vector GetQuadN2(const CollectionOfPoints& pts) +{ + const Vector A = GetQuadL1(pts); + const Vector B = GetQuadL2(pts); + const Vector N2 = vtkm::Cross(A, B); + return N2; +} + +/** + * Returns the N3 vector, as defined by the verdict manual. + * + * \param [in] pts The four points which define the quadrilateral. + * \return Returns the vector. + */ +template +VTKM_EXEC Vector GetQuadN3(const CollectionOfPoints& pts) +{ + const Vector A = GetQuadL2(pts); + const Vector B = GetQuadL3(pts); + const Vector N3 = vtkm::Cross(A, B); + return N3; +} + +/** + * Returns the normal center vector, as defined by the verdict manual. + * + * \param [in] pts The four points which define the quadrilateral. + * \return Returns the vector. + */ +template +VTKM_EXEC Vector GetQuadNc(const CollectionOfPoints& pts) +{ + const Vector A = GetQuadX0(pts); + const Vector B = GetQuadX1(pts); + const Vector Nc = vtkm::Cross(A, B); + return Nc; +} + +/** + * Returns the normalized N0 vector, as defined by the verdict manual. + * + * \param [in] pts The four points which define the quadrilateral. + * \return Returns the vector. + */ +template +VTKM_EXEC Vector GetQuadN0Normalized(const CollectionOfPoints& pts) +{ + return vtkm::Normal(GetQuadN0(pts)); +} + +/** + * Returns the normalized N1 vector, as defined by the verdict manual. + * + * \param [in] pts The four points which define the quadrilateral. + * \return Returns the vector. + */ +template +VTKM_EXEC Vector GetQuadN1Normalized(const CollectionOfPoints& pts) +{ + return vtkm::Normal(GetQuadN1(pts)); +} + +/** + * Returns the normalized N2 vector, as defined by the verdict manual. + * + * \param [in] pts The four points which define the quadrilateral. + * \return Returns the vector. + */ +template +VTKM_EXEC Vector GetQuadN2Normalized(const CollectionOfPoints& pts) +{ + return vtkm::Normal(GetQuadN2(pts)); +} + +/** + * Returns the normalized N3 vector, as defined by the verdict manual. + * + * \param [in] pts The four points which define the quadrilateral. + * \return Returns the vector. + */ +template +VTKM_EXEC Vector GetQuadN3Normalized(const CollectionOfPoints& pts) +{ + return vtkm::Normal(GetQuadN3(pts)); +} + +/** + * Returns the normalized Nc vector, as defined by the verdict manual. + * + * \param [in] pts The four points which define the quadrilateral. + * \return Returns the vector. + */ +template +VTKM_EXEC Vector GetQuadNcNormalized(const CollectionOfPoints& pts) +{ + return vtkm::Normal(GetQuadNc(pts)); +} + +/** + * Returns the alpha0 scalar, as defined by the verdict manual. + * + * \param [in] pts The four points which define the quadrilateral. + * \return Returns the scalar. + */ +template +VTKM_EXEC Scalar GetQuadAlpha0(const CollectionOfPoints& pts) +{ + const Vector normalizedCenterNormal = + GetQuadNcNormalized(pts); + const Vector N0 = GetQuadN0(pts); + return vtkm::Dot(normalizedCenterNormal, N0); +} + +/** + * Returns the alpha1 scalar, as defined by the verdict manual. + * + * \param [in] pts The four points which define the quadrilateral. + * \return Returns the scalar. + */ +template +VTKM_EXEC Scalar GetQuadAlpha1(const CollectionOfPoints& pts) +{ + const Vector normalizedCenterNormal = + GetQuadNcNormalized(pts); + const Vector N1 = GetQuadN1(pts); + return vtkm::Dot(normalizedCenterNormal, N1); +} + +/** + * Returns the alpha2 scalar, as defined by the verdict manual. + * + * \param [in] pts The four points which define the quadrilateral. + * \return Returns the scalar. + */ +template +VTKM_EXEC Scalar GetQuadAlpha2(const CollectionOfPoints& pts) +{ + const Vector normalizedCenterNormal = + GetQuadNcNormalized(pts); + const Vector N2 = GetQuadN2(pts); + return vtkm::Dot(normalizedCenterNormal, N2); +} + + +/** + * Returns the alpha3 scalar, as defined by the verdict manual. + * + * \param [in] pts The four points which define the quadrilateral. + * \return Returns the scalar. + */ +template +VTKM_EXEC Scalar GetQuadAlpha3(const CollectionOfPoints& pts) +{ + const Vector normalizedCenterNormal = + GetQuadNcNormalized(pts); + const Vector N3 = GetQuadN3(pts); + return vtkm::Dot(normalizedCenterNormal, N3); +} + + +/** + * Returns the area of the quad, as defined by the verdict manual. + * + * \param [in] pts The four points which define the quadrilateral. + * \return Returns the area. + */ +template +VTKM_EXEC Scalar GetQuadArea(const CollectionOfPoints& pts) +{ + const Scalar quarter(0.25); + const Scalar a0 = GetQuadAlpha0(pts); + const Scalar a1 = GetQuadAlpha1(pts); + const Scalar a2 = GetQuadAlpha2(pts); + const Scalar a3 = GetQuadAlpha3(pts); + return quarter * (a0 + a1 + a2 + a3); +} + +#endif diff --git a/vtkm/exec/cellmetrics/TypeOfCellTetrahedral.h b/vtkm/exec/cellmetrics/TypeOfCellTetrahedral.h new file mode 100644 index 000000000..81baf3aea --- /dev/null +++ b/vtkm/exec/cellmetrics/TypeOfCellTetrahedral.h @@ -0,0 +1,295 @@ +#ifndef vtk_m_exec_cellmetrics_TypeOfCellTetrahedral +#define vtk_m_exec_cellmetrics_TypeOfCellTetrahedral +/** + * Returns the L0 vector, as defined by the verdict manual. + * + * \param [in] pts The four points which define the quadrilateral. + * \return Returns the vector. + */ +template +VTKM_EXEC Vector GetTetraL0(const CollectionOfPoints& pts) +{ + const Vector L0(pts[1] - pts[0]); + return L0; +} + +/** + * Returns the L1 vector, as defined by the verdict manual. + * + * \param [in] pts The four points which define the quadrilateral. + * \return Returns the vector. + */ +template +VTKM_EXEC Vector GetTetraL1(const CollectionOfPoints& pts) +{ + const Vector L1(pts[2] - pts[1]); + return L1; +} + +/** + * Returns the L2 vector, as defined by the verdict manual. + * + * \param [in] pts The four points which define the quadrilateral. + * \return Returns the vector. + */ +template +VTKM_EXEC Vector GetTetraL2(const CollectionOfPoints& pts) +{ + const Vector L2(pts[0] - pts[2]); + return L2; +} + +/** + * Returns the L3 vector, as defined by the verdict manual. + * + * \param [in] pts The four points which define the quadrilateral. + * \return Returns the vector. + */ +template +VTKM_EXEC Vector GetTetraL3(const CollectionOfPoints& pts) +{ + const Vector L3(pts[3] - pts[0]); + return L3; +} + +/** + * Returns the L4 vector, as defined by the verdict manual. + * + * \param [in] pts The four points which define the quadrilateral. + * \return Returns the vector. + */ +template +VTKM_EXEC Vector GetTetraL4(const CollectionOfPoints& pts) +{ + const Vector L4(pts[3] - pts[1]); + return L4; +} + +/** + * Returns the L5 vector, as defined by the verdict manual. + * + * \param [in] pts The four points which define the quadrilateral. + * \return Returns the vector. + */ +template +VTKM_EXEC Vector GetTetraL5(const CollectionOfPoints& pts) +{ + const Vector L5(pts[3] - pts[2]); + return L5; +} + +/** + * Returns the L0 vector's magnitude, as defined by the verdict manual. + * + * \param [in] pts The four points which define the Tetra. + * \return Returns the magnitude of the vector. + */ +template +VTKM_EXEC Scalar GetTetraL0Magnitude(const CollectionOfPoints& pts) +{ + const Scalar l0 = + vtkm::Sqrt(vtkm::MagnitudeSquared(GetTetraL0(pts))); + return l0; +} + +/** + * Returns the L1 vector's magnitude, as defined by the verdict manual. + * + * \param [in] pts The four points which define the Tetra. + * \return Returns the magnitude of the vector. + */ +template +VTKM_EXEC Scalar GetTetraL1Magnitude(const CollectionOfPoints& pts) +{ + const Scalar l1 = + vtkm::Sqrt(vtkm::MagnitudeSquared(GetTetraL1(pts))); + return l1; +} + +/** + * Returns the L2 vector's magnitude, as defined by the verdict manual. + * + * \param [in] pts The four points which define the Tetra. + * \return Returns the magnitude of the vector. + */ +template +VTKM_EXEC Scalar GetTetraL2Magnitude(const CollectionOfPoints& pts) +{ + const Scalar l2 = + vtkm::Sqrt(vtkm::MagnitudeSquared(GetTetraL2(pts))); + return l2; +} + +/** + * Returns the L3 vector's magnitude, as defined by the verdict manual. + * + * \param [in] pts The four points which define the Tetra. + * \return Returns the magnitude of the vector. + */ +template +VTKM_EXEC Scalar GetTetraL3Magnitude(const CollectionOfPoints& pts) +{ + const Scalar l3 = + vtkm::Sqrt(vtkm::MagnitudeSquared(GetTetraL3(pts))); + return l3; +} + +/** + * Returns the L4 vector's magnitude, as defined by the verdict manual. + * + * \param [in] pts The four points which define the Tetra. + * \return Returns the magnitude of the vector. + */ +template +VTKM_EXEC Scalar GetTetraL4Magnitude(const CollectionOfPoints& pts) +{ + const Scalar l4 = + vtkm::Sqrt(vtkm::MagnitudeSquared(GetTetraL4(pts))); + return l4; +} + +/** + * Returns the L5 vector's magnitude, as defined by the verdict manual. + * + * \param [in] pts The four points which define the Tetra. + * \return Returns the magnitude of the vector. + */ +template +VTKM_EXEC Scalar GetTetraL5Magnitude(const CollectionOfPoints& pts) +{ + const Scalar l5 = + vtkm::Sqrt(vtkm::MagnitudeSquared(GetTetraL5(pts))); + return l5; +} + +/** + * Returns the Max of the magnitude of each vector which makes up the edges of the Tetra. + * + * \param [in] pts The four points which define the Tetra. + * \return Returns the max. + */ +template +VTKM_EXEC Scalar GetTetraLMax(const CollectionOfPoints& pts) +{ + const Scalar l0 = GetTetraL0Magnitude(pts); + const Scalar l1 = GetTetraL1Magnitude(pts); + const Scalar l2 = GetTetraL2Magnitude(pts); + const Scalar l3 = GetTetraL3Magnitude(pts); + const Scalar l4 = GetTetraL4Magnitude(pts); + const Scalar l5 = GetTetraL5Magnitude(pts); + const Scalar lmax = vtkm::Max(l0, vtkm::Max(l1, vtkm::Max(l2, vtkm::Max(l3, vtkm::Max(l4, l5))))); + return lmax; +} + +/** + * Returns the Min of the magnitude of each vector which makes up the sides of the Tetra. + * + * \param [in] pts The four points which define the Tetra. + * \return Returns the min. + */ +template +VTKM_EXEC Scalar GetTetraLMin(const CollectionOfPoints& pts) +{ + const Scalar l0 = GetTetraL0Magnitude(pts); + const Scalar l1 = GetTetraL1Magnitude(pts); + const Scalar l2 = GetTetraL2Magnitude(pts); + const Scalar l3 = GetTetraL3Magnitude(pts); + const Scalar l4 = GetTetraL4Magnitude(pts); + const Scalar l5 = GetTetraL5Magnitude(pts); + const Scalar lmin = vtkm::Min(l0, vtkm::Min(l1, vtkm::Min(l2, vtkm::Min(l3, vtkm::Min(l4, l5))))); + return lmin; +} + +/** + * Returns the surface area of the Tetra. + * + * \param [in] pts The four points which define the Tetra. + * \return Returns the area. + */ +template +VTKM_EXEC Scalar GetTetraArea(const CollectionOfPoints& pts) +{ + const Vector L0 = GetTetraL0(pts); + const Vector L1 = GetTetraL1(pts); + const Vector L2 = GetTetraL2(pts); + const Vector L3 = GetTetraL3(pts); + const Vector L4 = GetTetraL4(pts); + + const Scalar a = vtkm::Sqrt(vtkm::MagnitudeSquared(vtkm::Cross(L2, L0))); + const Scalar b = vtkm::Sqrt(vtkm::MagnitudeSquared(vtkm::Cross(L3, L0))); + const Scalar c = vtkm::Sqrt(vtkm::MagnitudeSquared(vtkm::Cross(L4, L1))); + const Scalar d = vtkm::Sqrt(vtkm::MagnitudeSquared(vtkm::Cross(L3, L2))); + const Scalar half(0.5); + + const Scalar area = half * (a + b + c + d); + return area; +} + +/** + * Returns the volume of the Tetra. + * + * \param [in] pts The four points which define the Tetra. + * \return Returns the volume. + */ +template +VTKM_EXEC Scalar GetTetraVolume(const CollectionOfPoints& pts) +{ + const Vector L0 = GetTetraL0(pts); + const Vector L2 = GetTetraL2(pts); + const Vector L3 = GetTetraL3(pts); + const Scalar six(6.0); + + return vtkm::Dot(vtkm::Cross(L2, L0), L3) / six; +} + +/** + * Returns the inradius of the Tetra. + * + * \param [in] pts The four points which define the Tetra. + * \return Returns the inradius. + */ +template +VTKM_EXEC Scalar GetTetraInradius(const CollectionOfPoints& pts) +{ + const Scalar three(3.0); + const Scalar volume = GetTetraVolume(pts); + const Scalar area = GetTetraArea(pts); + + return (three * volume) / area; +} + + +/** + * Returns the circumradius of the Tetra. + * + * \param [in] pts The four points which define the Tetra. + * \return Returns the circumradius. + */ +template +VTKM_EXEC Scalar GetTetraCircumradius(const CollectionOfPoints& pts) +{ + const Vector L0 = GetTetraL0(pts); + const Vector L1 = GetTetraL1(pts); + const Vector L2 = GetTetraL2(pts); + const Vector L3 = GetTetraL3(pts); + const Vector L4 = GetTetraL4(pts); + + const Scalar l0l0 = vtkm::MagnitudeSquared(L0); + const Scalar l2l2 = vtkm::MagnitudeSquared(L2); + const Scalar l3l3 = vtkm::MagnitudeSquared(L3); + + const Vector A = l3l3 * vtkm::Cross(L2, L0); + const Vector B = l2l2 * vtkm::Cross(L3, L0); + const Vector C = l0l0 * vtkm::Cross(L3, L2); + + const Vector D(A + B + C); + const Scalar d = vtkm::Sqrt(vtkm::MagnitudeSquared(D)); + + const Scalar twelve(12.0); + const Scalar volume = GetTetraVolume(pts); + + const Scalar circumradius = d / (twelve * volume); + return circumradius; +} + +#endif diff --git a/vtkm/exec/cellmetrics/TypeOfCellTriangle.h b/vtkm/exec/cellmetrics/TypeOfCellTriangle.h new file mode 100644 index 000000000..3a7bd329e --- /dev/null +++ b/vtkm/exec/cellmetrics/TypeOfCellTriangle.h @@ -0,0 +1,192 @@ +#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 +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 +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 +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 +VTKM_EXEC Scalar GetTriangleL0Magnitude(const CollectionOfPoints& pts) +{ + const Scalar l0 = + vtkm::Sqrt(vtkm::MagnitudeSquared(GetTriangleL0(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 +VTKM_EXEC Scalar GetTriangleL1Magnitude(const CollectionOfPoints& pts) +{ + const Scalar l1 = + vtkm::Sqrt(vtkm::MagnitudeSquared(GetTriangleL1(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 +VTKM_EXEC Scalar GetTriangleL2Magnitude(const CollectionOfPoints& pts) +{ + const Scalar l2 = + vtkm::Sqrt(vtkm::MagnitudeSquared(GetTriangleL2(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 +VTKM_EXEC Scalar GetTriangleLMax(const CollectionOfPoints& pts) +{ + const Scalar l0 = GetTriangleL0Magnitude(pts); + const Scalar l1 = GetTriangleL1Magnitude(pts); + const Scalar l2 = GetTriangleL2Magnitude(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 +VTKM_EXEC Scalar GetTriangleLMin(const CollectionOfPoints& pts) +{ + const Scalar l0 = GetTriangleL0Magnitude(pts); + const Scalar l1 = GetTriangleL1Magnitude(pts); + const Scalar l2 = GetTriangleL2Magnitude(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 +VTKM_EXEC Scalar GetTriangleArea(const CollectionOfPoints& pts) +{ + const Vector L0 = GetTriangleL0(pts); + const Vector L1 = GetTriangleL1(pts); + const Scalar half(0.5); + const Scalar crossProductMagnitude = vtkm::Sqrt(vtkm::MagnitudeSquared(vtkm::Cross(L0, L1))); + const Scalar area = half * 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 +VTKM_EXEC Scalar GetTriangleInradius(const CollectionOfPoints& pts) +{ + const Scalar two(2.0); + const Scalar area = GetTriangleArea(pts); + const Scalar l0 = GetTriangleL0Magnitude(pts); + const Scalar l1 = GetTriangleL1Magnitude(pts); + const Scalar l2 = GetTriangleL2Magnitude(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 +VTKM_EXEC Scalar GetTriangleCircumradius(const CollectionOfPoints& pts) +{ + const Scalar four(4.0); + const Scalar area = GetTriangleArea(pts); + const Scalar l0 = GetTriangleL0Magnitude(pts); + const Scalar l1 = GetTriangleL1Magnitude(pts); + const Scalar l2 = GetTriangleL2Magnitude(pts); + const Scalar circumradius = (l0 * l1 * l2) / (four * area); + return circumradius; +} + +#endif diff --git a/vtkm/filter/testing/UnitTestMeshQualityFilter.cxx b/vtkm/filter/testing/UnitTestMeshQualityFilter.cxx index 7b2a886c6..c85851309 100644 --- a/vtkm/filter/testing/UnitTestMeshQualityFilter.cxx +++ b/vtkm/filter/testing/UnitTestMeshQualityFilter.cxx @@ -31,8 +31,8 @@ #include //Adapted from vtkm/cont/testing/MakeTestDataSet.h -//Modified the content of the Make3DExplicitDataSetZoo() function -inline vtkm::cont::DataSet Make3DExplicitDataSet() +//Modified the content of the MakeExplicitDataSetZoo() function +inline vtkm::cont::DataSet MakeExplicitDataSet() { vtkm::cont::DataSet dataSet; vtkm::cont::DataSetBuilderExplicit dsb; @@ -40,14 +40,11 @@ inline vtkm::cont::DataSet Make3DExplicitDataSet() using CoordType = vtkm::Vec3f_64; std::vector coords = { - { 0.00, 0.00, 0.00 }, { 1.00, 0.00, 0.00 }, { 2.00, 0.00, 0.00 }, { 0.00, 0.00, 1.00 }, - { 1.00, 0.00, 1.00 }, { 2.00, 0.00, 1.00 }, { 0.00, 1.00, 0.00 }, { 1.00, 1.00, 0.00 }, - { 2.00, 1.00, 0.00 }, { 0.00, 1.00, 1.00 }, { 1.00, 1.00, 1.00 }, { 2.00, 1.00, 1.00 }, - { 0.00, 2.00, 0.00 }, { 1.00, 2.00, 0.00 }, { 2.00, 2.00, 0.00 }, { 0.00, 2.00, 1.00 }, - { 1.00, 2.00, 1.00 }, { 2.00, 2.00, 1.00 }, { 1.00, 3.00, 1.00 }, { 2.75, 0.00, 1.00 }, - { 3.00, 0.00, 0.75 }, { 3.00, 0.25, 1.00 }, { 3.00, 1.00, 1.00 }, { 3.00, 1.00, 0.00 }, - { 2.57, 2.00, 1.00 }, { 3.00, 1.75, 1.00 }, { 3.00, 1.75, 0.75 }, { 3.00, 0.00, 0.00 }, - { 2.57, 0.42, 0.57 }, { 2.59, 1.43, 0.71 } + { 0, 0, 0 }, { 3, 0, 0 }, { 2, 2, 0 }, { 4, 0, 0 }, { 7, 0, 0 }, { 7, 2, 0 }, + { 6, 2, 0 }, { 8, 0, 0 }, { 11, 0, 0 }, { 9, 2, 0 }, { 9, 1, 1 }, { 9, 3, 0 }, + { 11, 3, 0 }, { 11, 5, 0 }, { 9, 5, 0 }, { 10, 4, 1 }, { 12, 0, 0 }, { 12, 3, 0 }, + { 12, 2, 1 }, { 15, 0, 0 }, { 15, 3, 0 }, { 15, 1, 1 }, { 16, 0, 0 }, { 18, 0, 0 }, + { 18, 2, 0 }, { 16, 2, 0 }, { 17, 1, 1 }, { 19, 1, 1 }, { 19, 3, 1 }, { 17, 3, 1 } }; std::vector shapes; @@ -56,145 +53,64 @@ inline vtkm::cont::DataSet Make3DExplicitDataSet() //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); + shapes.push_back(vtkm::CELL_SHAPE_TRIANGLE); + numindices.push_back(3); conn.push_back(0); + conn.push_back(1); + conn.push_back(2); + + shapes.push_back(vtkm::CELL_SHAPE_QUAD); + numindices.push_back(4); conn.push_back(3); conn.push_back(4); - conn.push_back(1); + conn.push_back(5); conn.push_back(6); + + shapes.push_back(vtkm::CELL_SHAPE_TETRA); + numindices.push_back(4); + conn.push_back(7); + conn.push_back(8); conn.push_back(9); conn.push_back(10); - conn.push_back(7); + + shapes.push_back(vtkm::CELL_SHAPE_PYRAMID); + numindices.push_back(5); + conn.push_back(11); + conn.push_back(12); + conn.push_back(13); + conn.push_back(14); + conn.push_back(15); + + shapes.push_back(vtkm::CELL_SHAPE_WEDGE); + numindices.push_back(6); + conn.push_back(16); + conn.push_back(17); + conn.push_back(18); + conn.push_back(19); + conn.push_back(20); + conn.push_back(21); 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(23); + conn.push_back(24); + conn.push_back(25); + conn.push_back(26); + conn.push_back(27); + conn.push_back(28); 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"); return dataSet; } -template -std::vector TestMeshQualityFilter(const vtkm::cont::DataSet& input, - const std::vector& expectedVals, - const std::string& outputname, - T filter) +bool TestMeshQualityFilter(const vtkm::cont::DataSet& input, + const std::vector& expectedVals, + const std::string& outputname, + vtkm::filter::MeshQuality& filter) { - std::vector errors; vtkm::cont::DataSet output; try { @@ -202,168 +118,93 @@ std::vector TestMeshQualityFilter(const vtkm::cont::DataSet& input, } catch (vtkm::cont::ErrorExecution&) { - errors.push_back("Error occured while executing filter. Exiting..."); - return errors; + return true; } //Test the computed metric values (for all cells) and expected metric //values for equality. - vtkm::cont::testing::TestEqualResult result = vtkm::cont::testing::test_equal_ArrayHandles( - vtkm::cont::make_ArrayHandle(expectedVals), output.GetField(outputname).GetData()); - if (!result) - result.PushMessage(std::string("Data doesn't match")); - - return result.GetMessages(); -} - -//Either an error occurred during execution of the -//filter, or mismatches exist between the expected output -//and the computed filter output. -void CheckForErrors(const std::vector& messages) -{ - if (!messages.empty()) + vtkm::cont::ArrayHandle values; + output.GetField(outputname).GetData().CopyTo(values); + auto portal1 = values.GetPortalConstControl(); + if (portal1.GetNumberOfValues() != (vtkm::Id)expectedVals.size()) { - std::cout << "FAIL\n"; - for (std::string m : messages) - std::cout << m << "\n"; + printf("Number of expected values for %s does not match.\n", outputname.c_str()); + return true; } - else - std::cout << "SUCCESS\n"; + + std::vector cellTypes = { "triangle", "quadrilateral", "tetrahedron", + "pyramid", "wedge", "hexahedron" }; + bool anyFailures = false; + for (unsigned long i = 0; i < expectedVals.size(); i++) + { + vtkm::Id id = (vtkm::Id)i; + if (portal1.Get(id) != expectedVals[i]) + { + printf("Metric \"%s\" for cell type \"%s\" does not match. Expected %f and got %f\n", + outputname.c_str(), + cellTypes[i].c_str(), + expectedVals[i], + portal1.Get(id)); + anyFailures = true; + } + } + return anyFailures; } int TestMeshQuality() { using FloatVec = std::vector; - using PairVec = std::vector>; - using StringVec = std::vector; - using CharVec = std::vector; - using QualityFilter = vtkm::filter::MeshQuality; //Test variables - vtkm::cont::DataSet input = Make3DExplicitDataSet(); - std::unique_ptr filter; - std::string metricSuffix; - StringVec fieldNames; - CharVec testedShapes; - PairVec shapeMetricPairs; - FloatVec expectedValues; + vtkm::cont::DataSet input = MakeExplicitDataSet(); - /*************************************************** - * Test 1: Volume metric - ***************************************************/ + int numFailures = 0; + bool testFailed = false; - std::cout << "Testing MeshQuality filter: Volume metric" - << "\n++++++++++++++++++++++++++++++++++++++++++++++++++\n"; + std::vector expectedValues; + std::vector metrics; + std::vector metricName; - //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, - (vtkm::FloatDefault)0.0100042, - (vtkm::FloatDefault)0.0983333, - (vtkm::FloatDefault)0.0732667, - (vtkm::FloatDefault)0.0845833, - -0.5, - -0.5, - 0, - 0, - 1, - 1, - 1.5, - (vtkm::FloatDefault)0.7071068, - 2, - 1, - 0.5, - 1 }; + /* + FloatVec volumeExpectedValues = { 0, 0, 1, (float) 1.333333333, 4, 4 }; + expectedValues.push_back(volumeExpectedValues); + metrics.push_back(vtkm::filter::CellMetric::VOLUME); + metricName.push_back("volume"); + */ - filter.reset(new QualityFilter(vtkm::filter::CellMetric::VOLUME)); - std::vector errors = - TestMeshQualityFilter(input, expectedValues, "volume", *filter); - std::cout << "Volume metric test: "; - CheckForErrors(errors); + FloatVec jacobianExpectedValues = { 0, 2, 6, 0, 0, 4 }; + expectedValues.push_back(jacobianExpectedValues); + metrics.push_back(vtkm::filter::CellMetric::JACOBIAN); + metricName.push_back("jacobian"); - /*************************************************** - * Test 2: Diagonal Ratio metric - ***************************************************/ + /* + FloatVec diagonalRatioExpectedValues = { -1, -1, -1, -1, -1, (float) 0.39 }; + expectedValues.push_back(diagonalRatioExpectedValues); + metrics.push_back(vtkm::filter::CellMetric::DIAGONAL_RATIO); + metricName.push_back("diagonalRatio"); + */ - std::cout << "Testing MeshQuality filter: Diagonal Ratio metric" - << "\n++++++++++++++++++++++++++++++++++++++++++++++++++\n"; - - - expectedValues = { - 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, (vtkm::FloatDefault)2.23607 - }; - - filter.reset(new QualityFilter(vtkm::filter::CellMetric::DIAGONAL_RATIO)); - errors = TestMeshQualityFilter(input, expectedValues, "diagonalRatio", *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, (vtkm::FloatDefault)2.23607 - }; - - filter.reset(new QualityFilter(vtkm::filter::CellMetric::JACOBIAN)); - errors = TestMeshQualityFilter(input, expectedValues, "jacobian", *filter); - std::cout << "Jacobian metric test: "; - CheckForErrors(errors); - -#if 0 - /*************************************************** - * Test 5: Oddy metric - ***************************************************/ - - std::cout << "Testing MeshQuality filter: Oddy metric" - << "\n++++++++++++++++++++++++++++++++++++++++++++++++++\n"; - - testedShapes = {vtkm::CELL_SHAPE_HEXAHEDRON, vtkm::CELL_SHAPE_POLYGON, - vtkm::CELL_SHAPE_QUAD}; - - shapeMetricPairs.clear(); - for (auto s : testedShapes) - shapeMetricPairs.push_back(vtkm::make_Pair(s, vtkm::filter::CellMetric::ODDY)); - expectedValues = {0, 0, 0, 0, 0, 0, 0, - 0, 0, 0, 0, 0, 0, 0, 1.125, 0, - 0, 2.5}; /*(all cells, volume value)*/ - - filter.reset(new QualityFilter(shapeMetricPairs)); - errors = TestMeshQualityFilter(input, expectedValues, *filter); - std::cout << "Oddy metric test: "; - CheckForErrors(errors); - - /*************************************************** - * Test 4: Relative Size Squared metric - ***************************************************/ - - std::cout << "Testing MeshQuality filter: Relative Size Squared metric" - << "\n++++++++++++++++++++++++++++++++++++++++++++++++++\n"; - - testedShapes = {vtkm::CELL_SHAPE_HEXAHEDRON, vtkm::CELL_SHAPE_POLYGON, - vtkm::CELL_SHAPE_QUAD, vtkm::CELL_SHAPE_TRIANGLE, - vtkm::CELL_SHAPE_TETRA}; - - shapeMetricPairs.clear(); - for (auto s : testedShapes) - shapeMetricPairs.push_back(vtkm::make_Pair(s, vtkm::filter::CellMetric::RELATIVE_SIZE)); - expectedValues = {1, 1, 0.0341086, 0.303456, 0, 0, 0, - 0, 0, 0, 0, 0, 0.361898, 0.614047, 1, 1, - 0.307024, 1}; - - filter.reset(new QualityFilter(shapeMetricPairs)); - errors = TestMeshQualityFilter(input, expectedValues, *filter); - std::cout << "Relative Size Square metric test: "; - CheckForErrors(errors); -#endif + unsigned long numTests = metrics.size(); + for (unsigned long i = 0; i < numTests; i++) + { + printf("Testing metric %s\n", metricName[i].c_str()); + vtkm::filter::MeshQuality filter(metrics[i]); + testFailed = TestMeshQualityFilter(input, expectedValues[i], metricName[i], filter); + if (testFailed) + { + numFailures++; + printf("\ttest \"%s\" failed\n", metricName[i].c_str()); + } + else + printf("\t... passed\n"); + } + if (numFailures > 0) + { + printf("Number of failed metrics is %d\n", numFailures); + bool see_previous_messages = false; // this variable name plays well with macro + VTKM_TEST_ASSERT(see_previous_messages, "Failure occurred during test"); + } return 0; }