Shuenhoy/vtk-m
1
0
Fork 0
mirror of https://gitlab.kitware.com/vtk/vtk-m synced 2024-09-19 10:35:42 +00:00
Code Issues 2 Actions Packages Projects Releases Wiki Activity

finish copying over changes from failed firs git

Browse Source
This commit is contained in:
Hank Childs 2019-09-23 22:30:49 -07:00
parent 8af25aca48
commit fcc169fdf0
28 changed files with 3431 additions and 150 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
vtkm/filter/MeshQuality.h
View File

@ -34,9 +34,10 @@ namespace filter
//Names of the available cell metrics, for use in
//the output dataset fields
static const std::string MetricNames[] = {
"area", "aspectGamma", "aspectRatio", "condition", "diagonalRatio", "jacobian",
"minAngle", "maxAngle", "oddy", "relativeSize", "scaledJacobian", "shape",
"shear", "skew", "stretch", "taper", "volume", "warpage"
"area", "aspectGamma", "aspectRatio", "condition", "diagonalRatio", "dimension",
"jacobian", "maxAngle", "maxDiagonal", "minAngle", "minDiagonal", "oddy",
"relativeSize", "scaledJacobian", "shape", "shapeAndSize", "shear", "skew",
"stretch", "taper", "volume", "warpage"
};
//Different cell metrics available to use
@ -48,13 +49,17 @@ enum class CellMetric
ASPECT_RATIO,
CONDITION,
DIAGONAL_RATIO,
DIMENSION,
JACOBIAN,
MIN_ANGLE,
MAX_ANGLE,
MAX_DIAGONAL,
MIN_ANGLE,
MIN_DIAGONAL,
ODDY,
RELATIVE_SIZE,
SCALED_JACOBIAN,
SHAPE,
SHAPE_AND_SIZE,
SHEAR,
SKEW,
STRETCH,

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

@ -138,7 +138,9 @@ bool TestMeshQualityFilter(const vtkm::cont::DataSet& input,
for (unsigned long i = 0; i < expectedVals.size(); i++)
{
vtkm::Id id = (vtkm::Id)i;
if (portal1.Get(id) != expectedVals[i])
const vtkm::FloatDefault tol = (vtkm::FloatDefault)1.e-4;
const vtkm::FloatDefault diff = portal1.Get(id) - expectedVals[i];
if (vtkm::Abs(diff) > tol)
{
printf("Metric \"%s\" for cell type \"%s\" does not match. Expected %f and got %f\n",
outputname.c_str(),
@ -165,23 +167,121 @@ int TestMeshQuality()
std::vector<vtkm::filter::CellMetric> metrics;
std::vector<std::string> metricName;
/*
FloatVec volumeExpectedValues = { 0, 0, 1, (float) 1.333333333, 4, 4 };
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");
*/
FloatVec areaExpectedValues = { 3, 4, 0, 0, 0, 0 };
expectedValues.push_back(areaExpectedValues);
metrics.push_back(vtkm::filter::CellMetric::AREA);
metricName.push_back("area");
FloatVec aspectRatioExpectedValues = { 0, (float)1.118034, 0, 0, 0, (float)1.1547 };
expectedValues.push_back(aspectRatioExpectedValues);
metrics.push_back(vtkm::filter::CellMetric::ASPECT_RATIO);
metricName.push_back("aspectRatio");
FloatVec aspectGammaExpectedValues = { 0, 0, (float)1.52012, 0, 0, 0 };
expectedValues.push_back(aspectGammaExpectedValues);
metrics.push_back(vtkm::filter::CellMetric::ASPECT_GAMMA);
metricName.push_back("aspectGamma");
FloatVec conditionExpectedValues = {
(float)1.058475, 2.25, (float)1.354007, 0, 0, (float)1.563472
};
expectedValues.push_back(conditionExpectedValues);
metrics.push_back(vtkm::filter::CellMetric::CONDITION);
metricName.push_back("condition");
FloatVec minAngleExpectedValues = { 45, 45, -1, -1, -1, -1 };
expectedValues.push_back(minAngleExpectedValues);
metrics.push_back(vtkm::filter::CellMetric::MIN_ANGLE);
metricName.push_back("minAngle");
FloatVec maxAngleExpectedValues = { (float)71.56505, 135, -1, -1, -1, -1 };
expectedValues.push_back(maxAngleExpectedValues);
metrics.push_back(vtkm::filter::CellMetric::MAX_ANGLE);
metricName.push_back("maxAngle");
FloatVec minDiagonalExpectedValues = { -1, -1, -1, -1, -1, (float)1.73205 };
expectedValues.push_back(minDiagonalExpectedValues);
metrics.push_back(vtkm::filter::CellMetric::MIN_DIAGONAL);
metricName.push_back("minDiagonal");
FloatVec maxDiagonalExpectedValues = { -1, -1, -1, -1, -1, (float)4.3589 };
expectedValues.push_back(maxDiagonalExpectedValues);
metrics.push_back(vtkm::filter::CellMetric::MAX_DIAGONAL);
metricName.push_back("maxDiagonal");
FloatVec jacobianExpectedValues = { 0, 2, 6, 0, 0, 4 };
expectedValues.push_back(jacobianExpectedValues);
metrics.push_back(vtkm::filter::CellMetric::JACOBIAN);
metricName.push_back("jacobian");
/*
FloatVec diagonalRatioExpectedValues = { -1, -1, -1, -1, -1, (float) 0.39 };
FloatVec scaledJacobianExpectedValues = {
(float)0.816497, (float)0.707107, (float)0.408248, -2, -2, (float)0.57735
};
expectedValues.push_back(scaledJacobianExpectedValues);
metrics.push_back(vtkm::filter::CellMetric::SCALED_JACOBIAN);
metricName.push_back("scaledJacobian");
FloatVec oddyExpectedValues = { -1, 8.125, -1, -1, -1, (float)2.62484 };
expectedValues.push_back(oddyExpectedValues);
metrics.push_back(vtkm::filter::CellMetric::ODDY);
metricName.push_back("oddy");
FloatVec diagonalRatioExpectedValues = { -1, (float)0.620174, -1, -1, -1, (float)0.397360 };
expectedValues.push_back(diagonalRatioExpectedValues);
metrics.push_back(vtkm::filter::CellMetric::DIAGONAL_RATIO);
metricName.push_back("diagonalRatio");
FloatVec shapeExpectedValues = { (float)0.944755, (float)0.444444, (float)0.756394, -1, -1,
(float)0.68723 };
expectedValues.push_back(shapeExpectedValues);
metrics.push_back(vtkm::filter::CellMetric::SHAPE);
metricName.push_back("shape");
FloatVec shearExpectedValues = { (float)-1, (float).707107, -1, -1, -1, (float)0.57735 };
expectedValues.push_back(shearExpectedValues);
metrics.push_back(vtkm::filter::CellMetric::SHEAR);
metricName.push_back("shear");
FloatVec skewExpectedValues = { (float)-1, (float)0.447214, -1, -1, -1, (float)0.57735 };
expectedValues.push_back(skewExpectedValues);
metrics.push_back(vtkm::filter::CellMetric::SKEW);
metricName.push_back("skew");
FloatVec stretchExpectedValues = { -1, (float)0.392232, -1, -1, -1, (float)0.688247 };
expectedValues.push_back(stretchExpectedValues);
metrics.push_back(vtkm::filter::CellMetric::STRETCH);
metricName.push_back("stretch");
FloatVec taperExpectedValues = { -1, 0.5, -1, -1, -1, 0 };
expectedValues.push_back(taperExpectedValues);
metrics.push_back(vtkm::filter::CellMetric::TAPER);
metricName.push_back("taper");
FloatVec warpageExpectedValues = { -1, 1, -1, -1, -1, -1 };
expectedValues.push_back(warpageExpectedValues);
metrics.push_back(vtkm::filter::CellMetric::WARPAGE);
metricName.push_back("warpage");
/*
FloatVec dimensionExpectedValues = { -1, -1, -1, -1, -1, (float) 0.707107 };
expectedValues.push_back(dimensionExpectedValues);
metrics.push_back(vtkm::filter::CellMetric::DIMENSION);
metricName.push_back("dimension");
FloatVec relSizeExpectedValues = { 1, 1, 1, -1, -1, 1 };
expectedValues.push_back(relSizeExpectedValues);
metrics.push_back(vtkm::filter::CellMetric::RELATIVE_SIZE);
metricName.push_back("relativeSize");
FloatVec shapeAndSizeExpectedValues = { (float) 0.944755, (float) 0.444444, (float) 0.756394, -1, -1,(float) 0.68723 };
expectedValues.push_back(shapeAndSizeExpectedValues);
metrics.push_back(vtkm::filter::CellMetric::SHAPE_AND_SIZE);
metricName.push_back("shapeAndSize");
*/
unsigned long numTests = (unsigned long)metrics.size();

103
vtkm/worklet/MeshQuality.h
View File

@ -21,11 +21,25 @@
#ifndef vtk_m_worklet_MeshQuality_h
#define vtk_m_worklet_MeshQuality_h
#include "vtkm/exec/CellMeasure.h"
#include "vtkm/exec/cellmetrics/CellDiagonalRatioMetric.h"
#include "vtkm/exec/cellmetrics/CellEdgeRatioMetric.h"
#include "vtkm/exec/cellmetrics/CellJacobianMetric.h"
#include "vtkm/worklet/CellMeasure.h"
#include "vtkm/worklet/WorkletMapTopology.h"
#include "vtkm/worklet/cellmetrics/CellAspectGammaMetric.h"
#include "vtkm/worklet/cellmetrics/CellAspectRatioMetric.h"
#include "vtkm/worklet/cellmetrics/CellConditionMetric.h"
#include "vtkm/worklet/cellmetrics/CellDiagonalRatioMetric.h"
#include "vtkm/worklet/cellmetrics/CellJacobianMetric.h"
#include "vtkm/worklet/cellmetrics/CellMaxAngleMetric.h"
#include "vtkm/worklet/cellmetrics/CellMaxDiagonalMetric.h"
#include "vtkm/worklet/cellmetrics/CellMinAngleMetric.h"
#include "vtkm/worklet/cellmetrics/CellMinDiagonalMetric.h"
#include "vtkm/worklet/cellmetrics/CellOddyMetric.h"
#include "vtkm/worklet/cellmetrics/CellScaledJacobianMetric.h"
#include "vtkm/worklet/cellmetrics/CellShapeMetric.h"
#include "vtkm/worklet/cellmetrics/CellShearMetric.h"
#include "vtkm/worklet/cellmetrics/CellSkewMetric.h"
#include "vtkm/worklet/cellmetrics/CellStretchMetric.h"
#include "vtkm/worklet/cellmetrics/CellTaperMetric.h"
#include "vtkm/worklet/cellmetrics/CellWarpageMetric.h"
namespace vtkm
{
@ -94,23 +108,88 @@ protected:
{
switch (metric)
{
case MetricTagType::AREA:
metricValue = vtkm::exec::CellMeasure<OutType>(numPts, pts, tag, *this);
if (dims != 2)
metricValue = 0.;
break;
case MetricTagType::ASPECT_GAMMA:
metricValue =
vtkm::worklet::cellmetrics::CellAspectGammaMetric<OutType>(numPts, pts, tag, *this);
break;
case MetricTagType::ASPECT_RATIO:
metricValue =
vtkm::worklet::cellmetrics::CellAspectRatioMetric<OutType>(numPts, pts, tag, *this);
break;
case MetricTagType::CONDITION:
metricValue =
vtkm::worklet::cellmetrics::CellConditionMetric<OutType>(numPts, pts, tag, *this);
break;
case MetricTagType::DIAGONAL_RATIO:
metricValue =
vtkm::exec::cellmetrics::CellDiagonalRatioMetric<OutType>(numPts, pts, tag, *this);
vtkm::worklet::cellmetrics::CellDiagonalRatioMetric<OutType>(numPts, pts, tag, *this);
break;
/* DEPRECATING
case MetricTagType::EDGE_RATIO:
metricValue =
vtkm::exec::cellmetrics::CellEdgeRatioMetric<OutType>(numPts, pts, tag, *this);
break;
*/
case MetricTagType::JACOBIAN:
metricValue =
vtkm::exec::cellmetrics::CellJacobianMetric<OutType>(numPts, pts, tag, *this);
vtkm::worklet::cellmetrics::CellJacobianMetric<OutType>(numPts, pts, tag, *this);
break;
case MetricTagType::MAX_ANGLE:
metricValue =
vtkm::worklet::cellmetrics::CellMaxAngleMetric<OutType>(numPts, pts, tag, *this);
break;
case MetricTagType::MAX_DIAGONAL:
metricValue =
vtkm::worklet::cellmetrics::CellMaxDiagonalMetric<OutType>(numPts, pts, tag, *this);
break;
case MetricTagType::MIN_ANGLE:
metricValue =
vtkm::worklet::cellmetrics::CellMinAngleMetric<OutType>(numPts, pts, tag, *this);
break;
case MetricTagType::MIN_DIAGONAL:
metricValue =
vtkm::worklet::cellmetrics::CellMinDiagonalMetric<OutType>(numPts, pts, tag, *this);
break;
case MetricTagType::ODDY:
metricValue =
vtkm::worklet::cellmetrics::CellOddyMetric<OutType>(numPts, pts, tag, *this);
break;
case MetricTagType::SCALED_JACOBIAN:
metricValue =
vtkm::worklet::cellmetrics::CellScaledJacobianMetric<OutType>(numPts, pts, tag, *this);
break;
case MetricTagType::SHAPE:
metricValue =
vtkm::worklet::cellmetrics::CellShapeMetric<OutType>(numPts, pts, tag, *this);
break;
case MetricTagType::SHEAR:
metricValue =
vtkm::worklet::cellmetrics::CellShearMetric<OutType>(numPts, pts, tag, *this);
break;
case MetricTagType::SKEW:
metricValue =
vtkm::worklet::cellmetrics::CellSkewMetric<OutType>(numPts, pts, tag, *this);
break;
case MetricTagType::STRETCH:
metricValue =
vtkm::worklet::cellmetrics::CellStretchMetric<OutType>(numPts, pts, tag, *this);
break;
case MetricTagType::TAPER:
metricValue =
vtkm::worklet::cellmetrics::CellTaperMetric<OutType>(numPts, pts, tag, *this);
break;
case MetricTagType::VOLUME:
metricValue = vtkm::exec::CellMeasure<OutType>(numPts, pts, tag, *this);
if (dims != 3)
metricValue = 0.;
break;
case MetricTagType::WARPAGE:
metricValue =
vtkm::worklet::cellmetrics::CellWarpageMetric<OutType>(numPts, pts, tag, *this);
break;
case MetricTagType::DIMENSION:
case MetricTagType::RELATIVE_SIZE:
case MetricTagType::SHAPE_AND_SIZE:
this->RaiseError("Asked for unimplemented metric.");
case MetricTagType::EMPTY:
break;
default:

17
vtkm/worklet/cellmetrics/CMakeLists.txt
View File

@ -19,9 +19,26 @@
##============================================================================
set(headers
CellAspectFrobeniusMetric.h
CellAspectGammaMetric.h
CellAspectRatioMetric.h
CellConditionMetric.h
CellDiagonalRatioMetric.h
CellEdgeRatioMetric.h
CellJacobianMetric.h
CellMaxAngleMetric.h
CellMaxAspectFrobeniusMetric.h
CellMaxDiagonalMetric.h
CellMinAngleMetric.h
CellMinDiagonalMetric.h
CellOddyMetric.h
CellScaledJacobianMetric.h
CellShapeMetric.h
CellShearMetric.h
CellSkewMetric.h
CellStretchMetric.h
CellTaperMetric.h
CellWarpageMetric.h
TypeOfCellHexahedral.h
TypeOfCellQuadrilateral.h
TypeOfCellTetrahedral.h

244
vtkm/worklet/cellmetrics/CellAspectFrobeniusMetric.h Normal file
View File

@ -0,0 +1,244 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//
// Copyright 2014 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
// Copyright 2014 UT-Battelle, LLC.
// Copyright 2014 Los Alamos National Security.
//
// Under the terms of Contract DE-NA0003525 with NTESS,
// the U.S. Government retains certain rights in this software.
//
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#ifndef vtk_m_worklet_cellmetrics_CellAspectFrobeniusMetric_h
#define vtk_m_worklet_cellmetrics_CellAspectFrobeniusMetric_h
/*
* Mesh quality metric functions that compute the aspect frobenius of certain mesh cells.
* The aspect frobenius metric generally measures the degree of regularity of a cell, with
* a value of 1 representing a regular cell..
*
* These metric computations are adapted from the VTK implementation of the Verdict library,
* which provides a set of mesh/cell metrics for evaluating the geometric qualities of regions
* of mesh spaces.
*
* See: The Verdict Library Reference Manual (for per-cell-type metric formulae)
* See: vtk/ThirdParty/verdict/vtkverdict (for VTK code implementation of this metric)
*/
#include "vtkm/CellShape.h"
#include "vtkm/CellTraits.h"
#include "vtkm/VecTraits.h"
#include "vtkm/VectorAnalysis.h"
#include "vtkm/exec/FunctorBase.h"
#define UNUSED(expr) (void)(expr);
namespace vtkm
{
namespace worklet
{
namespace cellmetrics
{
using FloatType = vtkm::FloatDefault;
// ========================= Unsupported cells ==================================
// By default, cells have undefined aspect frobenius unless the shape type template is specialized below.
template <typename OutType, typename PointCoordVecType, typename CellShapeType>
VTKM_EXEC OutType CellAspectFrobeniusMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
CellShapeType shape,
const vtkm::exec::FunctorBase& worklet)
{
UNUSED(numPts);
UNUSED(pts);
UNUSED(shape);
worklet.RaiseError(
"Shape type template must be specialized to compute the aspect frobenius metric.");
return OutType(0.0);
} //If the polygon has 3 vertices or 4 vertices, then just call
//the functions for Triangle and Quad cell types. Otherwise,
//this metric is not supported for (n>4)-vertex polygons, such
//as pentagons or hexagons, or (n<3)-vertex polygons, such as lines or points.
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellAspectFrobeniusMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagPolygon,
const vtkm::exec::FunctorBase& worklet)
{
if (numPts == 3)
return CellAspectFrobeniusMetric<OutType>(numPts, pts, vtkm::CellShapeTagTriangle(), worklet);
else
{
worklet.RaiseError(
"Aspect frobenius metric is not supported for (n<3)- or (n>4)-vertex polygons.");
return OutType(0.0);
}
} //The aspect frobenius metric is not supported for lines/edges.
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellAspectFrobeniusMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagLine,
const vtkm::exec::FunctorBase& worklet)
{
UNUSED(numPts);
UNUSED(pts);
worklet.RaiseError("Aspect frobenius metric is not supported for lines/edges.");
return OutType(0.0);
} //The aspect frobenius metric is not uniquely defined for quads.
//Instead, use either the mean or max aspect frobenius metrics, which are
//defined in terms of the aspect frobenius of triangles.
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellAspectFrobeniusMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagQuad,
const vtkm::exec::FunctorBase& worklet)
{
UNUSED(numPts);
UNUSED(pts);
worklet.RaiseError("Aspect frobenius metric is not supported for quads.");
return OutType(0.0);
} //The aspect frobenius metric is not uniquely defined for hexahedrons.
//Instead, use either the mean or max aspect frobenius metrics, which are
//defined in terms of the aspect frobenius of tetrahedrons.
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellAspectFrobeniusMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagHexahedron,
const vtkm::exec::FunctorBase& worklet)
{
UNUSED(numPts);
UNUSED(pts);
worklet.RaiseError("Aspect frobenius metric is not supported for hexahedrons.");
return OutType(0.0);
} //The aspect frobenius metric is not uniquely defined for pyramids.
//Instead, use either the mean or max aspect frobenius metrics, which are
//defined in terms of the aspect frobenius of tetrahedrons.
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellAspectFrobeniusMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagPyramid,
const vtkm::exec::FunctorBase& worklet)
{
UNUSED(numPts);
UNUSED(pts);
worklet.RaiseError("Aspect frobenius metric is not supported for pyramids.");
return OutType(0.0);
} //The aspect frobenius metric is not uniquely defined for wedges.
//Instead, use either the mean or max aspect frobenius metrics, which are
//defined in terms of the aspect frobenius of tetrahedrons.
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellAspectFrobeniusMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagWedge,
const vtkm::exec::FunctorBase& worklet)
{
UNUSED(numPts);
UNUSED(pts);
worklet.RaiseError("Aspect frobenius metric is not supported for wedges.");
return OutType(0.0);
} // ========================= 2D cells ==================================
// Computes the aspect frobenius of a triangle.
// Formula: Sum of lengths of 3 edges, divided by a multiple of the triangle area.
// Equals 1 for an equilateral unit triangle.
// Acceptable range: [1,1.3]
// Full range: [1,FLOAT_MAX]
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellAspectFrobeniusMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagTriangle,
const vtkm::exec::FunctorBase& worklet)
{
if (numPts != 3)
{
worklet.RaiseError("Aspect frobenius metric(triangle) requires 3 points.");
return OutType(0.0);
}
//The 3 edges of a triangle
using Edge = typename PointCoordVecType::ComponentType;
const Edge TriEdges[3] = { pts[1] - pts[0], pts[2] - pts[1], pts[0] - pts[2] };
//Sum the length squared of each edge
FloatType sum = (FloatType)vtkm::MagnitudeSquared(TriEdges[0]) +
(FloatType)vtkm::MagnitudeSquared(TriEdges[1]) + (FloatType)vtkm::MagnitudeSquared(TriEdges[2]);
//Compute the length of the cross product of the triangle.
//The result is twice the area of the triangle.
FloatType crossLen = (FloatType)vtkm::Magnitude(vtkm::Cross(TriEdges[0], -TriEdges[2]));
if (crossLen == 0.0)
return vtkm::Infinity<OutType>();
OutType aspect_frobenius = (OutType)(sum / (vtkm::Sqrt(3.0) * 2 * crossLen));
if (aspect_frobenius > 0.0)
return vtkm::Min(aspect_frobenius, vtkm::Infinity<OutType>());
return vtkm::Max(aspect_frobenius, vtkm::NegativeInfinity<OutType>());
} // ============================= 3D Volume cells ==================================i
// Computes the aspect frobenius of a tetrahedron.
// Formula: Sum of lengths of 3 edges, divided by a multiple of the triangle area.
// Equals 1 for a right regular tetrahedron (4 equilateral triangles).
// Acceptable range: [1,1.3]
// Full range: [1,FLOAT_MAX]
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellAspectFrobeniusMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagTetra,
const vtkm::exec::FunctorBase& worklet)
{
if (numPts != 4)
{
worklet.RaiseError("Aspect frobenius metric(tetrahedron) requires 4 points.");
return OutType(0.0);
}
//Two base edges and one vertical edge, used to compute the tet volume
using Edge = typename PointCoordVecType::ComponentType;
const Edge TetEdges[3] = {
pts[1] - pts[0], //Base edge 1
pts[2] - pts[0], //Base edge 2
pts[3] - pts[0] //Vert edge 3
};
//Compute the tet volume
FloatType denominator = (FloatType)vtkm::Dot(TetEdges[0], vtkm::Cross(TetEdges[1], TetEdges[2]));
denominator *= denominator;
denominator *= 2.0f;
const FloatType normal_exp = 1.0f / 3.0f;
denominator = 3.0f * vtkm::Pow(denominator, normal_exp);
if (denominator < vtkm::NegativeInfinity<FloatType>())
return vtkm::Infinity<OutType>();
FloatType numerator = (FloatType)vtkm::Dot(TetEdges[0], TetEdges[0]);
numerator += (FloatType)vtkm::Dot(TetEdges[1], TetEdges[1]);
numerator += (FloatType)vtkm::Dot(TetEdges[2], TetEdges[2]);
numerator *= 1.5;
numerator -= (FloatType)vtkm::Dot(TetEdges[0], TetEdges[1]);
numerator -= (FloatType)vtkm::Dot(TetEdges[0], TetEdges[2]);
numerator -= (FloatType)vtkm::Dot(TetEdges[1], TetEdges[2]);
OutType aspect_frobenius = (OutType)(numerator / denominator);
if (aspect_frobenius > 0.0)
return vtkm::Min(aspect_frobenius, vtkm::Infinity<OutType>());
return vtkm::Max(aspect_frobenius, vtkm::NegativeInfinity<OutType>());
}
} // namespace cellmetrics
} // namespace worklet
} // namespace vtkm
#endif // vtk_m_worklet_cellmetrics_CellAspectFrobeniusMetric_h

