384 lines
16 KiB
C++
384 lines
16 KiB
C++
//============================================================================
|
|
// 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 = static_cast<OutType>(
|
|
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 =
|
|
static_cast<OutType>(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,
|
|
vtkm::ErrorCode& ec)
|
|
{
|
|
UNUSED(numPts);
|
|
UNUSED(pts);
|
|
UNUSED(shape);
|
|
ec = vtkm::ErrorCode::InvalidCellMetric;
|
|
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,
|
|
vtkm::ErrorCode& ec)
|
|
{
|
|
if (numPts == 3)
|
|
return vtkm::worklet::cellmetrics::CellAspectFrobeniusMetric<OutType>(
|
|
numPts, pts, vtkm::CellShapeTagTriangle(), ec);
|
|
if (numPts == 4)
|
|
return CellMaxAspectFrobeniusMetric<OutType>(numPts, pts, vtkm::CellShapeTagQuad(), ec);
|
|
else
|
|
{
|
|
ec = vtkm::ErrorCode::InvalidCellMetric;
|
|
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,
|
|
vtkm::ErrorCode& ec)
|
|
{
|
|
UNUSED(numPts);
|
|
UNUSED(pts);
|
|
ec = vtkm::ErrorCode::InvalidCellMetric;
|
|
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,
|
|
vtkm::ErrorCode& ec)
|
|
{
|
|
if (numPts != 3)
|
|
{
|
|
ec = vtkm::ErrorCode::InvalidNumberOfPoints;
|
|
return OutType(0.0);
|
|
}
|
|
return vtkm::worklet::cellmetrics::CellAspectFrobeniusMetric<OutType>(
|
|
numPts, pts, vtkm::CellShapeTagTriangle(), ec);
|
|
}
|
|
|
|
//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,
|
|
vtkm::ErrorCode& ec)
|
|
{
|
|
UNUSED(numPts);
|
|
UNUSED(pts);
|
|
ec = vtkm::ErrorCode::InvalidCellMetric;
|
|
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,
|
|
vtkm::ErrorCode& ec)
|
|
{
|
|
if (numPts != 4)
|
|
{
|
|
ec = vtkm::ErrorCode::InvalidNumberOfPoints;
|
|
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,
|
|
vtkm::ErrorCode& ec)
|
|
{
|
|
if (numPts != 4)
|
|
{
|
|
ec = vtkm::ErrorCode::InvalidNumberOfPoints;
|
|
return OutType(0.0);
|
|
}
|
|
|
|
return vtkm::worklet::cellmetrics::CellAspectFrobeniusMetric<OutType, PointCoordVecType>(
|
|
numPts, pts, vtkm::CellShapeTagTetra(), ec);
|
|
}
|
|
|
|
// 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,
|
|
vtkm::ErrorCode& ec)
|
|
{
|
|
if (numPts != 8)
|
|
{
|
|
ec = vtkm::ErrorCode::InvalidNumberOfPoints;
|
|
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 *= (OutType)0.3333333;
|
|
|
|
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,
|
|
vtkm::ErrorCode& ec)
|
|
{
|
|
if (numPts != 6)
|
|
{
|
|
ec = vtkm::ErrorCode::InvalidNumberOfPoints;
|
|
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(), ec);
|
|
OutType curr;
|
|
for (vtkm::Id i = 1; i < 6; i++)
|
|
{
|
|
curr = vtkm::worklet::cellmetrics::CellAspectFrobeniusMetric<OutType>(
|
|
tetPts, tetras[i], vtkm::CellShapeTagTetra(), ec);
|
|
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
|