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

Jacobian working, unit test working

Browse Source
This commit is contained in:
Hank Childs 2019-09-17 13:41:53 -07:00
parent 560aa1f673
commit ccedb5aac1
6 changed files with 1900 additions and 430 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

223
vtkm/exec/cellmetrics/CellJacobianMetric.h
View File

@ -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 <typename OutType, typename VecType>
VTKM_EXEC inline OutType CellJacobianMetricOfQuad(const VecType& edgeCalculations,
const VecType& axes)
{
const vtkm::Id numCalculations = edgeCalculations.GetNumberOfComponents();
//Compare partitions of quad to find min
using axesType = typename VecType::ComponentType;
axesType centerCalculation = vtkm::Cross(axes[0], axes[1]);
vtkm::Normalize(centerCalculation);
OutType currCalculation, minCalculation = vtkm::Infinity<OutType>();
for (vtkm::IdComponent i = 0; i < numCalculations; i++)
{
currCalculation = vtkm::Dot(edgeCalculations[i], centerCalculation);
if (currCalculation < minCalculation)
minCalculation = currCalculation;
}
if (minCalculation > 0)
return vtkm::Min(minCalculation, vtkm::Infinity<OutType>()); //normal case
return vtkm::Max(minCalculation, OutType(-1) * vtkm::Infinity<OutType>());
}
template <typename OutType, typename VecType>
VTKM_EXEC inline OutType CellJacobianMetricOfHex(const VecType& matrices)
{
const vtkm::IdComponent numMatrices = matrices.GetNumberOfComponents();
//Compare determinants to find min
OutType currDeterminant, minDeterminant;
//principle axes matrix computed outside of for loop to avoid un-necessary if statement
minDeterminant =
(OutType)vtkm::Dot(matrices[numMatrices - 1][0],
vtkm::Cross(matrices[numMatrices - 1][1], matrices[numMatrices - 1][2]));
minDeterminant /= 64.0;
for (vtkm::IdComponent i = 0; i < numMatrices - 1; i++)
{
currDeterminant =
(OutType)vtkm::Dot(matrices[i][0], vtkm::Cross(matrices[i][1], matrices[i][2]));
if (currDeterminant < minDeterminant)
minDeterminant = currDeterminant;
}
if (minDeterminant > 0)
return vtkm::Min(minDeterminant, vtkm::Infinity<OutType>()); //normal case
return vtkm::Max(minDeterminant, vtkm::NegativeInfinity<OutType>());
}
// ========================= Unsupported cells ==================================
// By default, cells have zero shape unless the shape type template is specialized below.
// 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,
@ -116,63 +66,9 @@ VTKM_EXEC OutType CellJacobianMetric(const vtkm::IdComponent& numPts,
return OutType(0.0);
}
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellJacobianMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagPolygon,
const vtkm::exec::FunctorBase& worklet)
{
switch (numPts)
{
case 4:
return CellJacobianMetric<OutType>(numPts, pts, vtkm::CellShapeTagQuad(), worklet);
default:
break;
}
return OutType(-1.0);
}
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellJacobianMetric(const vtkm::IdComponent&,
const PointCoordVecType&,
vtkm::CellShapeTagLine,
const vtkm::exec::FunctorBase& worklet)
{
UNUSED(worklet);
return OutType(-1.0);
}
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellJacobianMetric(const vtkm::IdComponent&,
const PointCoordVecType&,
vtkm::CellShapeTagTriangle,
const vtkm::exec::FunctorBase& worklet)
{
UNUSED(worklet);
return OutType(-1.0);
}
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellJacobianMetric(const vtkm::IdComponent&,
const PointCoordVecType&,
vtkm::CellShapeTagWedge,
const vtkm::exec::FunctorBase& worklet)
{
UNUSED(worklet);
return OutType(-1.0);
}
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellJacobianMetric(const vtkm::IdComponent&,
const PointCoordVecType&,
vtkm::CellShapeTagPyramid,
const vtkm::exec::FunctorBase& worklet)
{
UNUSED(worklet);
return OutType(-1.0);
}
// ========================= 2D cells ==================================
// Compute the jacobian of a quadrilateral.
// 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<OutType>(
vtkm::make_VecC(QuadEdgesToUse, 4), vtkm::make_VecC(QuadAxes, 2));
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}, 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<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 Edge hexMatrices[9][3] = { { HexEdges[0], HexEdges[3], HexEdges[4] },
{ HexEdges[1], -1 * HexEdges[0], HexEdges[5] },
{ HexEdges[2], -1 * HexEdges[1], HexEdges[6] },
{ -1 * HexEdges[3], -1 * HexEdges[2], HexEdges[7] },
{ HexEdges[11], HexEdges[8], -1 * HexEdges[4] },
{ -1 * HexEdges[8], HexEdges[9], -1 * HexEdges[5] },
{ -1 * HexEdges[9], HexEdges[10], -1 * HexEdges[6] },
{ -1 * HexEdges[10], -1 * HexEdges[11], -1 * HexEdges[7] },
{ principleXAxis, principleYAxis, principleZAxis } };
return vtkm::exec::cellmetrics::CellJacobianMetricOfHex<OutType>(
vtkm::make_VecC(hexMatrices, 12));
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<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_CellEdgeRatioMetric_h
#endif // vtk_m_exec_cellmetrics_CellJacobianMetric_h