111
vtkm/worklet/cellmetrics/CellAspectGammaMetric.h Normal file
View File

@ -0,0 +1,111 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//
// Copyright 2018 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
// Copyright 2018 UT-Battelle, LLC.
// Copyright 2018 Los Alamos National Security.
//
// Under the terms of Contract DE-NA0003525 with NTESS,
// the U.S. Government retains certain rights in this software.
//
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#ifndef vtk_m_worklet_cellmetrics_CellAspectGammaMetric_h
#define vtk_m_worklet_cellmetrics_CellAspectGammaMetric_h
/*
* Mesh quality metric functions that compute the aspect ratio of mesh cells.
** These metric computations are adapted from the VTK implementation of the Verdict library,
* which provides a set of mesh/cell metrics for evaluating the geometric qualities of regions
* of mesh spaces.
** See: The Verdict Library Reference Manual (for per-cell-type metric formulae)
* See: vtk/ThirdParty/verdict/vtkverdict (for VTK code implementation of this metric)
*/
#include "TypeOfCellHexahedral.h"
#include "TypeOfCellQuadrilateral.h"
#include "TypeOfCellTetrahedral.h"
#include "TypeOfCellTriangle.h"
#include "vtkm/CellShape.h"
#include "vtkm/CellTraits.h"
#include "vtkm/VecTraits.h"
#include "vtkm/VectorAnalysis.h"
#include "vtkm/exec/FunctorBase.h"
#define UNUSED(expr) (void)(expr);
namespace vtkm
{
namespace worklet
{
namespace cellmetrics
{
// ========================= Unsupported cells ==================================
// By default, cells have zero shape unless the shape type template is specialized below.
template <typename OutType, typename PointCoordVecType, typename CellShapeType>
VTKM_EXEC OutType CellAspectGammaMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
CellShapeType shape,
const vtkm::exec::FunctorBase&)
{
UNUSED(numPts);
UNUSED(pts);
UNUSED(shape);
return OutType(0);
}
// ============================= 3D Volume cells ==================================
// Compute the aspect ratio of a tetrahedron.
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellAspectGammaMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagTetra,
const vtkm::exec::FunctorBase& worklet)
{
if (numPts != 4)
{
worklet.RaiseError("Aspect gamma metric (tetrahedron) requires 4 points.");
return OutType(0.0);
}
using Scalar = OutType;
using CollectionOfPoints = PointCoordVecType;
using Vector = typename PointCoordVecType::ComponentType;
const Scalar volume = GetTetraVolume<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar vAbs = vtkm::Abs(volume);
if (vAbs <= Scalar(0.0))
{
return vtkm::Infinity<Scalar>();
}
const Scalar six = Scalar(6.0);
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 R = vtkm::Sqrt((vtkm::Pow(l0, 2) + vtkm::Pow(l1, 2) + vtkm::Pow(l2, 2) +
vtkm::Pow(l3, 2) + vtkm::Pow(l4, 2) + vtkm::Pow(l5, 2)) /
six);
const Scalar rootTwo(vtkm::Sqrt(Scalar(2.0)));
const Scalar twelve(12.0);
const Scalar q = (vtkm::Pow(R, 3) * rootTwo) / (twelve * vAbs);
return q;
}
} // namespace cellmetrics
} // namespace worklet
} // namespace vtkm
#endif // vtk_m_worklet_cellmetrics_CellAspectGammaMetric_h

221
vtkm/worklet/cellmetrics/CellAspectRatioMetric.h Normal file
View File

@ -0,0 +1,221 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//
// Copyright 2018 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
// Copyright 2018 UT-Battelle, LLC.
// Copyright 2018 Los Alamos National Security.
//
// Under the terms of Contract DE-NA0003525 with NTESS,
// the U.S. Government retains certain rights in this software.
//
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#ifndef vtk_m_worklet_cellmetrics_CellAspectRatioMetric_h
#define vtk_m_worklet_cellmetrics_CellAspectRatioMetric_h
/*
* Mesh quality metric functions that compute the aspect ratio of mesh cells.
** These metric computations are adapted from the VTK implementation of the Verdict library,
* which provides a set of mesh/cell metrics for evaluating the geometric qualities of regions
* of mesh spaces.
** See: The Verdict Library Reference Manual (for per-cell-type metric formulae)
* See: vtk/ThirdParty/verdict/vtkverdict (for VTK code implementation of this metric)
*/
#include "TypeOfCellHexahedral.h"
#include "TypeOfCellQuadrilateral.h"
#include "TypeOfCellTetrahedral.h"
#include "TypeOfCellTriangle.h"
#include "vtkm/CellShape.h"
#include "vtkm/CellTraits.h"
#include "vtkm/VecTraits.h"
#include "vtkm/VectorAnalysis.h"
#include "vtkm/exec/FunctorBase.h"
#define UNUSED(expr) (void)(expr);
namespace vtkm
{
namespace worklet
{
namespace cellmetrics
{
// The Verdict Manual and the Implementation have conflicting definitions.
// This duplicates the Verdict implementation in the VTKm Paradigm, with prior Manual
// definitions commented out when formerly coded.
// ========================= Unsupported cells ==================================
// By default, cells have zero shape unless the shape type template is specialized below.
template <typename OutType, typename PointCoordVecType, typename CellShapeType>
VTKM_EXEC OutType CellAspectRatioMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
CellShapeType shape,
const vtkm::exec::FunctorBase&)
{
UNUSED(numPts);
UNUSED(pts);
UNUSED(shape);
return OutType(0);
}
// ========================= 2D cells ==================================
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellAspectRatioMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagQuad,
const vtkm::exec::FunctorBase& worklet)
{
if (numPts != 4)
{
worklet.RaiseError("Aspect ratio metric (quad) requires 4 points.");
return OutType(0.0);
}
using Scalar = OutType;
using CollectionOfPoints = PointCoordVecType;
using Vector = typename PointCoordVecType::ComponentType;
const Vector X1 = GetQuadX0<Scalar, Vector, CollectionOfPoints>(pts);
const Vector X2 = GetQuadX1<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar x1 = vtkm::Sqrt(vtkm::MagnitudeSquared(X1));
const Scalar x2 = vtkm::Sqrt(vtkm::MagnitudeSquared(X2));
if (x1 <= Scalar(0.0) || x2 <= Scalar(0.0))
{
return vtkm::Infinity<Scalar>();
}
const Scalar q = vtkm::Max(x1 / x2, x2 / x1);
return q;
}
// ========================= 3D cells ==================================
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellAspectRatioMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagHexahedron,
const vtkm::exec::FunctorBase& worklet)
{
if (numPts != 8)
{
worklet.RaiseError("Aspect ratio metric (hex) requires 4 points.");
return OutType(0.0);
}
using Scalar = OutType;
using CollectionOfPoints = PointCoordVecType;
using Vector = typename PointCoordVecType::ComponentType;
const Vector X1 = GetHexX1<Scalar, Vector, CollectionOfPoints>(pts);
const Vector X2 = GetHexX2<Scalar, Vector, CollectionOfPoints>(pts);
const Vector X3 = GetHexX3<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar x1 = vtkm::Sqrt(vtkm::MagnitudeSquared(X1));
const Scalar x2 = vtkm::Sqrt(vtkm::MagnitudeSquared(X2));
const Scalar x3 = vtkm::Sqrt(vtkm::MagnitudeSquared(X3));
if (x1 <= Scalar(0.0) || x2 <= Scalar(0.0) || x3 <= Scalar(0.0))
{
return vtkm::Infinity<Scalar>();
}
const Scalar q = vtkm::Max(
x1 / x2,
vtkm::Max(x2 / x1, vtkm::Max(x1 / x3, vtkm::Max(x3 / x1, vtkm::Max(x3 / x2, x3 / x2)))));
return q;
}
///////////////////////////////////////////////////
////Definitions as found in the Verdict Manual/////
///////////////////////////////////////////////////
// ========================= 2D cells ==================================
// Compute the diagonal ratio of a triangle.
/*
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellAspectRatioMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagTriangle,
const vtkm::exec::FunctorBase& worklet)
{
if (numPts != 3)
{
worklet.RaiseError("Aspect ratio metric (triangle) requires 3 points.");
return OutType(0.0);
}
using Scalar = OutType;
using CollectionOfPoints = PointCoordVecType;
using Vector = typename PointCoordVecType::ComponentType;
const Scalar lmax = GetTriangleLMax <Scalar, Vector, CollectionOfPoints> (pts);
const Scalar r = GetTriangleInradius <Scalar, Vector, CollectionOfPoints> (pts);
const Scalar half(0.5);
const Scalar three(3.0);
const Scalar q = (lmax * half * vtkm::RSqrt(three)) / r;
return q;
}
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellAspectRatioMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagQuad,
const vtkm::exec::FunctorBase& worklet)
{
if (numPts != 4)
{
worklet.RaiseError("Aspect ratio metric (quad) requires 4 points.");
return OutType(0.0);
}
using Scalar = OutType;
using CollectionOfPoints = PointCoordVecType;
using Vector = typename PointCoordVecType::ComponentType;
const Scalar lmax = GetQuadLMax <Scalar, Vector, 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 four(4.0);
const Scalar area = GetQuadArea <Scalar, Vector, CollectionOfPoints> (pts);
const Scalar q = (lmax * (l0 + l1 + l2 + l3)) / (four * area);
return q;
}
// ============================= 3D Volume cells ==================================
// Compute the aspect ratio of a tetrahedron.
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellAspectRatioMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagTetra,
const vtkm::exec::FunctorBase& worklet)
{
if (numPts != 4)
{
worklet.RaiseError("Aspect ratio metric (tetrahedron) requires 4 points.");
return OutType(0.0);
}
using Scalar = OutType;
using CollectionOfPoints = PointCoordVecType;
using Vector = typename PointCoordVecType::ComponentType;
const Scalar rootSixInvert = vtkm::RSqrt(Scalar(6.0));
const Scalar half(0.5);
const Scalar lmax = GetTetraLMax <Scalar, Vector, CollectionOfPoints> (pts);
const Scalar r = GetTetraInradius <Scalar, Vector, CollectionOfPoints> (pts);
const Scalar q = (half * rootSixInvert * lmax) / r;
return q;
}
*/
} // namespace cellmetrics
} // namespace worklet
} // namespace vtkm
#endif // vtk_m_worklet_cellmetrics_CellAspectRatioMetric_h

198
vtkm/worklet/cellmetrics/CellConditionMetric.h Normal file
View File

@ -0,0 +1,198 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//
// Copyright 2014 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
// Copyright 2014 UT-Battelle, LLC.
// Copyright 2014 Los Alamos National Security.
//
// Under the terms of Contract DE-NA0003525 with NTESS,
// the U.S. Government retains certain rights in this software.
//
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#ifndef vtk_m_worklet_cellmetrics_CellConditionMetric_h
#define vtk_m_worklet_cellmetrics_CellConditionMetric_h
/*
* Mesh quality metric functions that compute the condition metric of mesh cells.
** These metric computations are adapted from the VTK implementation of the Verdict library,
* which provides a set of mesh/cell metrics for evaluating the geometric qualities of regions
* of mesh spaces.
** See: The Verdict Library Reference Manual (for per-cell-type metric formulae)
* See: vtk/ThirdParty/verdict/vtkverdict (for VTK code implementation of this metric)
*/
#include "CellMaxAspectFrobeniusMetric.h"
#include "TypeOfCellHexahedral.h"
#include "TypeOfCellQuadrilateral.h"
#include "TypeOfCellTetrahedral.h"
#include "TypeOfCellTriangle.h"
#include "vtkm/CellShape.h"
#include "vtkm/CellTraits.h"
#include "vtkm/VecTraits.h"
#include "vtkm/VectorAnalysis.h"
#include "vtkm/exec/FunctorBase.h"
#define UNUSED(expr) (void)(expr);
namespace vtkm
{
namespace worklet
{
namespace cellmetrics
{
// ========================= Unsupported cells ==================================
// By default, cells have zero shape unless the shape type template is specialized below.
template <typename OutType, typename PointCoordVecType, typename CellShapeType>
VTKM_EXEC OutType CellConditionMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
CellShapeType shape,
const vtkm::exec::FunctorBase&)
{
UNUSED(numPts);
UNUSED(pts);
UNUSED(shape);
return OutType(0);
}
// =============================== Condition metric cells ==================================
// Compute the condition quality metric of a triangular cell.
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellConditionMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagTriangle,
const vtkm::exec::FunctorBase& worklet)
{
if (numPts != 3)
{
worklet.RaiseError("Condition metric(triangle) requires 3 points.");
return OutType(0.0);
}
using Scalar = OutType;
using CollectionOfPoints = PointCoordVecType;
using Vector = typename PointCoordVecType::ComponentType;
const Scalar area = GetTriangleArea<Scalar, Vector, CollectionOfPoints>(pts);
if (area == Scalar(0.0))
{
return vtkm::Infinity<Scalar>();
}
const Scalar two(2.0);
const Scalar rootThree = vtkm::Sqrt(Scalar(3.0));
const Vector L1 = GetTriangleL1<Scalar, Vector, CollectionOfPoints>(pts);
const Vector L2 = GetTriangleL2<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar q =
(vtkm::Dot(L2, L2) + vtkm::Dot(L1, L1) + vtkm::Dot(L1, L2)) / (two * area * rootThree);
return q;
}
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellConditionMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagQuad,
const vtkm::exec::FunctorBase& worklet)
{
UNUSED(numPts);
UNUSED(worklet);
using Scalar = OutType;
using CollectionOfPoints = PointCoordVecType;
using Vector = typename PointCoordVecType::ComponentType;
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);
if (a0 < vtkm::NegativeInfinity<Scalar>() || a1 < vtkm::NegativeInfinity<Scalar>() ||
a2 < vtkm::NegativeInfinity<Scalar>() || a3 < vtkm::NegativeInfinity<Scalar>())
{
return vtkm::Infinity<Scalar>();
}
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 half(0.5);
const Scalar q0 = ((l0 * l0) + (l3 * l3)) / a0;
const Scalar q1 = ((l1 * l1) + (l0 * l0)) / a1;
const Scalar q2 = ((l2 * l2) + (l1 * l1)) / a2;
const Scalar q3 = ((l3 * l3) + (l2 * l2)) / a3;
const Scalar q = half * vtkm::Max(q0, vtkm::Max(q1, vtkm::Max(q2, q3)));
return q;
}
// ============================= 3D volumetric cells ==================================
/// Compute the condition metric of a tetrahedron.
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellConditionMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagTetra,
const vtkm::exec::FunctorBase& worklet)
{
if (numPts != 4)
{
worklet.RaiseError("Condition metric(tetrahedron) requires 4 points.");
return OutType(0.0);
}
using Scalar = OutType;
using CollectionOfPoints = PointCoordVecType;
using Vector = typename PointCoordVecType::ComponentType;
const Scalar negTwo(-2.0);
const Scalar three(3.0);
const Scalar root3 = vtkm::Sqrt(three);
const Scalar root6 = vtkm::Sqrt(Scalar(6.0));
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 Vector C1 = L0;
const Vector C2 = ((negTwo * L2) - L0) / root3;
const Vector C3 = ((three * L3) + L2 - L0) / root6;
const Scalar cDet = vtkm::Dot(C1, vtkm::Cross(C2, C3));
if (cDet <= Scalar(0.0))
{
return vtkm::Infinity<Scalar>();
}
const Vector C1xC2 = vtkm::Cross(C1, C2);
const Vector C2xC3 = vtkm::Cross(C2, C3);
const Vector C1xC3 = vtkm::Cross(C1, C3);
const Scalar t1 = vtkm::Dot(C1, C1) + vtkm::Dot(C2, C2) + vtkm::Dot(C3, C3);
const Scalar t2 = vtkm::Dot(C1xC2, C1xC2) + vtkm::Dot(C2xC3, C2xC3) + vtkm::Dot(C1xC3, C1xC3);
const Scalar q = vtkm::Sqrt(t1 * t2) / (three * cDet);
return q;
}
// Condition of a hex cell is a deprecated/legacy metric which is identical to the Max Aspect Frobenius metric.
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellConditionMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagHexahedron,
const vtkm::exec::FunctorBase& worklet)
{
return CellMaxAspectFrobeniusMetric<OutType, PointCoordVecType>(
numPts, pts, vtkm::CellShapeTagHexahedron(), worklet);
}
}
}
}
#endif // vtk_m_worklet_CellConditionMetric_h

125
vtkm/worklet/cellmetrics/CellDiagonalRatioMetric.h
View File

@ -17,24 +17,21 @@
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#ifndef vtk_m_exec_cellmetrics_CellDiagonalRatioMetric_h
#define vtk_m_exec_cellmetrics_CellDiagonalRatioMetric_h
#ifndef vtk_m_worklet_cellmetrics_CellDiagonalRatioMetric_h
#define vtk_m_worklet_cellmetrics_CellDiagonalRatioMetric_h
/*
* Mesh quality metric functions that compute the diagonal ratio of mesh cells.
* The diagonal ratio of a cell is defined as the length (magnitude) of the longest
* cell diagonal length divided by the length of the shortest cell diagonal length.
*
* These metric computations are adapted from the VTK implementation of the Verdict library,
* which provides a set of mesh/cell metrics for evaluating the geometric qualities of regions
* of mesh spaces.
*
* The edge ratio computations for a pyramid cell types is not defined in the
* VTK implementation, but is provided here.
*
* See: The Verdict Library Reference Manual (for per-cell-type metric formulae)
* See: vtk/ThirdParty/verdict/vtkverdict (for VTK code implementation of this metric)
*/
* Mesh quality metric functions that compute the diagonal ratio of mesh cells.
* The diagonal ratio of a cell is defined as the length (magnitude) of the longest
* cell diagonal length divided by the length of the shortest cell diagonal length.
** These metric computations are adapted from the VTK implementation of the Verdict library,
* which provides a set of mesh/cell metrics for evaluating the geometric qualities of regions
* of mesh spaces.
** The edge ratio computations for a pyramid cell types is not defined in the
* VTK implementation, but is provided here.
** See: The Verdict Library Reference Manual (for per-cell-type metric formulae)
* See: vtk/ThirdParty/verdict/vtkverdict (for VTK code implementation of this metric)
*/
#include "vtkm/CellShape.h"
#include "vtkm/CellTraits.h"
@ -46,14 +43,13 @@
namespace vtkm
{
namespace exec
namespace worklet
{
namespace cellmetrics
{
using FloatType = vtkm::FloatDefault;
template <typename OutType, typename VecType>
VTKM_EXEC inline OutType ComputeDiagonalRatio(const VecType& diagonals)
{
@ -73,20 +69,14 @@ VTKM_EXEC inline OutType ComputeDiagonalRatio(const VecType& diagonals)
maxLen = currLen;
}
if (minLen < vtkm::NegativeInfinity<FloatType>())
if (minLen <= OutType(0.0))
return vtkm::Infinity<OutType>();
//Take square root because we only did magnitude squared before
OutType diagonalRatio = (OutType)vtkm::Sqrt(maxLen / minLen);
if (diagonalRatio > 0)
return vtkm::Min(diagonalRatio, vtkm::Infinity<OutType>()); //normal case
return vtkm::Max(diagonalRatio, OutType(-1) * vtkm::Infinity<OutType>());
OutType diagonalRatio = (OutType)vtkm::Sqrt(minLen / maxLen);
return diagonalRatio;
}
// ========================= Unsupported cells ==================================
// By default, cells have zero shape unless the shape type template is specialized below.
template <typename OutType, typename PointCoordVecType, typename CellShapeType>
VTKM_EXEC OutType CellDiagonalRatioMetric(const vtkm::IdComponent& numPts,
@ -97,79 +87,8 @@ VTKM_EXEC OutType CellDiagonalRatioMetric(const vtkm::IdComponent& numPts,
UNUSED(numPts);
UNUSED(pts);
UNUSED(shape);
return OutType(0.0);
}
/*
//TODO: Should polygons be supported? Maybe call Quad or Triangle function...
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellDiagonalRatioMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagPolygon,
const vtkm::exec::FunctorBase& worklet)
{
switch (numPts)
{
case 4:
return CellDiagonalRatioMetric<OutType>(numPts, pts, vtkm::CellShapeTagQuad(), worklet);
default:
break;
}
return OutType(-1.0);
}
*/
/*
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellDiagonalRatioMetric(const vtkm::IdComponent&,
const PointCoordVecType&,
vtkm::CellShapeTagLine,
const vtkm::exec::FunctorBase& worklet)
{
UNUSED(worklet);
return OutType(-1.0);
}
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellDiagonalRatioMetric(const vtkm::IdComponent&,
const PointCoordVecType&,
vtkm::CellShapeTagTriangle,
const vtkm::exec::FunctorBase& worklet)
{
UNUSED(worklet);
return OutType(-1.0);
}
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellDiagonalRatioMetric(const vtkm::IdComponent&,
const PointCoordVecType&,
vtkm::CellShapeTagTetra,
const vtkm::exec::FunctorBase& worklet)
{
UNUSED(worklet);
return OutType(-1.0);
}
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellDiagonalRatioMetric(const vtkm::IdComponent&,
const PointCoordVecType&,
vtkm::CellShapeTagWedge,
const vtkm::exec::FunctorBase& worklet)
{
UNUSED(worklet);
return OutType(-1.0);
}
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellDiagonalRatioMetric(const vtkm::IdComponent&,
const PointCoordVecType&,
vtkm::CellShapeTagPyramid,
const vtkm::exec::FunctorBase& worklet)
{
UNUSED(worklet);
return OutType(-1.0);
}
*/
// ========================= 2D cells ==================================
// Compute the diagonal ratio of a quadrilateral.
@ -194,7 +113,7 @@ VTKM_EXEC OutType CellDiagonalRatioMetric(const vtkm::IdComponent& numPts,
using Diagonal = typename PointCoordVecType::ComponentType;
const Diagonal QuadDiagonals[2] = { pts[2] - pts[0], pts[3] - pts[1] };
return vtkm::exec::cellmetrics::ComputeDiagonalRatio<OutType>(
return vtkm::worklet::cellmetrics::ComputeDiagonalRatio<OutType>(
vtkm::make_VecC(QuadDiagonals, numDiagonals));
}
@ -225,12 +144,10 @@ VTKM_EXEC OutType CellDiagonalRatioMetric(const vtkm::IdComponent& numPts,
pts[6] - pts[0], pts[7] - pts[1], pts[4] - pts[2], pts[5] - pts[3]
};
return vtkm::exec::cellmetrics::ComputeDiagonalRatio<OutType>(
return vtkm::worklet::cellmetrics::ComputeDiagonalRatio<OutType>(
vtkm::make_VecC(HexDiagonals, numDiagonals));
}
} // namespace cellmetrics
} // namespace exec
} // namespace worklet
} // namespace vtkm
#endif // vtk_m_exec_cellmetrics_CellEdgeRatioMetric_h
#endif // vtk_m_worklet_cellmetrics_CellEdgeRatioMetric_h

24
vtkm/worklet/cellmetrics/CellEdgeRatioMetric.h
View File