765
vtkm/exec/cellmetrics/TypeOfCellHexahedral.h Normal file
View File

@ -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 <typename Scalar, typename Vector, typename CollectionOfPoints>
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 <typename Scalar, typename Vector, typename CollectionOfPoints>
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 <typename Scalar, typename Vector, typename CollectionOfPoints>
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 <typename Scalar, typename Vector, typename CollectionOfPoints>
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 <typename Scalar, typename Vector, typename CollectionOfPoints>
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 <typename Scalar, typename Vector, typename CollectionOfPoints>
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 <typename Scalar, typename Vector, typename CollectionOfPoints>
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 <typename Scalar, typename Vector, typename CollectionOfPoints>
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 <typename Scalar, typename Vector, typename CollectionOfPoints>
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 <typename Scalar, typename Vector, typename CollectionOfPoints>
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 <typename Scalar, typename Vector, typename CollectionOfPoints>
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 <typename Scalar, typename Vector, typename CollectionOfPoints>
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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetHexL0Magnitude(const CollectionOfPoints& pts)
{
const Scalar l0 =
vtkm::Sqrt(vtkm::MagnitudeSquared(GetHexL0<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetHexL1Magnitude(const CollectionOfPoints& pts)
{
const Scalar l1 =
vtkm::Sqrt(vtkm::MagnitudeSquared(GetHexL1<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetHexL2Magnitude(const CollectionOfPoints& pts)
{
const Scalar l2 =
vtkm::Sqrt(vtkm::MagnitudeSquared(GetHexL2<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetHexL3Magnitude(const CollectionOfPoints& pts)
{
const Scalar l3 =
vtkm::Sqrt(vtkm::MagnitudeSquared(GetHexL3<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetHexL4Magnitude(const CollectionOfPoints& pts)
{
const Scalar l4 =
vtkm::Sqrt(vtkm::MagnitudeSquared(GetHexL4<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetHexL5Magnitude(const CollectionOfPoints& pts)
{
const Scalar l5 =
vtkm::Sqrt(vtkm::MagnitudeSquared(GetHexL5<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetHexL6Magnitude(const CollectionOfPoints& pts)
{
const Scalar l6 =
vtkm::Sqrt(vtkm::MagnitudeSquared(GetHexL6<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetHexL7Magnitude(const CollectionOfPoints& pts)
{
const Scalar l7 =
vtkm::Sqrt(vtkm::MagnitudeSquared(GetHexL7<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetHexL8Magnitude(const CollectionOfPoints& pts)
{
const Scalar l8 =
vtkm::Sqrt(vtkm::MagnitudeSquared(GetHexL8<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetHexL9Magnitude(const CollectionOfPoints& pts)
{
const Scalar l9 =
vtkm::Sqrt(vtkm::MagnitudeSquared(GetHexL9<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetHexL10Magnitude(const CollectionOfPoints& pts)
{
const Scalar l10 =
vtkm::Sqrt(vtkm::MagnitudeSquared(GetHexL10<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetHexL11Magnitude(const CollectionOfPoints& pts)
{
const Scalar l11 =
vtkm::Sqrt(vtkm::MagnitudeSquared(GetHexL11<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetHexLMax(const CollectionOfPoints& pts)
{
const Scalar l0 = GetHexL0Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l1 = GetHexL1Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l2 = GetHexL2Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l3 = GetHexL3Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l4 = GetHexL4Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l5 = GetHexL5Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l6 = GetHexL6Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l7 = GetHexL7Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l8 = GetHexL8Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l9 = GetHexL9Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l10 = GetHexL10Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l11 = GetHexL11Magnitude<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetHexLMin(const CollectionOfPoints& pts)
{
const Scalar l0 = GetHexL0Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l1 = GetHexL1Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l2 = GetHexL2Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l3 = GetHexL3Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l4 = GetHexL4Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l5 = GetHexL5Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l6 = GetHexL6Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l7 = GetHexL7Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l8 = GetHexL8Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l9 = GetHexL9Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l10 = GetHexL10Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l11 = GetHexL11Magnitude<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
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 <typename Scalar, typename Vector, typename CollectionOfPoints>
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 <typename Scalar, typename Vector, typename CollectionOfPoints>
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 <typename Scalar, typename Vector, typename CollectionOfPoints>
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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetHexD0Magnitude(const CollectionOfPoints& pts)
{
const Scalar d0 =
vtkm::Sqrt(vtkm::MagnitudeSquared(GetHexD0<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetHexD1Magnitude(const CollectionOfPoints& pts)
{
const Scalar d1 =
vtkm::Sqrt(vtkm::MagnitudeSquared(GetHexD1<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetHexD2Magnitude(const CollectionOfPoints& pts)
{
const Scalar d2 =
vtkm::Sqrt(vtkm::MagnitudeSquared(GetHexD2<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetHexD3Magnitude(const CollectionOfPoints& pts)
{
const Scalar d3 =
vtkm::Sqrt(vtkm::MagnitudeSquared(GetHexD3<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetHexDMin(const CollectionOfPoints& pts)
{
const Scalar d0 = GetHexD0Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar d1 = GetHexD1Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar d2 = GetHexD2Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar d3 = GetHexD3Magnitude<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetHexDMax(const CollectionOfPoints& pts)
{
const Scalar d0 = GetHexD0Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar d1 = GetHexD1Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar d2 = GetHexD2Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar d3 = GetHexD3Magnitude<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
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 <typename Scalar, typename Vector, typename CollectionOfPoints>
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 <typename Scalar, typename Vector, typename CollectionOfPoints>
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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC vtkm::Vec<Vector, 3> GetHexAi(const CollectionOfPoints& pts, const vtkm::Id& index)
{
const Scalar neg1(-1.0);
if (index == 0)
{
const Vector v0 = GetHexL0<Scalar, Vector, CollectionOfPoints>(pts);
const Vector v1 = GetHexL3<Scalar, Vector, CollectionOfPoints>(pts);
const Vector v2 = GetHexL4<Scalar, Vector, CollectionOfPoints>(pts);
const vtkm::Vec<Vector, 3> A = { v0, v1, v2 };
return A;
}
else if (index == 1)
{
const Vector v0 = GetHexL1<Scalar, Vector, CollectionOfPoints>(pts);
const Vector v1 = neg1 * GetHexL0<Scalar, Vector, CollectionOfPoints>(pts);
const Vector v2 = GetHexL5<Scalar, Vector, CollectionOfPoints>(pts);
const vtkm::Vec<Vector, 3> A = { v0, v1, v2 };
return A;
}
else if (index == 2)
{
const Vector v0 = GetHexL2<Scalar, Vector, CollectionOfPoints>(pts);
const Vector v1 = neg1 * GetHexL1<Scalar, Vector, CollectionOfPoints>(pts);
const Vector v2 = GetHexL6<Scalar, Vector, CollectionOfPoints>(pts);
const vtkm::Vec<Vector, 3> A = { v0, v1, v2 };
return A;
}
else if (index == 3)
{
const Vector v0 = neg1 * GetHexL3<Scalar, Vector, CollectionOfPoints>(pts);
const Vector v1 = neg1 * GetHexL2<Scalar, Vector, CollectionOfPoints>(pts);
const Vector v2 = GetHexL7<Scalar, Vector, CollectionOfPoints>(pts);
const vtkm::Vec<Vector, 3> A = { v0, v1, v2 };
return A;
}
else if (index == 4)
{
const Vector v0 = GetHexL11<Scalar, Vector, CollectionOfPoints>(pts);
const Vector v1 = GetHexL8<Scalar, Vector, CollectionOfPoints>(pts);
const Vector v2 = neg1 * GetHexL4<Scalar, Vector, CollectionOfPoints>(pts);
const vtkm::Vec<Vector, 3> A = { v0, v1, v2 };
return A;
}
else if (index == 5)
{
const Vector v0 = neg1 * GetHexL8<Scalar, Vector, CollectionOfPoints>(pts);
const Vector v1 = GetHexL9<Scalar, Vector, CollectionOfPoints>(pts);
const Vector v2 = neg1 * GetHexL5<Scalar, Vector, CollectionOfPoints>(pts);
const vtkm::Vec<Vector, 3> A = { v0, v1, v2 };
return A;
}
else if (index == 6)
{
const Vector v0 = neg1 * GetHexL9<Scalar, Vector, CollectionOfPoints>(pts);
const Vector v1 = GetHexL10<Scalar, Vector, CollectionOfPoints>(pts);
const Vector v2 = neg1 * GetHexL6<Scalar, Vector, CollectionOfPoints>(pts);
const vtkm::Vec<Vector, 3> A = { v0, v1, v2 };
return A;
}
else if (index == 7)
{
const Vector v0 = neg1 * GetHexL10<Scalar, Vector, CollectionOfPoints>(pts);
const Vector v1 = neg1 * GetHexL11<Scalar, Vector, CollectionOfPoints>(pts);
const Vector v2 = neg1 * GetHexL7<Scalar, Vector, CollectionOfPoints>(pts);
const vtkm::Vec<Vector, 3> A = { v0, v1, v2 };
return A;
}
else
{
const Vector v0 = GetHexX1<Scalar, Vector, CollectionOfPoints>(pts);
const Vector v1 = GetHexX2<Scalar, Vector, CollectionOfPoints>(pts);
const Vector v2 = GetHexX3<Scalar, Vector, CollectionOfPoints>(pts);
const vtkm::Vec<Vector, 3> 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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetHexAiNormSquared(const CollectionOfPoints& pts, const vtkm::Id& index)
{
const vtkm::Vec<Vector, 3> Ai = GetHexAi<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetHexAiAdjNormSquared(const CollectionOfPoints& pts, const vtkm::Id& index)
{
const vtkm::Vec<Vector, 3> Ai = GetHexAi<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetHexAlphai(const CollectionOfPoints& pts, const vtkm::Id& index)
{
const vtkm::Vec<Vector, 3> Ai = GetHexAi<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC vtkm::Vec<Vector, 3> GetHexAiHat(const CollectionOfPoints& pts, const vtkm::Id& index)
{
const vtkm::Vec<Vector, 3> Ai = GetHexAi<Scalar, Vector, CollectionOfPoints>(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<Vector, 3> 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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetHexAlphaiHat(const CollectionOfPoints& pts, const vtkm::Id& index)
{
const vtkm::Vec<Vector, 3> Ai = GetHexAiHat<Scalar, Vector, CollectionOfPoints>(pts, index);
const Scalar hatAlpha_i = vtkm::Dot(Ai[0], vtkm::Cross(Ai[1], Ai[2]));
return hatAlpha_i;
}
#endif

468
vtkm/exec/cellmetrics/TypeOfCellQuadrilateral.h Normal file
View File

@ -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 <typename Scalar, typename Vector, typename CollectionOfPoints>
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 <typename Scalar, typename Vector, typename CollectionOfPoints>
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 <typename Scalar, typename Vector, typename CollectionOfPoints>
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 <typename Scalar, typename Vector, typename CollectionOfPoints>
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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetQuadL0Magnitude(const CollectionOfPoints& pts)
{
const Scalar l0 =
vtkm::Sqrt(vtkm::MagnitudeSquared(GetQuadL0<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetQuadL1Magnitude(const CollectionOfPoints& pts)
{
const Scalar l1 =
vtkm::Sqrt(vtkm::MagnitudeSquared(GetQuadL1<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetQuadL2Magnitude(const CollectionOfPoints& pts)
{
const Scalar l2 =
vtkm::Sqrt(vtkm::MagnitudeSquared(GetQuadL2<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetQuadL3Magnitude(const CollectionOfPoints& pts)
{
const Scalar l3 =
vtkm::Sqrt(vtkm::MagnitudeSquared(GetQuadL3<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetQuadLMax(const CollectionOfPoints& pts)
{
const Scalar l0 = GetQuadL0Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l1 = GetQuadL1Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l2 = GetQuadL2Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l3 = GetQuadL3Magnitude<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetQuadLMin(const CollectionOfPoints& pts)
{
const Scalar l0 = GetQuadL0Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l1 = GetQuadL1Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l2 = GetQuadL2Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l3 = GetQuadL3Magnitude<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
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 <typename Scalar, typename Vector, typename CollectionOfPoints>
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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetQuadD0Magnitude(const CollectionOfPoints& pts)
{
const Scalar d0 =
vtkm::Sqrt(vtkm::MagnitudeSquared(GetQuadD0<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetQuadD1Magnitude(const CollectionOfPoints& pts)
{
const Scalar d1 =
vtkm::Sqrt(vtkm::MagnitudeSquared(GetQuadD1<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetQuadDMax(const CollectionOfPoints& pts)
{
const Scalar d0 = GetQuadD0Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar d1 = GetQuadD1Magnitude<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
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 <typename Scalar, typename Vector, typename CollectionOfPoints>
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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Vector GetQuadN0(const CollectionOfPoints& pts)
{
const Vector A = GetQuadL3<Scalar, Vector, CollectionOfPoints>(pts);
const Vector B = GetQuadL0<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Vector GetQuadN1(const CollectionOfPoints& pts)
{
const Vector A = GetQuadL0<Scalar, Vector, CollectionOfPoints>(pts);
const Vector B = GetQuadL1<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Vector GetQuadN2(const CollectionOfPoints& pts)
{
const Vector A = GetQuadL1<Scalar, Vector, CollectionOfPoints>(pts);
const Vector B = GetQuadL2<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Vector GetQuadN3(const CollectionOfPoints& pts)
{
const Vector A = GetQuadL2<Scalar, Vector, CollectionOfPoints>(pts);
const Vector B = GetQuadL3<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Vector GetQuadNc(const CollectionOfPoints& pts)
{
const Vector A = GetQuadX0<Scalar, Vector, CollectionOfPoints>(pts);
const Vector B = GetQuadX1<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Vector GetQuadN0Normalized(const CollectionOfPoints& pts)
{
return vtkm::Normal(GetQuadN0<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Vector GetQuadN1Normalized(const CollectionOfPoints& pts)
{
return vtkm::Normal(GetQuadN1<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Vector GetQuadN2Normalized(const CollectionOfPoints& pts)
{
return vtkm::Normal(GetQuadN2<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Vector GetQuadN3Normalized(const CollectionOfPoints& pts)
{
return vtkm::Normal(GetQuadN3<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Vector GetQuadNcNormalized(const CollectionOfPoints& pts)
{
return vtkm::Normal(GetQuadNc<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetQuadAlpha0(const CollectionOfPoints& pts)
{
const Vector normalizedCenterNormal =
GetQuadNcNormalized<Scalar, Vector, CollectionOfPoints>(pts);
const Vector N0 = GetQuadN0<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetQuadAlpha1(const CollectionOfPoints& pts)
{
const Vector normalizedCenterNormal =
GetQuadNcNormalized<Scalar, Vector, CollectionOfPoints>(pts);
const Vector N1 = GetQuadN1<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetQuadAlpha2(const CollectionOfPoints& pts)
{
const Vector normalizedCenterNormal =
GetQuadNcNormalized<Scalar, Vector, CollectionOfPoints>(pts);
const Vector N2 = GetQuadN2<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetQuadAlpha3(const CollectionOfPoints& pts)
{
const Vector normalizedCenterNormal =
GetQuadNcNormalized<Scalar, Vector, CollectionOfPoints>(pts);
const Vector N3 = GetQuadN3<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetQuadArea(const CollectionOfPoints& pts)
{
const Scalar quarter(0.25);
const Scalar a0 = GetQuadAlpha0<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar a1 = GetQuadAlpha1<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar a2 = GetQuadAlpha2<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar a3 = GetQuadAlpha3<Scalar, Vector, CollectionOfPoints>(pts);
return quarter * (a0 + a1 + a2 + a3);
}
#endif

295
vtkm/exec/cellmetrics/TypeOfCellTetrahedral.h Normal file
View File

@ -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 <typename Scalar, typename Vector, typename CollectionOfPoints>
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 <typename Scalar, typename Vector, typename CollectionOfPoints>
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 <typename Scalar, typename Vector, typename CollectionOfPoints>
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 <typename Scalar, typename Vector, typename CollectionOfPoints>
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 <typename Scalar, typename Vector, typename CollectionOfPoints>
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 <typename Scalar, typename Vector, typename CollectionOfPoints>
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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetTetraL0Magnitude(const CollectionOfPoints& pts)
{
const Scalar l0 =
vtkm::Sqrt(vtkm::MagnitudeSquared(GetTetraL0<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetTetraL1Magnitude(const CollectionOfPoints& pts)
{
const Scalar l1 =
vtkm::Sqrt(vtkm::MagnitudeSquared(GetTetraL1<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetTetraL2Magnitude(const CollectionOfPoints& pts)
{
const Scalar l2 =
vtkm::Sqrt(vtkm::MagnitudeSquared(GetTetraL2<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetTetraL3Magnitude(const CollectionOfPoints& pts)
{
const Scalar l3 =
vtkm::Sqrt(vtkm::MagnitudeSquared(GetTetraL3<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetTetraL4Magnitude(const CollectionOfPoints& pts)
{
const Scalar l4 =
vtkm::Sqrt(vtkm::MagnitudeSquared(GetTetraL4<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetTetraL5Magnitude(const CollectionOfPoints& pts)
{
const Scalar l5 =
vtkm::Sqrt(vtkm::MagnitudeSquared(GetTetraL5<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetTetraLMax(const CollectionOfPoints& pts)
{
const Scalar l0 = GetTetraL0Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l1 = GetTetraL1Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l2 = GetTetraL2Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l3 = GetTetraL3Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l4 = GetTetraL4Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l5 = GetTetraL5Magnitude<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetTetraLMin(const CollectionOfPoints& pts)
{
const Scalar l0 = GetTetraL0Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l1 = GetTetraL1Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l2 = GetTetraL2Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l3 = GetTetraL3Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l4 = GetTetraL4Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l5 = GetTetraL5Magnitude<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetTetraArea(const CollectionOfPoints& pts)
{
const Vector L0 = GetTetraL0<Scalar, Vector, CollectionOfPoints>(pts);
const Vector L1 = GetTetraL1<Scalar, Vector, CollectionOfPoints>(pts);
const Vector L2 = GetTetraL2<Scalar, Vector, CollectionOfPoints>(pts);
const Vector L3 = GetTetraL3<Scalar, Vector, CollectionOfPoints>(pts);
const Vector L4 = GetTetraL4<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetTetraVolume(const CollectionOfPoints& pts)
{
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 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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetTetraInradius(const CollectionOfPoints& pts)
{
const Scalar three(3.0);
const Scalar volume = GetTetraVolume<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar area = GetTetraArea<Scalar, Vector, CollectionOfPoints>(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 <typename Scalar, typename Vector, typename CollectionOfPoints>
VTKM_EXEC Scalar GetTetraCircumradius(const CollectionOfPoints& pts)
{
const Vector L0 = GetTetraL0<Scalar, Vector, CollectionOfPoints>(pts);
const Vector L1 = GetTetraL1<Scalar, Vector, CollectionOfPoints>(pts);
const Vector L2 = GetTetraL2<Scalar, Vector, CollectionOfPoints>(pts);
const Vector L3 = GetTetraL3<Scalar, Vector, CollectionOfPoints>(pts);
const Vector L4 = GetTetraL4<Scalar, Vector, CollectionOfPoints>(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<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar circumradius = d / (twelve * volume);
return circumradius;
}
#endif

192
vtkm/exec/cellmetrics/TypeOfCellTriangle.h Normal file
View File

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

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

@ -31,8 +31,8 @@
#include <vtkm/io/writer/VTKDataSetWriter.h>
//Adapted from vtkm/cont/testing/MakeTestDataSet.h
//Modified the content of the Make3DExplicitDataSetZoo() function
inline vtkm::cont::DataSet Make3DExplicitDataSet()
//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<CoordType> coords = {
{ 0.00, 0.00, 0.00 }, { 1.00, 0.00, 0.00 }, { 2.00, 0.00, 0.00 }, { 0.00, 0.00, 1.00 },
{ 1.00, 0.00, 1.00 }, { 2.00, 0.00, 1.00 }, { 0.00, 1.00, 0.00 }, { 1.00, 1.00, 0.00 },
{ 2.00, 1.00, 0.00 }, { 0.00, 1.00, 1.00 }, { 1.00, 1.00, 1.00 }, { 2.00, 1.00, 1.00 },
{ 0.00, 2.00, 0.00 }, { 1.00, 2.00, 0.00 }, { 2.00, 2.00, 0.00 }, { 0.00, 2.00, 1.00 },
{ 1.00, 2.00, 1.00 }, { 2.00, 2.00, 1.00 }, { 1.00, 3.00, 1.00 }, { 2.75, 0.00, 1.00 },
{ 3.00, 0.00, 0.75 }, { 3.00, 0.25, 1.00 }, { 3.00, 1.00, 1.00 }, { 3.00, 1.00, 0.00 },
{ 2.57, 2.00, 1.00 }, { 3.00, 1.75, 1.00 }, { 3.00, 1.75, 0.75 }, { 3.00, 0.00, 0.00 },
{ 2.57, 0.42, 0.57 }, { 2.59, 1.43, 0.71 }
{ 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<vtkm::UInt8> 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 <typename T>
std::vector<std::string> TestMeshQualityFilter(const vtkm::cont::DataSet& input,
const std::vector<vtkm::FloatDefault>& expectedVals,
const std::string& outputname,
T filter)
bool TestMeshQualityFilter(const vtkm::cont::DataSet& input,
const std::vector<vtkm::FloatDefault>& expectedVals,
const std::string& outputname,
vtkm::filter::MeshQuality& filter)
{
std::vector<std::string> errors;
vtkm::cont::DataSet output;
try
{
@ -202,168 +118,93 @@ std::vector<std::string> 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<std::string>& messages)
{
if (!messages.empty())
vtkm::cont::ArrayHandle<vtkm::FloatDefault> 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<std::string> 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<vtkm::FloatDefault>;
using PairVec = std::vector<vtkm::Pair<vtkm::UInt8, vtkm::filter::CellMetric>>;
using StringVec = std::vector<std::string>;
using CharVec = std::vector<vtkm::UInt8>;
using QualityFilter = vtkm::filter::MeshQuality;
//Test variables
vtkm::cont::DataSet input = Make3DExplicitDataSet();
std::unique_ptr<QualityFilter> filter;
std::string metricSuffix;
StringVec fieldNames;
CharVec testedShapes;
PairVec shapeMetricPairs;
FloatVec expectedValues;
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<FloatVec> expectedValues;
std::vector<vtkm::filter::CellMetric> metrics;
std::vector<std::string> 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<std::string> errors =
TestMeshQualityFilter<QualityFilter>(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<QualityFilter>(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<QualityFilter>(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<QualityFilter>(input, expectedValues, *filter);
std::cout << "Oddy metric test: ";
CheckForErrors(errors);
/***************************************************
* Test 4: Relative Size Squared metric
***************************************************/
std::cout << "Testing MeshQuality filter: Relative Size Squared metric"
<< "\n++++++++++++++++++++++++++++++++++++++++++++++++++\n";
testedShapes = {vtkm::CELL_SHAPE_HEXAHEDRON, vtkm::CELL_SHAPE_POLYGON,
vtkm::CELL_SHAPE_QUAD, vtkm::CELL_SHAPE_TRIANGLE,
vtkm::CELL_SHAPE_TETRA};
shapeMetricPairs.clear();
for (auto s : testedShapes)
shapeMetricPairs.push_back(vtkm::make_Pair(s, vtkm::filter::CellMetric::RELATIVE_SIZE));
expectedValues = {1, 1, 0.0341086, 0.303456, 0, 0, 0,
0, 0, 0, 0, 0, 0.361898, 0.614047, 1, 1,
0.307024, 1};
filter.reset(new QualityFilter(shapeMetricPairs));
errors = TestMeshQualityFilter<QualityFilter>(input, expectedValues, *filter);
std::cout << "Relative Size Square metric test: ";
CheckForErrors(errors);
#endif
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;
}