@ -17,8 +17,8 @@
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#ifndef vtk_m_exec_cellmetrics_CellEdgeRatioMetric_h
#define vtk_m_exec_cellmetrics_CellEdgeRatioMetric_h
#ifndef vtk_m_worklet_cellmetrics_CellEdgeRatioMetric_h
#define vtk_m_worklet_cellmetrics_CellEdgeRatioMetric_h
/*
* Mesh quality metric functions that compute the edge ratio of mesh cells.
@ -46,7 +46,7 @@
namespace vtkm
{
namespace exec
namespace worklet
{
namespace cellmetrics
{
@ -143,7 +143,7 @@ VTKM_EXEC OutType CellEdgeRatioMetric(const vtkm::IdComponent& numPts,
//The 3 edges of a triangle
using Edge = typename PointCoordVecType::ComponentType;
const Edge TriEdges[3] = { pts[1] - pts[0], pts[2] - pts[1], pts[0] - pts[2] };
return vtkm::exec::cellmetrics::ComputeEdgeRatio<OutType>(vtkm::make_VecC(TriEdges, numEdges));
return vtkm::worklet::cellmetrics::ComputeEdgeRatio<OutType>(vtkm::make_VecC(TriEdges, numEdges));
}
@ -170,7 +170,8 @@ VTKM_EXEC OutType CellEdgeRatioMetric(const vtkm::IdComponent& numPts,
using Edge = typename PointCoordVecType::ComponentType;
const Edge QuadEdges[4] = { pts[1] - pts[0], pts[2] - pts[1], pts[3] - pts[2], pts[0] - pts[3] };
return vtkm::exec::cellmetrics::ComputeEdgeRatio<OutType>(vtkm::make_VecC(QuadEdges, numEdges));
return vtkm::worklet::cellmetrics::ComputeEdgeRatio<OutType>(
vtkm::make_VecC(QuadEdges, numEdges));
}
@ -201,7 +202,7 @@ VTKM_EXEC OutType CellEdgeRatioMetric(const vtkm::IdComponent& numPts,
const Edge TetEdges[6] = { pts[1] - pts[0], pts[2] - pts[1], pts[0] - pts[2],
pts[3] - pts[0], pts[3] - pts[1], pts[3] - pts[2] };
return vtkm::exec::cellmetrics::ComputeEdgeRatio<OutType>(vtkm::make_VecC(TetEdges, numEdges));
return vtkm::worklet::cellmetrics::ComputeEdgeRatio<OutType>(vtkm::make_VecC(TetEdges, numEdges));
}
@ -229,7 +230,7 @@ VTKM_EXEC OutType CellEdgeRatioMetric(const vtkm::IdComponent& numPts,
pts[5] - pts[4], pts[6] - pts[5], pts[7] - pts[6], pts[4] - pts[7],
pts[4] - pts[0], pts[5] - pts[1], pts[6] - pts[2], pts[7] - pts[3] };
return vtkm::exec::cellmetrics::ComputeEdgeRatio<OutType>(vtkm::make_VecC(HexEdges, numEdges));
return vtkm::worklet::cellmetrics::ComputeEdgeRatio<OutType>(vtkm::make_VecC(HexEdges, numEdges));
}
@ -257,7 +258,8 @@ VTKM_EXEC OutType CellEdgeRatioMetric(const vtkm::IdComponent& numPts,
pts[4] - pts[3], pts[5] - pts[4], pts[3] - pts[5],
pts[3] - pts[0], pts[4] - pts[1], pts[5] - pts[2] };
return vtkm::exec::cellmetrics::ComputeEdgeRatio<OutType>(vtkm::make_VecC(WedgeEdges, numEdges));
return vtkm::worklet::cellmetrics::ComputeEdgeRatio<OutType>(
vtkm::make_VecC(WedgeEdges, numEdges));
}
// Compute the edge ratio of a pyramid.
@ -286,14 +288,14 @@ VTKM_EXEC OutType CellEdgeRatioMetric(const vtkm::IdComponent& numPts,
pts[4] - pts[0], pts[4] - pts[1], pts[4] - pts[2], pts[4] - pts[3]
};
return vtkm::exec::cellmetrics::ComputeEdgeRatio<OutType>(
return vtkm::worklet::cellmetrics::ComputeEdgeRatio<OutType>(
vtkm::make_VecC(PyramidEdges, numEdges));
}
} // namespace cellmetrics
} // namespace exec
} // namespace worklet
} // namespace vtkm
#endif // vtk_m_exec_cellmetrics_CellEdgeRatioMetric_h
#endif // vtk_m_worklet_cellmetrics_CellEdgeRatioMetric_h

10
vtkm/worklet/cellmetrics/CellJacobianMetric.h
View File

@ -18,8 +18,8 @@
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#ifndef vtk_m_exec_cellmetrics_Jacobian_h
#define vtk_m_exec_cellmetrics_Jacobian_h
#ifndef vtk_m_worklet_cellmetrics_Jacobian_h
#define vtk_m_worklet_cellmetrics_Jacobian_h
/*
* Mesh quality metric functions that computes the Jacobian of mesh cells.
@ -46,7 +46,7 @@
namespace vtkm
{
namespace exec
namespace worklet
{
namespace cellmetrics
{
@ -181,7 +181,7 @@ VTKM_EXEC OutType CellJacobianMetric(const vtkm::IdComponent& numPts,
return q;
}
} // namespace cellmetrics
} // namespace exec
} // namespace worklet
} // namespace vtkm
#endif // vtk_m_exec_cellmetrics_CellJacobianMetric_h
#endif // vtk_m_worklet_cellmetrics_CellJacobianMetric_h

185
vtkm/worklet/cellmetrics/CellMaxAngleMetric.h Normal file
View File

@ -0,0 +1,185 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//
// Copyright 2018 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
// Copyright 2018 UT-Battelle, LLC.
// Copyright 2018 Los Alamos National Security.
//
// Under the terms of Contract DE-NA0003525 with NTESS,
// the U.S. Government retains certain rights in this software.
//
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#ifndef vtk_m_worklet_cellmetrics_CellMaxAngleMetric_h
#define vtk_m_worklet_cellmetrics_CellMaxAngleMetric_h
/*
* Mesh quality metric functions that compute the maximum angle of cell in a mesh.
** These metric computations are adapted from the VTK implementation of the Verdict library,
* which provides a set of mesh/cell metrics for evaluating the geometric qualities of regions
* of mesh spaces.
** See: The Verdict Library Reference Manual (for per-cell-type metric formulae)
* See: vtk/ThirdParty/verdict/vtkverdict (for VTK code implementation of this metric)
*/
#include "TypeOfCellHexahedral.h"
#include "TypeOfCellQuadrilateral.h"
#include "TypeOfCellTetrahedral.h"
#include "TypeOfCellTriangle.h"
#include "vtkm/CellShape.h"
#include "vtkm/CellTraits.h"
#include "vtkm/VecTraits.h"
#include "vtkm/VectorAnalysis.h"
#include "vtkm/exec/FunctorBase.h"
#define UNUSED(expr) (void)(expr);
namespace vtkm
{
namespace worklet
{
namespace cellmetrics
{
// ========================= Unsupported cells ==================================
// By default, cells have zero shape unless the shape type template is specialized below.
template <typename OutType, typename PointCoordVecType, typename CellShapeType>
VTKM_EXEC OutType CellMaxAngleMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
CellShapeType shape,
const vtkm::exec::FunctorBase&)
{
UNUSED(numPts);
UNUSED(pts);
UNUSED(shape);
return OutType(-1.0);
}
// ========================= 2D cells ==================================
// Compute the maximum angle of a triangle.
// Formula: q = max( arccos((Ln dot Ln+1)/(||Ln|| * ||Ln+1||))(180º/π) for n 0,1, and 2 )
// - L3 = L0
// - if any edge has length 0, return q = 360º
// - All angle measurements are in degrees
// q equals 60 for a unit triangle
// Acceptable range: [30º, 60º]
// Normal Range: [0º, 360º]
// Full range: [0º, 360º]
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellMaxAngleMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagTriangle,
const vtkm::exec::FunctorBase& worklet)
{
if (numPts != 3)
{
worklet.RaiseError("Minimum angle metric (triangle) requires 3 points.");
return OutType(0.0);
}
using Scalar = OutType;
using CollectionOfPoints = PointCoordVecType;
using Vector = typename PointCoordVecType::ComponentType;
const Scalar l0 = GetTriangleL0Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l1 = GetTriangleL1Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l2 = GetTriangleL2Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
if (l0 <= Scalar(0.0) || l1 <= Scalar(0.0) || l2 <= Scalar(0.0))
{
return Scalar(0.0);
}
const Scalar oneEightyOverPi = (Scalar)57.2957795131;
const Scalar two(2.0);
const Scalar q0 = vtkm::ACos(((l1 * l1) + (l2 * l2) - (l0 * l0)) / (two * l1 * l2));
const Scalar q1 = vtkm::ACos(((l2 * l2) + (l0 * l0) - (l1 * l1)) / (two * l2 * l0));
const Scalar q2 = vtkm::ACos(((l0 * l0) + (l1 * l1) - (l2 * l2)) / (two * l0 * l1));
const Scalar q = vtkm::Max(q0, vtkm::Max(q1, q2));
const Scalar qInDegrees = q * oneEightyOverPi;
return qInDegrees;
}
// Compute the max angle of a quadrilateral.
// Formula: q = max( Ai for i 0,1,2, and 3 )
// - L4 = L0
// - Ai = -1^Si arccos(-1(Li dot Li+1)/(||Li||||Li+1||) )(180/π) + 360º*Si
// - if ||Li|| <= FLOAT_MIN or ||Li+1|| <= FLOAT_MIN, return q = 360º
// q = 90º for a unit square
// Acceptable range: [45º, 90º]
// Normal Range: [0º, 90º]
// Full range: [0º, 360º]
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellMaxAngleMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagQuad,
const vtkm::exec::FunctorBase& worklet)
{
if (numPts != 4)
{
worklet.RaiseError("Minimum angle metric(quad) requires 4 points.");
return OutType(0.0);
}
using Scalar = OutType;
using CollectionOfPoints = PointCoordVecType;
using Vector = typename PointCoordVecType::ComponentType;
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);
if (l0 <= Scalar(0.0) || l1 <= Scalar(0.0) || l2 <= Scalar(0.0) || l3 <= Scalar(0.0))
{
return Scalar(0.0);
}
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 s0 = alpha0 < Scalar(0.0) ? Scalar(1.0) : Scalar(0.0);
const Scalar s1 = alpha1 < Scalar(0.0) ? Scalar(1.0) : Scalar(0.0);
const Scalar s2 = alpha2 < Scalar(0.0) ? Scalar(1.0) : Scalar(0.0);
const Scalar s3 = alpha3 < Scalar(0.0) ? Scalar(1.0) : Scalar(0.0);
const Vector L0 = GetQuadL0<Scalar, Vector, CollectionOfPoints>(pts);
const Vector L1 = GetQuadL1<Scalar, Vector, CollectionOfPoints>(pts);
const Vector L2 = GetQuadL2<Scalar, Vector, CollectionOfPoints>(pts);
const Vector L3 = GetQuadL3<Scalar, Vector, CollectionOfPoints>(pts);
// This angle is given in degrees, not radians. The verdict definition (1) converts to degrees, (2) gives co(terminal) angles, and (3) takes the min/max.
// Further, it combines steps 1 & 2 into a single expression using clever (-1)^power flags.
const Scalar neg1(-1.0);
const Scalar oneEightyOverPi = (Scalar)57.2957795131; // ~ 180/pi
const Scalar threeSixty(360.0);
const Scalar q0 =
(vtkm::Pow(neg1, s0) * vtkm::ACos(neg1 * ((vtkm::Dot(L0, L1)) / (l0 * l1))) * oneEightyOverPi) +
(threeSixty * s0);
const Scalar q1 =
(vtkm::Pow(neg1, s1) * vtkm::ACos(neg1 * ((vtkm::Dot(L1, L2)) / (l1 * l2))) * oneEightyOverPi) +
(threeSixty * s1);
const Scalar q2 =
(vtkm::Pow(neg1, s2) * vtkm::ACos(neg1 * ((vtkm::Dot(L2, L3)) / (l2 * l3))) * oneEightyOverPi) +
(threeSixty * s2);
const Scalar q3 =
(vtkm::Pow(neg1, s3) * vtkm::ACos(neg1 * ((vtkm::Dot(L3, L0)) / (l3 * l0))) * oneEightyOverPi) +
(threeSixty * s3);
const Scalar q = vtkm::Max(q0, vtkm::Max(q1, vtkm::Max(q2, q3)));
return q;
}
} // namespace cellmetrics
} // namespace worklet
} // namespace vtkm
#endif // vtk_m_worklet_cellmetrics_CellMaxAngleMetric_h

384
vtkm/worklet/cellmetrics/CellMaxAspectFrobeniusMetric.h Normal file
View File

@ -0,0 +1,384 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//
// Copyright 2014 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
// Copyright 2014 UT-Battelle, LLC.
// Copyright 2014 Los Alamos National Security.
//
// Under the terms of Contract DE-NA0003525 with NTESS,
// the U.S. Government retains certain rights in this software.
//
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#ifndef vtk_m_worklet_cellmetrics_CellMaxAspectFrobeniusMetric_h
#define vtk_m_worklet_cellmetrics_CellMaxAspectFrobeniusMetric_h
/*
* Mesh quality metric functions that compute the maximum aspect frobenius of certain mesh cells,
* each of which are composed of two or more triangles or tetrahedrons. The output aspect metric
* value is the maximum among all triangles or tetrahedrons.
* The aspect frobenius metric generally measures the degree of regularity of a cell, with
* a value of 1 representing a regular cell.
** These metric computations are adapted from the VTK implementation of the Verdict library,
* which provides a set of mesh/cell metrics for evaluating the geometric qualities of regions
* of mesh spaces.
** The maximum aspect frobenious computations for pyramid and wedge cell types are not defined in the
* VTK implementation, but are provided here.
** See: The Verdict Library Reference Manual (for per-cell-type metric formulae)
* See: vtk/ThirdParty/verdict/vtkverdict (for VTK code implementation of this metric)
*/
#include "vtkm/CellShape.h"
#include "vtkm/CellTraits.h"
#include "vtkm/VecTraits.h"
#include "vtkm/VectorAnalysis.h"
#include "vtkm/exec/FunctorBase.h"
#include "vtkm/worklet/cellmetrics/CellAspectFrobeniusMetric.h"
#define UNUSED(expr) (void)(expr);
namespace vtkm
{
namespace worklet
{
namespace cellmetrics
{
using FloatType = vtkm::FloatDefault;
//This approximates the aspect frobenius of a tetrahedron, except for slight
//mathematical differences. In the standard aspect frobenius metric, a tetrahedron
//is compared to a reference right equilateral tetrahedron. However, in the max
//aspect frobenius metric of hexahedrons, the component tetrahedrons are compared
//to reference right isoceles tetrahedrons. Thus, some of the calculations differ
//to account for the change in reference tetrahedron. This condition computation
//is not to be confused with the separate CellConditionMetric metric, but is similar
//in computation.
template <typename OutType, typename VecType>
VTKM_EXEC OutType ComputeTetCondition(const VecType edges[])
{
//Compute the determinant/volume of the reference tet.
//(right isosceles tet for max aspect frobenius of hexs, pyramids, and wedges)
OutType det = (OutType)vtkm::Dot(edges[0], vtkm::Cross(edges[1], edges[2]));
if (det <= vtkm::NegativeInfinity<OutType>())
return vtkm::Infinity<OutType>();
OutType term1 =
vtkm::Dot(edges[0], edges[0]) + vtkm::Dot(edges[1], edges[1]) + vtkm::Dot(edges[2], edges[2]);
VecType crosses[3] = { vtkm::Cross(edges[0], edges[1]),
vtkm::Cross(edges[1], edges[2]),
vtkm::Cross(edges[2], edges[0]) };
OutType term2 = vtkm::Dot(crosses[0], crosses[0]) + vtkm::Dot(crosses[1], crosses[1]) +
vtkm::Dot(crosses[2], crosses[2]);
return vtkm::Sqrt(term1 * term2) / det;
}
// ========================= Unsupported cells ==================================
// By default, cells have undefined aspect frobenius unless the shape type template is specialized below.
template <typename OutType, typename PointCoordVecType, typename CellShapeType>
VTKM_EXEC OutType CellMaxAspectFrobeniusMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
CellShapeType shape,
const vtkm::exec::FunctorBase& worklet)
{
UNUSED(numPts);
UNUSED(pts);
UNUSED(shape);
worklet.RaiseError(
"Shape type template must be specialized to compute the max aspect frobenius metric.");
return OutType(0.0);
}
//If the polygon has 3 vertices or 4 vertices, then just call
//the functions for Triangle and Quad cell types. Otherwise,
//this metric is not supported for (n>4)-vertex polygons, such
//as pentagons or hexagons, or (n<3)-vertex polygons, such as lines or points.
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellMaxAspectFrobeniusMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagPolygon,
const vtkm::exec::FunctorBase& worklet)
{
if (numPts == 3)
return vtkm::worklet::cellmetrics::CellAspectFrobeniusMetric<OutType>(
numPts, pts, vtkm::CellShapeTagTriangle(), worklet);
if (numPts == 4)
return CellMaxAspectFrobeniusMetric<OutType>(numPts, pts, vtkm::CellShapeTagQuad(), worklet);
else
{
worklet.RaiseError(
"Max aspect frobenius metric is not supported for (n<3)- or (n>4)-vertex polygons.");
return OutType(0.0);
}
}
//The max aspect frobenius metric is not supported for lines/edges.
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellMaxAspectFrobeniusMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagLine,
const vtkm::exec::FunctorBase& worklet)
{
UNUSED(numPts);
UNUSED(pts);
worklet.RaiseError("Max aspect frobenius metric is not supported for lines/edges.");
return OutType(0.0);
}
//The max aspect frobenius metric is not uniquely defined for triangles,
//since the standard aspect frobenius metric is used for triangles.
//Thus, this implementation simply calls the aspect frobenius metric.
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellMaxAspectFrobeniusMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagTriangle,
const vtkm::exec::FunctorBase& worklet)
{
if (numPts != 3)
{
worklet.RaiseError("Max aspect frobenius metric(triangle) requires 3 points.");
return OutType(0.0);
}
return vtkm::worklet::cellmetrics::CellAspectFrobeniusMetric<OutType>(
numPts, pts, vtkm::CellShapeTagTriangle(), worklet);
}
//The max aspect frobenius metric is not supported for pyramids.
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellMaxAspectFrobeniusMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagPyramid,
const vtkm::exec::FunctorBase& worklet)
{
UNUSED(numPts);
UNUSED(pts);
worklet.RaiseError("Max aspect frobenius metric is not supported for pyramids.");
return OutType(0.0);
}
// ========================= 2D cells ==================================
// Computes the max aspect frobenius of a quad.
// Formula: The maximum aspect frobenius metric among the four triangles formed
// at the four corner points of the quad. Given a corner point, two other points are
// chosen in a uniform, counter-clockwise manner to form a triangle. The aspect frobenius metric
// is computed on this triangle. The maximum among this four computed triangle metrics
// is returned as output.
// Equals 1 for a unit square.
// Acceptable range: [1,1.3]
// Full range: [1,FLOAT_MAX]
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellMaxAspectFrobeniusMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagQuad,
const vtkm::exec::FunctorBase& worklet)
{
if (numPts != 4)
{
worklet.RaiseError("Max aspect frobenius metric(quad) requires 4 points.");
return OutType(0.0);
}
//The 4 edges of a quad
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] };
FloatType a2 = vtkm::MagnitudeSquared(QuadEdges[0]);
FloatType b2 = vtkm::MagnitudeSquared(QuadEdges[1]);
FloatType c2 = vtkm::MagnitudeSquared(QuadEdges[2]);
FloatType d2 = vtkm::MagnitudeSquared(QuadEdges[3]);
//Compute the length of the cross product for each of the 4 reference triangles.
//The result is twice the area of the triangle.
FloatType ab = vtkm::Magnitude(vtkm::Cross(QuadEdges[0], QuadEdges[1]));
FloatType bc = vtkm::Magnitude(vtkm::Cross(QuadEdges[1], QuadEdges[2]));
FloatType cd = vtkm::Magnitude(vtkm::Cross(QuadEdges[2], QuadEdges[3]));
FloatType da = vtkm::Magnitude(vtkm::Cross(QuadEdges[3], QuadEdges[0]));
if (ab < vtkm::NegativeInfinity<OutType>() || bc < vtkm::NegativeInfinity<OutType>() ||
cd < vtkm::NegativeInfinity<OutType>() || da < vtkm::NegativeInfinity<OutType>())
return vtkm::Infinity<OutType>();
FloatType qmax = (a2 + b2) / ab; //Initial max aspect frobenius (triangle 0)
FloatType qcur = (b2 + c2) / bc; //Keep checking/updating max (triangles 1 - 3)
qmax = qmax > qcur ? qmax : qcur;
qcur = (c2 + d2) / cd;
qmax = qmax > qcur ? qmax : qcur;
qcur = (d2 + a2) / da;
qmax = qmax > qcur ? qmax : qcur;
OutType max_aspect_frobenius = 0.5f * (OutType)qmax;
if (max_aspect_frobenius > 0)
vtkm::Min(max_aspect_frobenius, vtkm::Infinity<OutType>());
return vtkm::Max(max_aspect_frobenius, vtkm::NegativeInfinity<OutType>());
}
// ============================= 3D Volume cells ==================================i
// Computes the aspect frobenius of a tetrahedron.
// Formula: Sum of lengths of 3 edges, divided by a multiple of the triangle area.
// Equals 1 for a right regular tetrahedron (4 equilateral triangles).
// Acceptable range: [1,1.3]
// Full range: [1,FLOAT_MAX]
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellMaxAspectFrobeniusMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagTetra,
const vtkm::exec::FunctorBase& worklet)
{
if (numPts != 4)
{
worklet.RaiseError("Max aspect frobenius metric(tetrahedron) requires 4 points.");
return OutType(0.0);
}
return vtkm::worklet::cellmetrics::CellAspectFrobeniusMetric<OutType, PointCoordVecType>(
numPts, pts, vtkm::CellShapeTagTetra(), worklet);
}
// Computes the maximum aspect frobenius of a hexahedron.
// Formula: The maximum aspect frobenius metric among the eight tetrahedrons formed
// at the eight corner points of the hex. Given a corner point, three other points are
// chosen in a uniform, counter-clockwise manner to form a tetrahedron. The aspect frobenius metric
// is computed on this tet, with respect to a reference right isosceles tet. The maximum among
// these eight computed tet metrics is returned as output.
// Equals 1 for a unit cube (right isosceles tet formed at all 8 corner points).
// Acceptable range: [1,1.3]
// Full range: [1,FLOAT_MAX]
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellMaxAspectFrobeniusMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagHexahedron,
const vtkm::exec::FunctorBase& worklet)
{
if (numPts != 8)
{
worklet.RaiseError("Max aspect frobenius metric(hexahedron) requires 8 points.");
return OutType(0.0);
}
using Edge = typename PointCoordVecType::ComponentType;
//8 tets: one constructed at each different corner of the hex
//For each tet: Two base edges and one vertical edge, used to compute the tet volume
const Edge TetEdges[8][3] = {
{
pts[1] - pts[0], //Base edge 1
pts[3] - pts[0], //Base edge 2
pts[4] - pts[0] //Vertical edge 3
}, //tet 0
{ pts[2] - pts[1], pts[0] - pts[1], pts[5] - pts[1] }, //tet 1
{ pts[3] - pts[2], pts[1] - pts[2], pts[6] - pts[2] }, //tet 2
{ pts[0] - pts[3], pts[2] - pts[3], pts[7] - pts[3] }, //tet 3
{ pts[7] - pts[4], pts[5] - pts[4], pts[0] - pts[4] }, //tet 4
{ pts[4] - pts[5], pts[6] - pts[5], pts[1] - pts[5] }, //tet 5
{ pts[5] - pts[6], pts[7] - pts[6], pts[2] - pts[6] }, //tet 6
{ pts[6] - pts[7], pts[4] - pts[7], pts[3] - pts[7] } //tet 7
};
//For each tet, compute the condition metric, which approximates the deviation of the
//tet's volume to that of a right isoceles tetrahedron. Return the maximum metric value
//among all 8 tets as the maximum aspect frobenius.
OutType max_aspect_frobenius = ComputeTetCondition<OutType, Edge>(TetEdges[0]);
OutType curr;
for (vtkm::Id i = 1; i < 8; i++)
{
curr = ComputeTetCondition<OutType, Edge>(TetEdges[i]);
if (curr <= 0.0)
return vtkm::Infinity<OutType>();
if (curr > max_aspect_frobenius)
max_aspect_frobenius = curr;
}
max_aspect_frobenius /= 3.0;
if (max_aspect_frobenius > 0)
return vtkm::Min(max_aspect_frobenius, vtkm::Infinity<OutType>());
return vtkm::Max(max_aspect_frobenius, OutType(0.0));
}
// Computes the maximum aspect frobenius of a wedge.
// Formula: The maximum aspect frobenius metric among the six tetrahedrons formed
// from the six corner points of the two triangular faces. Given a corner point, three other points are
// chosen in a uniform, counter-clockwise manner to form a tetrahedron. The aspect frobenius metric
// is computed on this tet, with respect to an equilateral tet. The maximum among
// these six computed tet metrics is returned as output.
// Equals 1 for a unit wedge (two equilateral triangles of unit edge length and 3 unit squares).
// Acceptable range: [1,1.3]
// Full range: [1,FLOAT_MAX]
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellMaxAspectFrobeniusMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagWedge,
const vtkm::exec::FunctorBase& worklet)
{
if (numPts != 6)
{
worklet.RaiseError("Max aspect frobenius metric(wedge) requires 6 points.");
return OutType(0.0);
}
using Vector = typename PointCoordVecType::ComponentType;
//Four positively-oriented points of each tet
const vtkm::Vec<Vector, 4> tetras[6] = {
vtkm::Vec<Vector, 4>(pts[0], pts[1], pts[2], pts[3]), //tet 0
vtkm::Vec<Vector, 4>(pts[1], pts[2], pts[0], pts[4]), //tet 1
vtkm::Vec<Vector, 4>(pts[2], pts[0], pts[1], pts[5]), //tet 2
vtkm::Vec<Vector, 4>(pts[3], pts[5], pts[4], pts[0]), //tet 3
vtkm::Vec<Vector, 4>(pts[4], pts[3], pts[5], pts[1]), //tet 4
vtkm::Vec<Vector, 4>(pts[5], pts[4], pts[3], pts[2]) //tet 5
};
//For each tet, call the aspect frobenius metric.
//Return the maximum metric value among all 6 tets as the maximum aspect frobenius.
const vtkm::IdComponent tetPts = 4;
OutType max_aspect_frobenius = vtkm::worklet::cellmetrics::CellAspectFrobeniusMetric<OutType>(
tetPts, tetras[0], vtkm::CellShapeTagTetra(), worklet);
OutType curr;
for (vtkm::Id i = 1; i < 6; i++)
{
curr = vtkm::worklet::cellmetrics::CellAspectFrobeniusMetric<OutType>(
tetPts, tetras[i], vtkm::CellShapeTagTetra(), worklet);
if (curr > max_aspect_frobenius)
max_aspect_frobenius = curr;
}
//Divide by metric value of unit wedge (normalization)
max_aspect_frobenius /= 1.16477;
if (max_aspect_frobenius > 0)
return vtkm::Min(max_aspect_frobenius, vtkm::Infinity<OutType>());
return vtkm::Max(max_aspect_frobenius, vtkm::NegativeInfinity<OutType>());
return OutType(0.0);
}
} // namespace cellmetrics
} // namespace worklet
} // namespace vtkm
#endif // vtk_m_worklet_cellmetrics_CellMaxAspectFrobeniusMetric_h

122
vtkm/worklet/cellmetrics/CellMaxDiagonalMetric.h Normal file
View File

@ -0,0 +1,122 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//
// Copyright 2018 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
// Copyright 2018 UT-Battelle, LLC.
// Copyright 2018 Los Alamos National Security.
//
// Under the terms of Contract DE-NA0003525 with NTESS,
// the U.S. Government retains certain rights in this software.
//
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#ifndef vtk_m_worklet_cellmetrics_Max_Diagonal_h
#define vtk_m_worklet_cellmetrics_Max_Diagonal_h
/*
* Mesh quality metric functions that compute the Oddy of mesh cells.
** These metric computations are adapted from the VTK implementation of the Verdict library,
* which provides a set of mesh/cell metrics for evaluating the geometric qualities of regions
* of mesh spaces.
** See: The Verdict Library Reference Manual (for per-cell-type metric formulae)
* See: vtk/ThirdParty/verdict/vtkverdict (for VTK code implementation of this metric)
*/
#include "TypeOfCellHexahedral.h"
#include "TypeOfCellQuadrilateral.h"
#include "TypeOfCellTetrahedral.h"
#include "TypeOfCellTriangle.h"
#include "vtkm/CellShape.h"
#include "vtkm/CellTraits.h"
#include "vtkm/VecTraits.h"
#include "vtkm/VectorAnalysis.h"
#include "vtkm/exec/FunctorBase.h"
#define UNUSED(expr) (void)(expr);
namespace vtkm
{
namespace worklet
{
namespace cellmetrics
{
// ========================= Unsupported cells ==================================
// By default, cells have zero shape unless the shape type template is specialized below.
template <typename OutType, typename PointCoordVecType, typename CellShapeType>
VTKM_EXEC OutType CellMaxDiagonalMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
CellShapeType shape,
const vtkm::exec::FunctorBase&)
{
UNUSED(numPts);
UNUSED(pts);
UNUSED(shape);
return OutType(-1.0);
}
// ============================= 3D Volume cells ==================================
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellMaxDiagonalMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagHexahedron,
const vtkm::exec::FunctorBase& worklet)
{
if (numPts != 8)
{
worklet.RaiseError("Max diagonal metric(hexahedron) requires 8 points.");
return OutType(0.0);
}
using Scalar = OutType;
Scalar temp[3], diag[4];
vtkm::IdComponent i(0);
//lengths^2 f diag nals
for (i = 0; i < 3; i++)
{
temp[i] = pts[6][i] - pts[0][i];
temp[i] = temp[i] * temp[i];
}
diag[0] = vtkm::Sqrt(temp[0] + temp[1] + temp[2]);
for (i = 0; i < 3; i++)
{
temp[i] = pts[4][i] - pts[2][i];
temp[i] = temp[i] * temp[i];
}
diag[1] = vtkm::Sqrt(temp[0] + temp[1] + temp[2]);
for (i = 0; i < 3; i++)
{
temp[i] = pts[7][i] - pts[1][i];
temp[i] = temp[i] * temp[i];
}
diag[2] = vtkm::Sqrt(temp[0] + temp[1] + temp[2]);
for (i = 0; i < 3; i++)
{
temp[i] = pts[5][i] - pts[3][i];
temp[i] = temp[i] * temp[i];
}
diag[3] = vtkm::Sqrt(temp[0] + temp[1] + temp[2]);
Scalar diagonal = diag[0];
for (i = 1; i < 4; i++)
{
diagonal = vtkm::Max(diagonal, diag[i]);
}
return Scalar(diagonal);
}
} // namespace cellmetrics
} // namespace worklet
} // namespace vtkm
#endif

182
vtkm/worklet/cellmetrics/CellMinAngleMetric.h Normal file
View File

@ -0,0 +1,182 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//
// Copyright 2018 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
// Copyright 2018 UT-Battelle, LLC.
// Copyright 2018 Los Alamos National Security.
//
// Under the terms of Contract DE-NA0003525 with NTESS,
// the U.S. Government retains certain rights in this software.
//
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#ifndef vtk_m_worklet_cellmetrics_CellMinAngleMetric_h
#define vtk_m_worklet_cellmetrics_CellMinAngleMetric_h
/*
* Mesh quality metric functions that compute the minimum angle of cell in a mesh.
** These metric computations are adapted from the VTK implementation of the Verdict library,
* which provides a set of mesh/cell metrics for evaluating the geometric qualities of regions
* of mesh spaces.
** See: The Verdict Library Reference Manual (for per-cell-type metric formulae)
* See: vtk/ThirdParty/verdict/vtkverdict (for VTK code implementation of this metric)
*/
#include "vtkm/CellShape.h"
#include "vtkm/CellTraits.h"
#include "vtkm/VecTraits.h"
#include "vtkm/VectorAnalysis.h"
#include "vtkm/exec/FunctorBase.h"
#define UNUSED(expr) (void)(expr);
namespace vtkm
{
namespace worklet
{
namespace cellmetrics
{
// ========================= Unsupported cells ==================================
// By default, cells have zero shape unless the shape type template is specialized below.
template <typename OutType, typename PointCoordVecType, typename CellShapeType>
VTKM_EXEC OutType CellMinAngleMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
CellShapeType shape,
const vtkm::exec::FunctorBase&)
{
UNUSED(numPts);
UNUSED(pts);
UNUSED(shape);
return OutType(-1.0);
}
// ========================= 2D cells ==================================
// Compute the minimum angle of a triangle.
// Formula: q = min( arccos((Ln dot Ln+1)/(||Ln|| * ||Ln+1||))(180º/π) for n 0,1, and 2 )
// - L3 = L0
// - if any edge has length 0, return q = 360º
// - All angle measurements are in degrees
// q equals 60 for a unit triangle
// Acceptable range: [30º, 60º]
// Normal Range: [0º, 360º]
// Full range: [0º, 360º]
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellMinAngleMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagTriangle,
const vtkm::exec::FunctorBase& worklet)
{
if (numPts != 3)
{
worklet.RaiseError("Minimum angle metric(quad) requires 3 points.");
return OutType(0.0);
}
using Scalar = OutType;
using CollectionOfPoints = PointCoordVecType;
using Vector = typename PointCoordVecType::ComponentType;
const Scalar l0 = GetTriangleL0Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l1 = GetTriangleL1Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l2 = GetTriangleL2Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
if (l0 <= Scalar(0.0) || l1 <= Scalar(0.0) || l2 <= Scalar(0.0))
{
return Scalar(0.0);
}
const Scalar oneEightyOverPi = (Scalar)57.2957795131;
const Scalar two(2.0);
const Scalar q0 = vtkm::ACos(((l1 * l1) + (l2 * l2) - (l0 * l0)) / (two * l1 * l2));
const Scalar q1 = vtkm::ACos(((l2 * l2) + (l0 * l0) - (l1 * l1)) / (two * l2 * l0));
const Scalar q2 = vtkm::ACos(((l0 * l0) + (l1 * l1) - (l2 * l2)) / (two * l0 * l1));
const Scalar q = vtkm::Min(q0, vtkm::Min(q1, q2));
const Scalar qInDegrees = q * oneEightyOverPi;
return qInDegrees;
}
// Compute the min angle of a quadrilateral.
// Formula: q = min( Ai for i 0,1,2, and 3 )
// - L4 = L0
// - Ai = -1^Si arccos(-1(Li dot Li+1)/(||Li||||Li+1||) )(180/π) + 360º*Si
// - if ||Li|| <= FLOAT_MIN or ||Li+1|| <= FLOAT_MIN, return q = 360º
// q = 90º for a unit square
// Acceptable range: [45º, 90º]
// Normal Range: [0º, 90º]
// Full range: [0º, 360º]
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellMinAngleMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagQuad,
const vtkm::exec::FunctorBase& worklet)
{
if (numPts != 4)
{
worklet.RaiseError("Minimum angle metric(quad) requires 4 points.");
return OutType(0.0);
}
using Scalar = OutType;
using CollectionOfPoints = PointCoordVecType;
using Vector = typename PointCoordVecType::ComponentType;
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);
if (l0 <= Scalar(0.0) || l1 <= Scalar(0.0) || l2 <= Scalar(0.0) || l3 <= Scalar(0.0))
{
return Scalar(0.0);
}
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 s0 = alpha0 < Scalar(0.0) ? Scalar(1.0) : Scalar(0.0);
const Scalar s1 = alpha1 < Scalar(0.0) ? Scalar(1.0) : Scalar(0.0);
const Scalar s2 = alpha2 < Scalar(0.0) ? Scalar(1.0) : Scalar(0.0);
const Scalar s3 = alpha3 < Scalar(0.0) ? Scalar(1.0) : Scalar(0.0);
const Vector L0 = GetQuadL0<Scalar, Vector, CollectionOfPoints>(pts);
const Vector L1 = GetQuadL1<Scalar, Vector, CollectionOfPoints>(pts);
const Vector L2 = GetQuadL2<Scalar, Vector, CollectionOfPoints>(pts);
const Vector L3 = GetQuadL3<Scalar, Vector, CollectionOfPoints>(pts);
// This angle is given in degrees, not radians. The verdict definition (1) converts to degrees, (2) gives co(terminal) angles, and (3) takes the min/max.
// Further, it combines steps 1 & 2 into a single expression using clever (-1)^power flags.
const Scalar neg1(-1.0);
const Scalar oneEightyOverPi = (Scalar)57.2957795131; // ~ 180/pi
const Scalar threeSixty(360.0);
const Scalar q0 =
(vtkm::Pow(neg1, s0) * vtkm::ACos(neg1 * ((vtkm::Dot(L0, L1)) / (l0 * l1))) * oneEightyOverPi) +
(threeSixty * s0);
const Scalar q1 =
(vtkm::Pow(neg1, s1) * vtkm::ACos(neg1 * ((vtkm::Dot(L1, L2)) / (l1 * l2))) * oneEightyOverPi) +
(threeSixty * s1);
const Scalar q2 =
(vtkm::Pow(neg1, s2) * vtkm::ACos(neg1 * ((vtkm::Dot(L2, L3)) / (l2 * l3))) * oneEightyOverPi) +
(threeSixty * s2);
const Scalar q3 =
(vtkm::Pow(neg1, s3) * vtkm::ACos(neg1 * ((vtkm::Dot(L3, L0)) / (l3 * l0))) * oneEightyOverPi) +
(threeSixty * s3);
const Scalar q = vtkm::Min(q0, vtkm::Min(q1, vtkm::Min(q2, q3)));
return q;
}
} // namespace cellmetrics
} // namespace worklet
} // namespace vtkm
#endif // vtk_m_worklet_cellmetrics_CellEdgeRatioMetric_h

122
vtkm/worklet/cellmetrics/CellMinDiagonalMetric.h Normal file
View File

@ -0,0 +1,122 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//
// Copyright 2018 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
// Copyright 2018 UT-Battelle, LLC.
// Copyright 2018 Los Alamos National Security.
//
// Under the terms of Contract DE-NA0003525 with NTESS,
// the U.S. Government retains certain rights in this software.
//
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#ifndef vtk_m_worklet_cellmetrics_Min_Diagonal_h
#define vtk_m_worklet_cellmetrics_Min_Diagonal_h
/*
* Mesh quality metric functions that compute the Oddy of mesh cells.
** These metric computations are adapted from the VTK implementation of the Verdict library,
* which provides a set of mesh/cell metrics for evaluating the geometric qualities of regions
* of mesh spaces.
** See: The Verdict Library Reference Manual (for per-cell-type metric formulae)
* See: vtk/ThirdParty/verdict/vtkverdict (for VTK code implementation of this metric)
*/
#include "TypeOfCellHexahedral.h"
#include "TypeOfCellQuadrilateral.h"
#include "TypeOfCellTetrahedral.h"
#include "TypeOfCellTriangle.h"
#include "vtkm/CellShape.h"
#include "vtkm/CellTraits.h"
#include "vtkm/VecTraits.h"
#include "vtkm/VectorAnalysis.h"
#include "vtkm/exec/FunctorBase.h"
#define UNUSED(expr) (void)(expr);
namespace vtkm
{
namespace worklet
{
namespace cellmetrics
{
// ========================= Unsupported cells ==================================
// By default, cells have zero shape unless the shape type template is specialized below.
template <typename OutType, typename PointCoordVecType, typename CellShapeType>
VTKM_EXEC OutType CellMinDiagonalMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
CellShapeType shape,
const vtkm::exec::FunctorBase&)
{
UNUSED(numPts);
UNUSED(pts);
UNUSED(shape);
return OutType(-1.0);
}
// ============================= 3D Volume cells ==================================
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellMinDiagonalMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagHexahedron,
const vtkm::exec::FunctorBase& worklet)
{
if (numPts != 8)
{
worklet.RaiseError("Min diagonal metric(hexahedron) requires 8 points.");
return OutType(0.0);
}
using Scalar = OutType;
Scalar temp[3], diag[4];
vtkm::IdComponent i(0);
//lengths^2 f diag nals
for (i = 0; i < 3; i++)
{
temp[i] = pts[6][i] - pts[0][i];
temp[i] = temp[i] * temp[i];
}
diag[0] = vtkm::Sqrt(temp[0] + temp[1] + temp[2]);
for (i = 0; i < 3; i++)
{
temp[i] = pts[4][i] - pts[2][i];
temp[i] = temp[i] * temp[i];
}
diag[1] = vtkm::Sqrt(temp[0] + temp[1] + temp[2]);
for (i = 0; i < 3; i++)
{
temp[i] = pts[7][i] - pts[1][i];
temp[i] = temp[i] * temp[i];
}
diag[2] = vtkm::Sqrt(temp[0] + temp[1] + temp[2]);
for (i = 0; i < 3; i++)
{
temp[i] = pts[5][i] - pts[3][i];
temp[i] = temp[i] * temp[i];
}
diag[3] = vtkm::Sqrt(temp[0] + temp[1] + temp[2]);
Scalar diagonal = diag[0];
for (i = 1; i < 4; i++)
{
diagonal = vtkm::Min(diagonal, diag[i]);
}
return Scalar(diagonal);
}
} // namespace cellmetrics
} // namespace worklet
} // namespace vtkm
#endif

234
vtkm/worklet/cellmetrics/CellOddyMetric.h Normal file
View File

@ -0,0 +1,234 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//
// Copyright 2018 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
// Copyright 2018 UT-Battelle, LLC.
// Copyright 2018 Los Alamos National Security.
//
// Under the terms of Contract DE-NA0003525 with NTESS,
// the U.S. Government retains certain rights in this software.
//
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#ifndef vtk_m_worklet_cellmetrics_Oddy_h
#define vtk_m_worklet_cellmetrics_Oddy_h
/*
* Mesh quality metric functions that compute the Oddy of mesh cells.
** These metric computations are adapted from the VTK implementation of the Verdict library,
* which provides a set of mesh/cell metrics for evaluating the geometric qualities of regions
* of mesh spaces.
** See: The Verdict Library Reference Manual (for per-cell-type metric formulae)
* See: vtk/ThirdParty/verdict/vtkverdict (for VTK code implementation of this metric)
*/
#include "TypeOfCellHexahedral.h"
#include "TypeOfCellQuadrilateral.h"
#include "TypeOfCellTetrahedral.h"
#include "TypeOfCellTriangle.h"
#include "vtkm/CellShape.h"
#include "vtkm/CellTraits.h"
#include "vtkm/VecTraits.h"
#include "vtkm/VectorAnalysis.h"
#include "vtkm/exec/FunctorBase.h"
#define UNUSED(expr) (void)(expr);
namespace vtkm
{
namespace worklet
{
namespace cellmetrics
{
// ========================= Unsupported cells ==================================
// By default, cells have zero shape unless the shape type template is specialized below.
template <typename OutType, typename PointCoordVecType, typename CellShapeType>
VTKM_EXEC OutType CellOddyMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
CellShapeType shape,
const vtkm::exec::FunctorBase&)
{
UNUSED(numPts);
UNUSED(pts);
UNUSED(shape);
return OutType(-1.0);
}
// ========================= 2D cells ==================================
// Compute the Oddy of a quadrilateral.
// Formula: for i 0 to 3: max{[(||Li||^2 - ||Li+1||^2)^2 + 4((Li * Li+1)^2)] / (2||Ni+1||^2)}
// - L4 = L0
// - '*' symbolizes the dot product of two vectors
// - Ni is the normal vector associated with each point
// Equals 0 for a unit square
// Acceptable range: [0,0.5]
// Normal range: [0,FLOAT_MAX]
// Full range: [0,FLOAT_MAX]
// Note, for loop avoided because L4 = L0
template <typename Scalar, typename Vector>
VTKM_EXEC Scalar GetQuadOddyQi(const Vector& Li, const Vector& LiPlus1, const Vector& NiPlus1)
{
const Scalar two(2.0);
const Scalar four(4.0);
const Scalar liMagnitudeSquared = vtkm::MagnitudeSquared(Li);
const Scalar liPlus1MagnitudeSquared = vtkm::MagnitudeSquared(LiPlus1);
const Scalar niPlus1MagnitudeSquared = vtkm::MagnitudeSquared(NiPlus1);
const Scalar q = (((liMagnitudeSquared - liPlus1MagnitudeSquared) *
(liMagnitudeSquared - liPlus1MagnitudeSquared)) +
(four * (vtkm::Dot(Li, LiPlus1) * vtkm::Dot(Li, LiPlus1)))) /
(two * niPlus1MagnitudeSquared);
return q;
}
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellOddyMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagQuad,
const vtkm::exec::FunctorBase& worklet)
{
if (numPts != 4)
{
worklet.RaiseError("Oddy metric(quad) requires 4 points.");
return OutType(0.0);
}
using Scalar = OutType;
using CollectionOfPoints = PointCoordVecType;
using Vector = typename PointCoordVecType::ComponentType;
const Vector L0 = GetQuadL0<Scalar, Vector, CollectionOfPoints>(pts);
const Vector L1 = GetQuadL1<Scalar, Vector, CollectionOfPoints>(pts);
const Vector L2 = GetQuadL2<Scalar, Vector, CollectionOfPoints>(pts);
const Vector L3 = GetQuadL3<Scalar, Vector, CollectionOfPoints>(pts);
const Vector N0 = GetQuadN0<Scalar, Vector, CollectionOfPoints>(pts);
const Vector N1 = GetQuadN1<Scalar, Vector, CollectionOfPoints>(pts);
const Vector N2 = GetQuadN2<Scalar, Vector, CollectionOfPoints>(pts);
const Vector N3 = GetQuadN3<Scalar, Vector, CollectionOfPoints>(pts);
if (vtkm::MagnitudeSquared(N0) <= Scalar(0.0) || vtkm::MagnitudeSquared(N1) <= Scalar(0.0) ||
vtkm::MagnitudeSquared(N2) <= Scalar(0.0) || vtkm::MagnitudeSquared(N3) <= Scalar(0.0))
{
return vtkm::Infinity<Scalar>();
}
const Scalar q0 = GetQuadOddyQi<Scalar, Vector>(L0, L1, N1);
const Scalar q1 = GetQuadOddyQi<Scalar, Vector>(L1, L2, N2);
const Scalar q2 = GetQuadOddyQi<Scalar, Vector>(L2, L3, N3);
const Scalar q3 = GetQuadOddyQi<Scalar, Vector>(L3, L0, N0);
const Scalar q = vtkm::Max(q0, vtkm::Max(q1, vtkm::Max(q2, q3)));
return q;
}
// ============================= 3D Volume cells ==================================
// Compute the Oddy of a hexahedron.
// Formula:
// Equals 0 for a unit cube
// Acceptable Range: [0, 0.5]
// Normal range: [0,FLOAT_MAX]
// Full range: [0,FLOAT_MAX]
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellOddyMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagHexahedron,
const vtkm::exec::FunctorBase& worklet)
{
if (numPts != 8)
{
worklet.RaiseError("Oddy metric(hexahedron) requires 8 points.");
return OutType(0.0);
}
//The 12 points of a hexahedron
using Edge = typename PointCoordVecType::ComponentType;
const Edge HexEdges[12] = {
pts[1] - pts[0], // 0
pts[2] - pts[1], pts[3] - pts[2],
pts[3] - pts[0], // 3
pts[4] - pts[0], pts[5] - pts[1],
pts[6] - pts[2], // 6
pts[7] - pts[3], pts[5] - pts[4],
pts[6] - pts[5], // 9
pts[7] - pts[6],
pts[7] - pts[4] // 11
};
Edge principleXAxis = HexEdges[0] + (pts[2] - pts[3]) + HexEdges[8] + (pts[6] - pts[7]);
Edge principleYAxis = (pts[3] - pts[0]) + HexEdges[1] + (pts[7] - pts[4]) + HexEdges[9];
Edge principleZAxis = HexEdges[4] + HexEdges[5] + HexEdges[6] + HexEdges[7];
Edge hexJacobianMatrices[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 } };
OutType third = (OutType)(1.0 / 3.0);
OutType negativeInfinity = vtkm::NegativeInfinity<OutType>();
OutType tempMatrix1_1, tempMatrix1_2, tempMatrix1_3, tempMatrix2_2, tempMatrix2_3, tempMatrix3_3,
determinant;
OutType firstPartOfNumerator, secondPartOfNumerator, currentOddy, maxOddy = negativeInfinity;
for (vtkm::IdComponent matrixNumber = 0; matrixNumber < 9; matrixNumber++)
{
/*
// these computations equal the value at X_Y for the matrix B,
// where B = matrix A multiplied by its transpose.
// the values at X_X are also used for the matrix multiplication AA.
// Note that the values 1_2 = 2_1, 1_3 = 3_1, and 2_3 = 3_2.
// This fact is used to optimize the computation
*/
tempMatrix1_1 =
vtkm::Dot(hexJacobianMatrices[matrixNumber][0], hexJacobianMatrices[matrixNumber][0]);
tempMatrix1_2 =
vtkm::Dot(hexJacobianMatrices[matrixNumber][0], hexJacobianMatrices[matrixNumber][1]);
tempMatrix1_3 =
vtkm::Dot(hexJacobianMatrices[matrixNumber][0], hexJacobianMatrices[matrixNumber][2]);
tempMatrix2_2 =
vtkm::Dot(hexJacobianMatrices[matrixNumber][1], hexJacobianMatrices[matrixNumber][1]);
tempMatrix2_3 =
vtkm::Dot(hexJacobianMatrices[matrixNumber][1], hexJacobianMatrices[matrixNumber][2]);
tempMatrix3_3 =
vtkm::Dot(hexJacobianMatrices[matrixNumber][2], hexJacobianMatrices[matrixNumber][2]);
determinant = vtkm::Dot(
hexJacobianMatrices[matrixNumber][0],
vtkm::Cross(hexJacobianMatrices[matrixNumber][1], hexJacobianMatrices[matrixNumber][2]));
if (determinant <= OutType(0.0))
{
return vtkm::Infinity<OutType>();
}
else
{
firstPartOfNumerator = (tempMatrix1_1 * tempMatrix1_1) + 2 * (tempMatrix1_2 * tempMatrix1_2) +
2 * (tempMatrix1_3 * tempMatrix1_3) + (tempMatrix2_2 * tempMatrix2_2) +
2 * (tempMatrix2_3 * tempMatrix2_3) + (tempMatrix3_3 * tempMatrix3_3);
secondPartOfNumerator = tempMatrix1_1 + tempMatrix2_2 + tempMatrix3_3;
secondPartOfNumerator *= secondPartOfNumerator;
secondPartOfNumerator *= third;
currentOddy = (firstPartOfNumerator - secondPartOfNumerator) /
(vtkm::Pow(determinant, OutType(4.0 * third)));
maxOddy = currentOddy > maxOddy ? currentOddy : maxOddy;
}
}
if (maxOddy > 0)
{
return vtkm::Min(maxOddy, vtkm::Infinity<OutType>()); //normal case
}
return vtkm::Max(maxOddy, negativeInfinity);
}
} // namespace cellmetrics
} // namespace worklet
} // namespace vtkm
#endif // vtk_m_worklet_cellmetrics_Oddy_h

310
vtkm/worklet/cellmetrics/CellScaledJacobianMetric.h Normal file
View File

@ -0,0 +1,310 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//
// Copyright 2018 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
// Copyright 2018 UT-Battelle, LLC.
// Copyright 2018 Los Alamos National Security.
//
// Under the terms of Contract DE-NA0003525 with NTESS,
// the U.S. Government retains certain rights in this software.
//
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#ifndef vtk_m_worklet_cellmetrics_ScaledJacobian_h
#define vtk_m_worklet_cellmetrics_ScaledJacobian_h
/*
* Mesh quality metric functions that computes the scaled jacobian of mesh cells.
* The jacobian of a cell is defined as the determinant of the Jociabian matrix
** These metric computations are adapted from the VTK implementation of the Verdict library,
* which provides a set of mesh/cell metrics for evaluating the geometric qualities of regions
* of mesh spaces.
** See: The Verdict Library Reference Manual (for per-cell-type metric formulae)
* See: vtk/ThirdParty/verdict/vtkverdict (for VTK code implementation of this metric)
*/
#include "TypeOfCellHexahedral.h"
#include "TypeOfCellQuadrilateral.h"
#include "TypeOfCellTetrahedral.h"
#include "TypeOfCellTriangle.h"
#include "vtkm/CellShape.h"
#include "vtkm/CellTraits.h"
#include "vtkm/VecTraits.h"
#include "vtkm/VectorAnalysis.h"
#include "vtkm/exec/FunctorBase.h"
#define UNUSED(expr) (void)(expr);
namespace vtkm
{
namespace worklet
{
namespace cellmetrics
{
// ========================= Unsupported cells ==================================
// By default, cells have zero shape unless the shape type template is specialized below.
template <typename OutType, typename PointCoordVecType, typename CellShapeType>
VTKM_EXEC OutType CellScaledJacobianMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
CellShapeType shape,
const vtkm::exec::FunctorBase&)
{
UNUSED(numPts);
UNUSED(pts);
UNUSED(shape);
return OutType(-2.0);
}
// ========================= 2D cells ==================================
//Compute the scaled jacobian of a triangle
//Formula: q = ((2*sqrt(3))/3) * (J/Lmax)
// - J -> jacobian, if surface normal N is center of triangle and
// and N*L2*L1 < 0, then -jacobian
// - Lmax -> max{ |L0| * |L1|, |L0| * |L2|, |L1| * |L2| }
//Equals 1 for equilateral unit triangle
//Acceptable Range: [0.5, 2*sqrt(3)/3]
//Normal Range : [ -(2*sqrt(3)/3) , 2*sqrt(3)/3]
//Full Range : [-FLOAT_MAX, FLOAT_MAX]
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellScaledJacobianMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagTriangle,
const vtkm::exec::FunctorBase& worklet)
{
if (numPts != 3)
{
worklet.RaiseError("ScaledJacobian metric(tri) requires 3 points.");
return OutType(0.0);
}
using Scalar = OutType;
using CollectionOfPoints = PointCoordVecType;
using Vector = typename PointCoordVecType::ComponentType;
const Vector l0 = GetTriangleL0<Scalar, Vector, CollectionOfPoints>(pts);
const Vector l1 = GetTriangleL1<Scalar, Vector, CollectionOfPoints>(pts);
const Vector l2 = GetTriangleL2<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar modifier = (Scalar)(2 * vtkm::Sqrt(3) / 3);
const Scalar l0_magnitude = GetTriangleL0Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l1_magnitude = GetTriangleL1Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l2_magnitude = GetTriangleL2Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l0l1_product = l0_magnitude * l1_magnitude;
const Scalar l0l2_product = l0_magnitude * l2_magnitude;
const Scalar l1l2_product = l1_magnitude * l2_magnitude;
const Scalar productMax = vtkm::Max(l0l1_product, vtkm::Max(l0l2_product, l1l2_product));
if (productMax <= Scalar(0.0))
{
return Scalar(0.0);
}
// compute jacobian of triangle
Vector TriCross = vtkm::Cross(l2, l1);
Scalar scaledJacobian = vtkm::Magnitude(TriCross);
//add all pieces together
//TODO change
Vector normalVector = (Scalar(1.0) / Scalar(3.0)) * (l0 + l1 + l2);
Vector surfaceNormalAtCenter = vtkm::TriangleNormal(normalVector, normalVector, normalVector);
if (vtkm::Dot(surfaceNormalAtCenter, TriCross) < 0)
{
scaledJacobian *= -1;
}
scaledJacobian *= modifier;
const Scalar q = scaledJacobian / productMax;
return q;
}
// Compute the scaled jacobian of a quadrilateral.
// Formula: min{J0/(L0*L3), J1/(L1*L0), J2/(L2*L1), J3/(L3*L2)}
// -Ji -> Jacobian at corner i, the intersection of the edge vectors
// it is divided by
// Equals 1 for a unit square
//Acceptable Range: [0.3, 1]
//Normal Range : [-1, 1]
// Full range : [-1, 1]
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellScaledJacobianMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagQuad,
const vtkm::exec::FunctorBase& worklet)
{
if (numPts != 4)
{
worklet.RaiseError("ScaledJacobian metric(quad) requires 4 points.");
return OutType(0.0);
}
using Scalar = OutType;
using CollectionOfPoints = PointCoordVecType;
using Vector = typename PointCoordVecType::ComponentType;
//The 4 edges of a quadrilateral
const Scalar l0_magnitude = GetQuadL0Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l1_magnitude = GetQuadL1Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l2_magnitude = GetQuadL2Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar l3_magnitude = GetQuadL3Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar negativeInfinity = vtkm::NegativeInfinity<Scalar>();
if (l0_magnitude < negativeInfinity || l1_magnitude < negativeInfinity ||
l2_magnitude < negativeInfinity || l3_magnitude < negativeInfinity)
{
return Scalar(0.0);
}
const Scalar l0l3_product = l0_magnitude * l3_magnitude;
const Scalar l1l0_product = l1_magnitude * l0_magnitude;
const Scalar l2l1_product = l2_magnitude * l1_magnitude;
const Scalar l3l2_product = l3_magnitude * l2_magnitude;
/*
3 * 0
0 * 1
1 * 2
2 * 3
*/
Scalar alpha0Scaled = GetQuadAlpha0<Scalar, Vector, CollectionOfPoints>(pts);
Scalar alpha1Scaled = GetQuadAlpha1<Scalar, Vector, CollectionOfPoints>(pts);
Scalar alpha2Scaled = GetQuadAlpha2<Scalar, Vector, CollectionOfPoints>(pts);
Scalar alpha3Scaled = GetQuadAlpha3<Scalar, Vector, CollectionOfPoints>(pts);
alpha0Scaled /= l0l3_product;
alpha1Scaled /= l1l0_product;
alpha2Scaled /= l2l1_product;
alpha3Scaled /= l3l2_product;
const Scalar q =
vtkm::Min(alpha0Scaled, vtkm::Min(alpha1Scaled, vtkm::Min(alpha2Scaled, alpha3Scaled)));
return q;
}
// ============================= 3D Volume cells ==================================
// Compute the scaled jacobian of a hexahedron.
// Formula: q = min{Ai}
// Ai -> for i 1...8 (Jacobian determinant at respective corner, divided by corresponding edge lengths
// Equals 1 for a unit cube
// Acceptable Range: [0.5, 1]
// Normal Range : [-1, 1]
// Full range : [1,FLOAT_MAX]
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellScaledJacobianMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagHexahedron,
const vtkm::exec::FunctorBase& worklet)
{
if (numPts != 8)
{
worklet.RaiseError("ScaledJacobian metric(hexahedron) requires 8 points.");
return OutType(0.0);
}
//The 12 edges of a hexahedron
using Edge = typename PointCoordVecType::ComponentType;
Edge HexEdges[12] = { pts[1] - pts[0], pts[2] - pts[1], pts[3] - pts[2], pts[3] - pts[0],
pts[4] - pts[0], pts[5] - pts[1], pts[6] - pts[2], pts[7] - pts[3],
pts[5] - pts[4], pts[6] - pts[5], pts[7] - pts[6], pts[7] - pts[4] };
Edge principleXAxis = HexEdges[0] + (pts[2] - pts[3]) + HexEdges[8] + (pts[6] - pts[7]);
Edge principleYAxis = (pts[3] - pts[0]) + HexEdges[1] + (pts[7] - pts[4]) + HexEdges[9];
Edge principleZAxis = HexEdges[4] + HexEdges[5] + HexEdges[6] + HexEdges[7];
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 } };
OutType currDeterminant, minDeterminant = vtkm::Infinity<OutType>();
FloatType lenSquared1, lenSquared2, lenSquared3, minLengthSquared = vtkm::Infinity<FloatType>();
vtkm::IdComponent matrixIndex;
for (matrixIndex = 0; matrixIndex < 9; matrixIndex++)
{
lenSquared1 = (FloatType)vtkm::MagnitudeSquared(hexMatrices[matrixIndex][0]);
minLengthSquared = lenSquared1 < minLengthSquared ? lenSquared1 : minLengthSquared;
lenSquared2 = (FloatType)vtkm::MagnitudeSquared(hexMatrices[matrixIndex][1]);
minLengthSquared = lenSquared2 < minLengthSquared ? lenSquared2 : minLengthSquared;
lenSquared3 = (FloatType)vtkm::MagnitudeSquared(hexMatrices[matrixIndex][2]);
minLengthSquared = lenSquared3 < minLengthSquared ? lenSquared3 : minLengthSquared;
vtkm::Normalize(hexMatrices[matrixIndex][0]);
vtkm::Normalize(hexMatrices[matrixIndex][1]);
vtkm::Normalize(hexMatrices[matrixIndex][2]);
currDeterminant =
(OutType)vtkm::Dot(hexMatrices[matrixIndex][0],
vtkm::Cross(hexMatrices[matrixIndex][1], hexMatrices[matrixIndex][2]));
if (currDeterminant < minDeterminant)
{
minDeterminant = currDeterminant;
}
}
if (minLengthSquared < vtkm::NegativeInfinity<FloatType>())
{
return vtkm::Infinity<OutType>();
}
OutType toReturn = minDeterminant;
if (toReturn > 0)
return vtkm::Min(toReturn, vtkm::Infinity<OutType>()); //normal case
return vtkm::Max(toReturn, vtkm::NegativeInfinity<OutType>());
}
// Compute the scaled jacobian of a tetrahedron.
// Formula: q = J*sqrt(2)/Lamda_max
// J -> jacobian,(L2 * L0) * L3
// Lamda_max -> max{ L0*L2*L3, L0*L1*L4, L1*L2*L5, L3*L4*L5}
// Equals Sqrt(2) / 2 for unit equilateral tetrahedron
// Acceptable Range: [0, FLOAT_MAX]
// Normal Range: [0, FLOAT_MAX]
// Full range: [FLOAT_MIN,FLOAT_MAX]
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellScaledJacobianMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagTetra,
const vtkm::exec::FunctorBase& worklet)
{
if (numPts != 4)
{
worklet.RaiseError("ScaledJacobian metric requires 4 points");
return OutType(0.0);
}
//the edge and side sets
using Edge = typename PointCoordVecType::ComponentType;
Edge Edges[6] = { pts[1] - pts[0], pts[2] - pts[1], pts[0] - pts[2],
pts[3] - pts[0], pts[3] - pts[1], pts[3] - pts[2] };
OutType EdgesSquared[6];
OutType jacobian = vtkm::Dot(vtkm::Cross(Edges[2], Edges[0]), Edges[3]);
// compute the scaled jacobian
OutType currSide, maxSide = vtkm::NegativeInfinity<OutType>();
vtkm::IdComponent edgeIndex, sideIndex;
for (edgeIndex = 0; edgeIndex < 6; edgeIndex++)
{
EdgesSquared[edgeIndex] = vtkm::MagnitudeSquared(Edges[edgeIndex]);
}
OutType Sides[4] = { EdgesSquared[0] * EdgesSquared[2] * EdgesSquared[3],
EdgesSquared[0] * EdgesSquared[1] * EdgesSquared[4],
EdgesSquared[1] * EdgesSquared[2] * EdgesSquared[5],
EdgesSquared[3] * EdgesSquared[4] * EdgesSquared[5] };
for (sideIndex = 0; sideIndex < 4; sideIndex++)
{
currSide = Sides[sideIndex];
maxSide = currSide > maxSide ? currSide : maxSide;
}
maxSide = vtkm::Sqrt(maxSide);
OutType toUseInCalculation = jacobian > maxSide ? jacobian : maxSide;
if (toUseInCalculation < vtkm::NegativeInfinity<OutType>())
{
return vtkm::Infinity<OutType>();
}
return (vtkm::Sqrt<OutType>(2) * jacobian) / toUseInCalculation;
}
} // namespace cellmetrics
} // namespace worklet
} // namespace vtkm
#endif // vtk_m_worklet_cellmetrics_CellEdgeRatioMetric_h

267
vtkm/worklet/cellmetrics/CellShapeMetric.h Normal file
View File

@ -0,0 +1,267 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//
// Copyright 2014 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
// Copyright 2014 UT-Battelle, LLC.
// Copyright 2014 Los Alamos National Security.
//
// Under the terms of Contract DE-NA0003525 with NTESS,
// the U.S. Government retains certain rights in this software.
//
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#ifndef vtk_m_worklet_CellShapeMetric_h
#define vtk_m_worklet_CellShapeMetric_h
/*
* Mesh quality metric functions that compute the shape, or weighted Jacobian, of mesh cells.
* The Jacobian of a cell is weighted by the condition metric value of the cell.
** These metric computations are adapted from the VTK implementation of the Verdict library,
* which provides a set of cell metrics for evaluating the geometric qualities of regions of mesh spaces.
** See: The Verdict Library Reference Manual (for per-cell-type metric formulae)
* See: vtk/ThirdParty/verdict/vtkverdict (for VTK code implementation of this metric)
*/
#include "TypeOfCellHexahedral.h"
#include "TypeOfCellQuadrilateral.h"
#include "TypeOfCellTetrahedral.h"
#include "TypeOfCellTriangle.h"
#include "vtkm/CellShape.h"
#include "vtkm/CellTraits.h"
#include "vtkm/VecTraits.h"
#include "vtkm/VectorAnalysis.h"
#include "vtkm/exec/FunctorBase.h"
#include "vtkm/worklet/cellmetrics/CellConditionMetric.h"
#include "vtkm/worklet/cellmetrics/CellJacobianMetric.h"
namespace vtkm
{
namespace worklet
{
namespace cellmetrics
{
// By default, cells have no shape unless the shape type template is specialized below.
template <typename OutType, typename PointCoordVecType, typename CellShapeType>
VTKM_EXEC OutType CellShapeMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
CellShapeType shape,
const vtkm::exec::FunctorBase&)
{
UNUSED(numPts);
UNUSED(pts);
UNUSED(shape);
return OutType(-1.0);
}
// =============================== Shape metric cells ==================================
// Compute the shape quality metric of a triangular cell.
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellShapeMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagTriangle,
const vtkm::exec::FunctorBase& worklet)
{
if (numPts != 3)
{
worklet.RaiseError("Shape metric(triangle) requires 3 points.");
return OutType(0.0);
}
using Scalar = OutType;
using CollectionOfPoints = PointCoordVecType;
const Scalar condition =
vtkm::worklet::cellmetrics::CellConditionMetric<Scalar, CollectionOfPoints>(
numPts, pts, vtkm::CellShapeTagTriangle(), worklet);
const Scalar q(1 / condition);
return q;
}
/// Compute the shape of a quadrilateral.
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellShapeMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagQuad,
const vtkm::exec::FunctorBase& worklet)
{
if (numPts != 4)
{
worklet.RaiseError("Area(quad) requires 4 points.");
return OutType(0.0);
}
using Scalar = OutType;
using CollectionOfPoints = PointCoordVecType;
using Vector = typename PointCoordVecType::ComponentType;
const Scalar two(2.0);
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 l0Squared =
vtkm::Pow(GetQuadL0Magnitude<Scalar, Vector, CollectionOfPoints>(pts), 2);
const Scalar l1Squared =
vtkm::Pow(GetQuadL1Magnitude<Scalar, Vector, CollectionOfPoints>(pts), 2);
const Scalar l2Squared =
vtkm::Pow(GetQuadL2Magnitude<Scalar, Vector, CollectionOfPoints>(pts), 2);
const Scalar l3Squared =
vtkm::Pow(GetQuadL3Magnitude<Scalar, Vector, CollectionOfPoints>(pts), 2);
const Scalar min = vtkm::Min(
(alpha0 / (l0Squared + l3Squared)),
vtkm::Min((alpha1 / (l1Squared + l0Squared)),
vtkm::Min((alpha2 / (l2Squared + l1Squared)), (alpha3 / (l3Squared + l2Squared)))));
const Scalar q(two * min);
return q;
}
// ============================= 3D cells ==================================
/// Compute the shape of a tetrahedron.
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellShapeMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagTetra,
const vtkm::exec::FunctorBase& worklet)
{
if (numPts != 4)
{
worklet.RaiseError("Shape metric(tetrahedron) requires 4 points.");
return OutType(0.0);
}
using Scalar = OutType;
using CollectionOfPoints = PointCoordVecType;
using Vector = typename PointCoordVecType::ComponentType;
const Scalar zero(0.0);
const Scalar twoThirds = (Scalar)(2.0 / 3.0);
const Scalar threeHalves(1.5);
const Scalar rtTwo = (Scalar)(vtkm::Sqrt(2.0));
const Scalar three(3.0);
const Scalar jacobian =
vtkm::worklet::cellmetrics::CellJacobianMetric<Scalar, CollectionOfPoints>(
numPts, pts, vtkm::CellShapeTagTetra(), worklet);
if (jacobian <= zero)
{
return zero;
}
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 Vector negl2 = -1 * l2;
const Scalar l0l0 = vtkm::Dot(l0, l0);
const Scalar l2l2 = vtkm::Dot(l2, l2);
const Scalar l3l3 = vtkm::Dot(l3, l3);
const Scalar l0negl2 = vtkm::Dot(l0, negl2);
const Scalar l0l3 = vtkm::Dot(l0, l3);
const Scalar negl2l3 = vtkm::Dot(negl2, l3);
const Scalar numerator = three * vtkm::Pow(jacobian * rtTwo, twoThirds);
Scalar denominator = (threeHalves * (l0l0 + l2l2 + l3l3)) - (l0negl2 + l0l3 + negl2l3);
if (denominator <= zero)
{
return zero;
}
Scalar q(numerator / denominator);
return q;
}
/// Compute the shape of a hexahedral cell.
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellShapeMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagHexahedron,
const vtkm::exec::FunctorBase& worklet)
{
if (numPts != 8)
{
worklet.RaiseError("Volume(hexahedron) requires 8 points.");
return OutType(0.0);
}
using Scalar = OutType;
using CollectionOfPoints = PointCoordVecType;
using Vector = typename PointCoordVecType::ComponentType;
const Scalar zero(0.0);
const Scalar twoThirds = (Scalar)(2.0 / 3.0);
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 A0squared =
GetHexAiNormSquared<Scalar, Vector, CollectionOfPoints>(pts, vtkm::Id(0));
const Scalar A1squared =
GetHexAiNormSquared<Scalar, Vector, CollectionOfPoints>(pts, vtkm::Id(1));
const Scalar A2squared =
GetHexAiNormSquared<Scalar, Vector, CollectionOfPoints>(pts, vtkm::Id(2));
const Scalar A3squared =
GetHexAiNormSquared<Scalar, Vector, CollectionOfPoints>(pts, vtkm::Id(3));
const Scalar A4squared =
GetHexAiNormSquared<Scalar, Vector, CollectionOfPoints>(pts, vtkm::Id(4));
const Scalar A5squared =
GetHexAiNormSquared<Scalar, Vector, CollectionOfPoints>(pts, vtkm::Id(5));
const Scalar A6squared =
GetHexAiNormSquared<Scalar, Vector, CollectionOfPoints>(pts, vtkm::Id(6));
const Scalar A7squared =
GetHexAiNormSquared<Scalar, Vector, CollectionOfPoints>(pts, vtkm::Id(7));
const Scalar A8squared =
GetHexAiNormSquared<Scalar, Vector, CollectionOfPoints>(pts, vtkm::Id(8));
if (alpha0 <= zero || alpha1 <= zero || alpha2 <= zero || alpha3 <= zero || alpha4 <= zero ||
alpha5 <= zero || alpha6 <= zero || alpha7 <= zero || alpha8 <= zero || A0squared <= zero ||
A1squared <= zero || A2squared <= zero || A3squared <= zero || A4squared <= zero ||
A5squared <= zero || A6squared <= zero || A7squared <= zero || A8squared <= zero)
{
return zero;
}
// find min according to verdict manual
Scalar min = vtkm::Pow(alpha0, twoThirds) / A0squared;
Scalar temp = vtkm::Pow(alpha1, twoThirds) / A1squared;
min = temp < min ? temp : min;
temp = vtkm::Pow(alpha2, twoThirds) / A2squared;
min = temp < min ? temp : min;
temp = vtkm::Pow(alpha3, twoThirds) / A3squared;
min = temp < min ? temp : min;
temp = vtkm::Pow(alpha4, twoThirds) / A4squared;
min = temp < min ? temp : min;
temp = vtkm::Pow(alpha5, twoThirds) / A5squared;
min = temp < min ? temp : min;
temp = vtkm::Pow(alpha6, twoThirds) / A6squared;
min = temp < min ? temp : min;
temp = vtkm::Pow(alpha7, twoThirds) / A7squared;
min = temp < min ? temp : min;
temp = vtkm::Pow(alpha8, twoThirds) / A8squared;
min = temp < min ? temp : min;
// determine the shape
Scalar q(3 * min);
return q;
}
} // cellmetrics
} // worklet
} // vtkm
#endif // vtk_m_worklet_CellShapeMetric_h

145
vtkm/worklet/cellmetrics/CellShearMetric.h Normal file
View File

@ -0,0 +1,145 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//
// Copyright 2018 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
// Copyright 2018 UT-Battelle, LLC.
// Copyright 2018 Los Alamos National Security.
//
// Under the terms of Contract DE-NA0003525 with NTESS,
// the U.S. Government retains certain rights in this software.
//
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#ifndef vtk_m_worklet_cellmetrics_CellShearMetric_h
#define vtk_m_worklet_cellmetrics_CellShearMetric_h
/*
* Mesh quality metric function that computes the shear of a cell. The shear is
* defined as the minimum of areas of the faces divided by the length of the
* indexed face multiplied by the index - 1 mod num_faces.
*
* These metric computations are adapted from the VTK implementation of the
* Verdict library, which provides a set of mesh/cell metrics for evaluating the
* geometric qualities of regions of mesh spaces.
*
* See: The Verdict Library Reference Manual (for per-cell-type metric formulae)
* See: vtk/ThirdParty/verdict/vtkverdict (for VTK code implementation of this
* metric)
*/
#include "TypeOfCellHexahedral.h"
#include "TypeOfCellQuadrilateral.h"
#include "vtkm/CellShape.h"
#include "vtkm/CellTraits.h"
#include "vtkm/VecTraits.h"
#include "vtkm/VectorAnalysis.h"
#include "vtkm/exec/FunctorBase.h"
#define UNUSED(expr) (void)(expr);
namespace vtkm
{
namespace worklet
{
namespace cellmetrics
{
// ========================= Unsupported cells ==================================
// By default, cells have zero shape unless the shape type template is specialized below.
template <typename OutType, typename PointCoordVecType, typename CellShapeType>
VTKM_EXEC OutType CellShearMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
CellShapeType shape,
const vtkm::exec::FunctorBase&)
{
UNUSED(numPts);
UNUSED(pts);
UNUSED(shape);
return OutType(-1.0);
}
// ========================= 2D cells ==================================
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellShearMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagQuad,
const vtkm::exec::FunctorBase& worklet)
{
if (numPts != 4)
{
worklet.RaiseError("Diagonal ratio metric(quad) requires 4 points.");
return OutType(0.0);
}
using Scalar = OutType;
using CollectionOfPoints = PointCoordVecType;
using Vector = typename PointCoordVecType::ComponentType;
const Scalar alpha0 = GetQuadAlpha0<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar alpha1 = GetQuadAlpha1<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar alpha2 = GetQuadAlpha2<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar alpha3 = GetQuadAlpha3<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar 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 x0 = alpha0 / (L0 * L3);
const Scalar x1 = alpha1 / (L1 * L0);
const Scalar x2 = alpha2 / (L2 * L1);
const Scalar x3 = alpha3 / (L3 * L2);
const Scalar q = vtkm::Min(x0, vtkm::Min(x1, vtkm::Min(x2, x3)));
return q;
}
// ========================= 3D cells ==================================
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellShearMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagHexahedron,
const vtkm::exec::FunctorBase& worklet)
{
if (numPts != 8)
{
worklet.RaiseError("Diagonal ratio metric(hex) requires 8 points.");
return OutType(-1.0);
}
using Scalar = OutType;
using CollectionOfPoints = PointCoordVecType;
using Vector = typename PointCoordVecType::ComponentType;
const Scalar a0 = GetHexAlphaiHat<Scalar, Vector, CollectionOfPoints>(pts, vtkm::Id(0));
const Scalar a1 = GetHexAlphaiHat<Scalar, Vector, CollectionOfPoints>(pts, vtkm::Id(1));
const Scalar a2 = GetHexAlphaiHat<Scalar, Vector, CollectionOfPoints>(pts, vtkm::Id(2));
const Scalar a3 = GetHexAlphaiHat<Scalar, Vector, CollectionOfPoints>(pts, vtkm::Id(3));
const Scalar a4 = GetHexAlphaiHat<Scalar, Vector, CollectionOfPoints>(pts, vtkm::Id(4));
const Scalar a5 = GetHexAlphaiHat<Scalar, Vector, CollectionOfPoints>(pts, vtkm::Id(5));
const Scalar a6 = GetHexAlphaiHat<Scalar, Vector, CollectionOfPoints>(pts, vtkm::Id(6));
const Scalar a7 = GetHexAlphaiHat<Scalar, Vector, CollectionOfPoints>(pts, vtkm::Id(7));
const Scalar q = vtkm::Min(
a0,
vtkm::Min(a1, vtkm::Min(a2, vtkm::Min(a3, vtkm::Min(a4, vtkm::Min(a5, vtkm::Min(a6, a7)))))));
return q;
}
} // namespace cellmetrics
} // namespace worklet
} // namespace vtkm
#endif // vtk_m_worklet_cellmetrics_CellEdgeRatioMetric_h

109
vtkm/worklet/cellmetrics/CellSkewMetric.h Normal file
View File

@ -0,0 +1,109 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//
// Copyright 2014 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
// Copyright 2014 UT-Battelle, LLC.
// Copyright 2014 Los Alamos National Security.
//
// Under the terms of Contract DE-NA0003525 with NTESS,
// the U.S. Government retains certain rights in this software.
//
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#ifndef vtk_m_worklet_CellSkewMetric_h
#define vtk_m_worklet_CellSkewMetric_h
/*
*/
#include "TypeOfCellHexahedral.h"
#include "TypeOfCellQuadrilateral.h"
#include "TypeOfCellTetrahedral.h"
#include "TypeOfCellTriangle.h"
#include "vtkm/CellShape.h"
#include "vtkm/CellTraits.h"
#include "vtkm/VecTraits.h"
#include "vtkm/VectorAnalysis.h"
#include "vtkm/exec/FunctorBase.h"
#include "vtkm/worklet/cellmetrics/CellConditionMetric.h"
namespace vtkm
{
namespace worklet
{
namespace cellmetrics
{
template <typename OutType, typename PointCoordVecType, typename CellShapeType>
VTKM_EXEC OutType CellSkewMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
CellShapeType shape,
const vtkm::exec::FunctorBase& worklet)
{
UNUSED(numPts);
UNUSED(pts);
UNUSED(shape);
UNUSED(worklet);
return OutType(-1.0);
}
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellSkewMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagHexahedron,
const vtkm::exec::FunctorBase& worklet)
{
UNUSED(numPts);
UNUSED(worklet);
using Scalar = OutType;
using Vector = typename PointCoordVecType::ComponentType;
Vector X1 = (pts[1] - pts[0]) + (pts[2] - pts[3]) + (pts[5] - pts[4]) + (pts[6] - pts[7]);
Scalar X1_mag = vtkm::Magnitude(X1);
if (X1_mag <= Scalar(0.0))
return vtkm::Infinity<Scalar>();
Vector x1 = X1 / X1_mag;
Vector X2 = (pts[3] - pts[0]) + (pts[2] - pts[1]) + (pts[7] - pts[4]) + (pts[6] - pts[5]);
Scalar X2_mag = vtkm::Magnitude(X2);
if (X2_mag <= Scalar(0.0))
return vtkm::Infinity<Scalar>();
Vector x2 = X2 / X2_mag;
Vector X3 = (pts[4] - pts[0]) + (pts[5] - pts[1]) + (pts[6] - pts[2]) + (pts[7] - pts[3]);
Scalar X3_mag = vtkm::Magnitude(X3);
if (X3_mag <= Scalar(0.0))
return vtkm::Infinity<Scalar>();
Vector x3 = X3 / X3_mag;
return vtkm::Max(vtkm::Dot(x1, x2), vtkm::Max(vtkm::Dot(x1, x3), vtkm::Dot(x2, x3)));
}
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellSkewMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagQuad,
const vtkm::exec::FunctorBase& worklet)
{
UNUSED(numPts);
UNUSED(worklet);
using Scalar = OutType;
using CollectionOfPoints = PointCoordVecType;
using Vector = typename PointCoordVecType::ComponentType;
const Vector X0 = GetQuadX0<Scalar, Vector, CollectionOfPoints>(pts);
const Vector X1 = GetQuadX1<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar X0Mag = vtkm::Magnitude(X0);
const Scalar X1Mag = vtkm::Magnitude(X1);
if (X0Mag < Scalar(0.0) || X1Mag < Scalar(0.0))
return Scalar(0.0);
const Vector x0Normalized = X0 / X0Mag;
const Vector x1Normalized = X1 / X1Mag;
const Scalar dot = vtkm::Dot(x0Normalized, x1Normalized);
return vtkm::Abs(dot);
}
}
} // worklet
} // vtkm
#endif

118
vtkm/worklet/cellmetrics/CellStretchMetric.h Normal file
View File

@ -0,0 +1,118 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//
// Copyright 2014 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
// Copyright 2014 UT-Battelle, LLC.
// Copyright 2014 Los Alamos National Security.
//
// Under the terms of Contract DE-NA0003525 with NTESS,
// the U.S. Government retains certain rights in this software.
//
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#ifndef vtk_m_worklet_cellmetrics_CellStretchMetric_h
#define vtk_m_worklet_cellmetrics_CellStretchMetric_h
/*
* Mesh quality metric functions that compute the aspect frobenius of certain mesh cells.
* The aspect frobenius metric generally measures the degree of regularity of a cell, with
* a value of 1 representing a regular cell..
** These metric computations are adapted from the VTK implementation of the Verdict library,
* which provides a set of mesh/cell metrics for evaluating the geometric qualities of regions
* of mesh spaces.
** See: The Verdict Library Reference Manual (for per-cell-type metric formulae)
* See: vtk/ThirdParty/verdict/vtkverdict (for VTK code implementation of this metric)
*/
#include "TypeOfCellHexahedral.h"
#include "TypeOfCellQuadrilateral.h"
#include "TypeOfCellTetrahedral.h"
#include "TypeOfCellTriangle.h"
#include "vtkm/CellShape.h"
#include "vtkm/CellTraits.h"
#include "vtkm/VecTraits.h"
#include "vtkm/VectorAnalysis.h"
#include "vtkm/exec/FunctorBase.h"
namespace vtkm
{
namespace worklet
{
namespace cellmetrics
{
template <typename OutType, typename PointCoordVecType, typename CellShapeType>
VTKM_EXEC OutType CellStretchMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
CellShapeType shape,
const vtkm::exec::FunctorBase& worklet)
{
UNUSED(numPts);
UNUSED(pts);
UNUSED(shape);
UNUSED(worklet);
return OutType(-1.0);
}
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellStretchMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagQuad,
const vtkm::exec::FunctorBase& worklet)
{
UNUSED(numPts);
UNUSED(worklet);
using Scalar = OutType;
using CollectionOfPoints = PointCoordVecType;
using Vector = typename PointCoordVecType::ComponentType;
const Scalar root2 = vtkm::Sqrt(Scalar(2.0));
const Scalar LMin = GetQuadLMin<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar DMax = GetQuadDMax<Scalar, Vector, CollectionOfPoints>(pts);
if (DMax <= Scalar(0.0))
{
return vtkm::Infinity<Scalar>();
}
const Scalar q = root2 * (LMin / DMax);
return q;
}
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellStretchMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagHexahedron,
const vtkm::exec::FunctorBase& worklet)
{
UNUSED(numPts);
UNUSED(worklet);
using Scalar = OutType;
using CollectionOfPoints = PointCoordVecType;
using Vector = typename PointCoordVecType::ComponentType;
const Scalar root3 = vtkm::Sqrt(Scalar(3.0));
const Scalar LMin = GetHexLMin<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar DMax = GetHexDMax<Scalar, Vector, CollectionOfPoints>(pts);
if (DMax <= Scalar(0.0))
{
return vtkm::Infinity<Scalar>();
}
const Scalar q = root3 * (LMin / DMax);
return q;
}
}
}
}
#endif // vtk_m_worklet_cellmetrics_CellStretchMetric_h

129
vtkm/worklet/cellmetrics/CellTaperMetric.h Normal file
View File

@ -0,0 +1,129 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//
// Copyright 2014 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
// Copyright 2014 UT-Battelle, LLC.
// Copyright 2014 Los Alamos National Security.
//
// Under the terms of Contract DE-NA0003525 with NTESS,
// the U.S. Government retains certain rights in this software.
//
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#ifndef vtk_m_worklet_CellTaperMetric_h
#define vtk_m_worklet_CellTaperMetric_h
/*
* Mesh quality metric functions that compute the shape, or weighted Jacobian, of mesh cells.
* The Jacobian of a cell is weighted by the condition metric value of the cell.
** These metric computations are adapted from the VTK implementation of the Verdict library,
* which provides a set of cell metrics for evaluating the geometric qualities of regions of mesh spaces.
** See: The Verdict Library Reference Manual (for per-cell-type metric formulae)
* See: vtk/ThirdParty/verdict/vtkverdict (for VTK code implementation of this metric)
*/
#include "TypeOfCellHexahedral.h"
#include "TypeOfCellQuadrilateral.h"
#include "TypeOfCellTetrahedral.h"
#include "TypeOfCellTriangle.h"
#include "vtkm/CellShape.h"
#include "vtkm/CellTraits.h"
#include "vtkm/VecTraits.h"
#include "vtkm/VectorAnalysis.h"
#include "vtkm/exec/FunctorBase.h"
namespace vtkm
{
namespace worklet
{
namespace cellmetrics
{
// ========================= Unsupported cells ==================================
// By default, cells have zero shape unless the shape type template is specialized below.
template <typename OutType, typename PointCoordVecType, typename CellShapeType>
VTKM_EXEC OutType CellTaperMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
CellShapeType shape,
const vtkm::exec::FunctorBase& worklet)
{
UNUSED(numPts);
UNUSED(pts);
UNUSED(shape);
UNUSED(worklet);
return OutType(-1.0);
}
// ========================= 2D cells ==================================
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellTaperMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagQuad,
const vtkm::exec::FunctorBase& worklet)
{
UNUSED(numPts);
UNUSED(worklet);
using Scalar = OutType;
using CollectionOfPoints = PointCoordVecType;
using Vector = typename PointCoordVecType::ComponentType;
const Vector X12 = Vector((pts[0] - pts[1]) + (pts[2] - pts[3]));
const Vector X1 = GetQuadX0<Scalar, Vector, CollectionOfPoints>(pts);
const Vector X2 = GetQuadX1<Scalar, Vector, CollectionOfPoints>(pts);
const Scalar x12 = vtkm::Sqrt(vtkm::MagnitudeSquared(X12));
const Scalar x1 = vtkm::Sqrt(vtkm::MagnitudeSquared(X1));
const Scalar x2 = vtkm::Sqrt(vtkm::MagnitudeSquared(X2));
const Scalar minLength = vtkm::Min(x1, x2);
if (minLength <= Scalar(0.0))
{
return vtkm::Infinity<Scalar>();
}
const Scalar q = x12 / minLength;
return q;
}
// ========================= 3D cells ==================================
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellTaperMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagHexahedron,
const vtkm::exec::FunctorBase& worklet)
{
UNUSED(numPts);
UNUSED(worklet);
using Scalar = OutType;
Scalar X1 = vtkm::Sqrt(vtkm::MagnitudeSquared((pts[1] - pts[0]) + (pts[2] - pts[3]) +
(pts[5] - pts[4]) + (pts[6] - pts[7])));
Scalar X2 = vtkm::Sqrt(vtkm::MagnitudeSquared((pts[3] - pts[0]) + (pts[2] - pts[1]) +
(pts[7] - pts[4]) + (pts[6] - pts[5])));
Scalar X3 = vtkm::Sqrt(vtkm::MagnitudeSquared((pts[4] - pts[0]) + (pts[5] - pts[1]) +
(pts[6] - pts[2]) + (pts[7] - pts[3])));
if ((X1 <= Scalar(0.0)) || (X2 <= Scalar(0.0)) || (X3 <= Scalar(0.0)))
{
return vtkm::Infinity<Scalar>();
}
Scalar X12 = vtkm::Sqrt(vtkm::MagnitudeSquared(((pts[2] - pts[3]) - (pts[1] - pts[0])) +
((pts[6] - pts[7]) - (pts[5] - pts[4]))));
Scalar X13 = vtkm::Sqrt(vtkm::MagnitudeSquared(((pts[5] - pts[1]) - (pts[4] - pts[0])) +
((pts[6] - pts[2]) - (pts[7] - pts[3]))));
Scalar X23 = vtkm::Sqrt(vtkm::MagnitudeSquared(((pts[7] - pts[4]) - (pts[3] - pts[0])) +
((pts[6] - pts[5]) - (pts[2] - pts[1]))));
Scalar T12 = X12 / vtkm::Min(X1, X2);
Scalar T13 = X13 / vtkm::Min(X1, X3);
Scalar T23 = X23 / vtkm::Min(X2, X3);
return vtkm::Max(T12, vtkm::Max(T13, T23));
}
}
} // worklet
} // vtkm
#endif // vtk_m_worklet_CellTaper_Metric_h

80
vtkm/worklet/cellmetrics/CellWarpageMetric.h Normal file
View File

@ -0,0 +1,80 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//
// Copyright 2014 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
// Copyright 2014 UT-Battelle, LLC.
// Copyright 2014 Los Alamos National Security.
//
// Under the terms of Contract DE-NA0003525 with NTESS,
// the U.S. Government retains certain rights in this software.
//
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#ifndef vtk_m_worklet_CellWarpageMetric_h
#define vtk_m_worklet_CellWarpageMetric_h
#include "vtkm/CellShape.h"
#include "vtkm/CellTraits.h"
#include "vtkm/VecTraits.h"
#include "vtkm/VectorAnalysis.h"
#include "vtkm/exec/FunctorBase.h"
namespace vtkm
{
namespace worklet
{
namespace cellmetrics
{
template <typename OutType, typename PointCoordVecType, typename CellShapeType>
VTKM_EXEC OutType CellWarpageMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
CellShapeType shape,
const vtkm::exec::FunctorBase& worklet)
{
UNUSED(numPts);
UNUSED(pts);
UNUSED(shape);
UNUSED(worklet);
//worklet.RaiseError("Shape type template must be Quad to compute warpage");
return OutType(-1.0);
}
template <typename OutType, typename PointCoordVecType>
VTKM_EXEC OutType CellWarpageMetric(const vtkm::IdComponent& numPts,
const PointCoordVecType& pts,
vtkm::CellShapeTagQuad,
const vtkm::exec::FunctorBase& worklet)
{
UNUSED(numPts);
UNUSED(worklet);
using Scalar = OutType;
using CollectionOfPoints = PointCoordVecType;
using Vector = typename PointCoordVecType::ComponentType;
const Vector N0Mag = GetQuadN0Normalized<Scalar, Vector, CollectionOfPoints>(pts);
const Vector N1Mag = GetQuadN1Normalized<Scalar, Vector, CollectionOfPoints>(pts);
const Vector N2Mag = GetQuadN2Normalized<Scalar, Vector, CollectionOfPoints>(pts);
const Vector N3Mag = GetQuadN3Normalized<Scalar, Vector, CollectionOfPoints>(pts);
if (N0Mag < Scalar(0.0) || N1Mag < Scalar(0.0) || N2Mag < Scalar(0.0) || N3Mag < Scalar(0.0))
return vtkm::Infinity<Scalar>();
const Scalar n0dotn2 = vtkm::Dot(N0Mag, N2Mag);
const Scalar n1dotn3 = vtkm::Dot(N1Mag, N3Mag);
const Scalar min = vtkm::Min(n0dotn2, n1dotn3);
const Scalar minCubed = vtkm::Pow(min, 3);
//return Scalar(1.0) - minCubed; // AS DEFINED IN THE MANUAL
return minCubed; // AS DEFINED IN VISIT SOURCE CODE
}
}
} // worklet
} // vtkm
#endif // vtk_m_worklet_CellWarpage_Metric_h

4
vtkm/worklet/cellmetrics/TypeOfCellHexahedral.h
View File

@ -17,8 +17,8 @@
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#ifndef vtk_m_exec_cellmetrics_TypeOfCellHexahedral
#define vtk_m_exec_cellmetrics_TypeOfCellHexahedral
#ifndef vtk_m_worklet_cellmetrics_TypeOfCellHexahedral
#define vtk_m_worklet_cellmetrics_TypeOfCellHexahedral
/**
* The Verdict manual defines a set of commonly
* used components of a hexahedra (hex). For example,

4
vtkm/worklet/cellmetrics/TypeOfCellQuadrilateral.h
View File

@ -17,8 +17,8 @@
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#ifndef vtk_m_exec_cellmetrics_TypeOfCellQuadrilateral
#define vtk_m_exec_cellmetrics_TypeOfCellQuadrilateral
#ifndef vtk_m_worklet_cellmetrics_TypeOfCellQuadrilateral
#define vtk_m_worklet_cellmetrics_TypeOfCellQuadrilateral
/**
* The Verdict manual defines a set of commonly
* used components of a quadrilateral (quad). For example,

4
vtkm/worklet/cellmetrics/TypeOfCellTetrahedral.h
View File

@ -17,8 +17,8 @@
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#ifndef vtk_m_exec_cellmetrics_TypeOfCellTetrahedral
#define vtk_m_exec_cellmetrics_TypeOfCellTetrahedral
#ifndef vtk_m_worklet_cellmetrics_TypeOfCellTetrahedral
#define vtk_m_worklet_cellmetrics_TypeOfCellTetrahedral
/**
* Returns the L0 vector, as defined by the verdict manual.
*

4
vtkm/worklet/cellmetrics/TypeOfCellTriangle.h
View File

@ -17,8 +17,8 @@
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#ifndef vtk_m_exec_cellmetrics_TypeOfCellTriangle
#define vtk_m_exec_cellmetrics_TypeOfCellTriangle
#ifndef vtk_m_worklet_cellmetrics_TypeOfCellTriangle
#define vtk_m_worklet_cellmetrics_TypeOfCellTriangle
/**
* The Verdict manual defines a set of commonly
* used components of a triangle. For example,