VTK-m  2.0
CellAspectFrobeniusMetric.h
Go to the documentation of this file.
1 //============================================================================
2 // Copyright (c) Kitware, Inc.
3 // All rights reserved.
4 // See LICENSE.txt for details.
5 // This software is distributed WITHOUT ANY WARRANTY; without even
6 // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
7 // PURPOSE. See the above copyright notice for more information.
8 //
9 // Copyright 2014 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
10 // Copyright 2014 UT-Battelle, LLC.
11 // Copyright 2014 Los Alamos National Security.
12 //
13 // Under the terms of Contract DE-NA0003525 with NTESS,
14 // the U.S. Government retains certain rights in this software.
15 //
16 // Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
17 // Laboratory (LANL), the U.S. Government retains certain rights in
18 // this software.
19 //============================================================================
20 #ifndef vtk_m_worklet_cellmetrics_CellAspectFrobeniusMetric_h
21 #define vtk_m_worklet_cellmetrics_CellAspectFrobeniusMetric_h
22 
23 /*
24 * Mesh quality metric functions that compute the aspect frobenius of certain mesh cells.
25 * The aspect frobenius metric generally measures the degree of regularity of a cell, with
26 * a value of 1 representing a regular cell..
27 *
28 * These metric computations are adapted from the VTK implementation of the Verdict library,
29 * which provides a set of mesh/cell metrics for evaluating the geometric qualities of regions
30 * of mesh spaces.
31 *
32 * See: The Verdict Library Reference Manual (for per-cell-type metric formulae)
33 * See: vtk/ThirdParty/verdict/vtkverdict (for VTK code implementation of this metric)
34 */
35 
36 #include <vtkm/CellShape.h>
37 #include <vtkm/CellTraits.h>
38 #include <vtkm/ErrorCode.h>
39 #include <vtkm/VecTraits.h>
40 #include <vtkm/VectorAnalysis.h>
41 
42 #define UNUSED(expr) (void)(expr);
43 
44 namespace vtkm
45 {
46 namespace worklet
47 {
48 namespace cellmetrics
49 {
51 
52 // ========================= Unsupported cells ==================================
53 
54 // By default, cells have undefined aspect frobenius unless the shape type template is specialized below.
55 template <typename OutType, typename PointCoordVecType, typename CellShapeType>
57  const PointCoordVecType& pts,
58  CellShapeType shape,
59  vtkm::ErrorCode& ec)
60 {
61  UNUSED(numPts);
62  UNUSED(pts);
63  UNUSED(shape);
65  return OutType(0.0);
66 }
67 
68 //If the polygon has 3 vertices or 4 vertices, then just call
69 //the functions for Triangle and Quad cell types. Otherwise,
70 //this metric is not supported for (n>4)-vertex polygons, such
71 //as pentagons or hexagons, or (n<3)-vertex polygons, such as lines or points.
72 template <typename OutType, typename PointCoordVecType>
74  const PointCoordVecType& pts,
75  vtkm::CellShapeTagPolygon,
76  vtkm::ErrorCode& ec)
77 {
78  if (numPts == 3)
79  return CellAspectFrobeniusMetric<OutType>(numPts, pts, vtkm::CellShapeTagTriangle(), ec);
80  else
81  {
83  return OutType(0.0);
84  }
85 }
86 
87 //The aspect frobenius metric is not supported for lines/edges.
88 template <typename OutType, typename PointCoordVecType>
90  const PointCoordVecType& pts,
91  vtkm::CellShapeTagLine,
92  vtkm::ErrorCode& ec)
93 {
94  UNUSED(numPts);
95  UNUSED(pts);
97  return OutType(0.0);
98 }
99 
100 //The aspect frobenius metric is not uniquely defined for quads.
101 //Instead, use either the mean or max aspect frobenius metrics, which are
102 //defined in terms of the aspect frobenius of triangles.
103 template <typename OutType, typename PointCoordVecType>
105  const PointCoordVecType& pts,
106  vtkm::CellShapeTagQuad,
107  vtkm::ErrorCode& ec)
108 {
109  UNUSED(numPts);
110  UNUSED(pts);
112  return OutType(0.0);
113 }
114 
115 //The aspect frobenius metric is not uniquely defined for hexahedrons.
116 //Instead, use either the mean or max aspect frobenius metrics, which are
117 //defined in terms of the aspect frobenius of tetrahedrons.
118 template <typename OutType, typename PointCoordVecType>
120  const PointCoordVecType& pts,
121  vtkm::CellShapeTagHexahedron,
122  vtkm::ErrorCode& ec)
123 {
124  UNUSED(numPts);
125  UNUSED(pts);
127  return OutType(0.0);
128 }
129 
130 //The aspect frobenius metric is not uniquely defined for pyramids.
131 //Instead, use either the mean or max aspect frobenius metrics, which are
132 //defined in terms of the aspect frobenius of tetrahedrons.
133 template <typename OutType, typename PointCoordVecType>
135  const PointCoordVecType& pts,
136  vtkm::CellShapeTagPyramid,
137  vtkm::ErrorCode& ec)
138 {
139  UNUSED(numPts);
140  UNUSED(pts);
142  return OutType(0.0);
143 }
144 
145 //The aspect frobenius metric is not uniquely defined for wedges.
146 //Instead, use either the mean or max aspect frobenius metrics, which are
147 //defined in terms of the aspect frobenius of tetrahedrons.
148 template <typename OutType, typename PointCoordVecType>
150  const PointCoordVecType& pts,
151  vtkm::CellShapeTagWedge,
152  vtkm::ErrorCode& ec)
153 {
154  UNUSED(numPts);
155  UNUSED(pts);
157  return OutType(0.0);
158 }
159 
160 // ========================= 2D cells ==================================
161 
162 // Computes the aspect frobenius of a triangle.
163 // Formula: Sum of lengths of 3 edges, divided by a multiple of the triangle area.
164 // Equals 1 for an equilateral unit triangle.
165 // Acceptable range: [1,1.3]
166 // Full range: [1,FLOAT_MAX]
167 template <typename OutType, typename PointCoordVecType>
169  const PointCoordVecType& pts,
170  vtkm::CellShapeTagTriangle,
171  vtkm::ErrorCode& ec)
172 {
173  if (numPts != 3)
174  {
176  return OutType(0.0);
177  }
178 
179  //The 3 edges of a triangle
180  using Edge = typename PointCoordVecType::ComponentType;
181  const Edge TriEdges[3] = { pts[1] - pts[0], pts[2] - pts[1], pts[0] - pts[2] };
182 
183  //Sum the length squared of each edge
184  FloatType sum = (FloatType)vtkm::MagnitudeSquared(TriEdges[0]) +
186 
187  //Compute the length of the cross product of the triangle.
188  //The result is twice the area of the triangle.
189  FloatType crossLen = (FloatType)vtkm::Magnitude(vtkm::Cross(TriEdges[0], -TriEdges[2]));
190 
191  if (crossLen == 0.0)
192  return vtkm::Infinity<OutType>();
193 
194  OutType aspect_frobenius = (OutType)(sum / (vtkm::Sqrt(3.0) * 2 * crossLen));
195 
196  if (aspect_frobenius > 0.0)
197  return vtkm::Min(aspect_frobenius, vtkm::Infinity<OutType>());
198 
199  return vtkm::Max(aspect_frobenius, vtkm::NegativeInfinity<OutType>());
200 } // ============================= 3D Volume cells ==================================i
201 
202 // Computes the aspect frobenius of a tetrahedron.
203 // Formula: Sum of lengths of 3 edges, divided by a multiple of the triangle area.
204 // Equals 1 for a right regular tetrahedron (4 equilateral triangles).
205 // Acceptable range: [1,1.3]
206 // Full range: [1,FLOAT_MAX]
207 template <typename OutType, typename PointCoordVecType>
209  const PointCoordVecType& pts,
210  vtkm::CellShapeTagTetra,
211  vtkm::ErrorCode& ec)
212 {
213  if (numPts != 4)
214  {
216  return OutType(0.0);
217  }
218 
219  //Two base edges and one vertical edge, used to compute the tet volume
220  using Edge = typename PointCoordVecType::ComponentType;
221 
222  const Edge TetEdges[3] = {
223  pts[1] - pts[0], //Base edge 1
224  pts[2] - pts[0], //Base edge 2
225  pts[3] - pts[0] //Vert edge 3
226  };
227 
228  //Compute the tet volume
229  FloatType denominator = (FloatType)vtkm::Dot(TetEdges[0], vtkm::Cross(TetEdges[1], TetEdges[2]));
230  denominator *= denominator;
231  denominator *= 2.0f;
232  const FloatType normal_exp = 1.0f / 3.0f;
233  denominator = 3.0f * vtkm::Pow(denominator, normal_exp);
234 
235  if (denominator < vtkm::NegativeInfinity<FloatType>())
236  return vtkm::Infinity<OutType>();
237 
238  FloatType numerator = (FloatType)vtkm::Dot(TetEdges[0], TetEdges[0]);
239  numerator += (FloatType)vtkm::Dot(TetEdges[1], TetEdges[1]);
240  numerator += (FloatType)vtkm::Dot(TetEdges[2], TetEdges[2]);
241  numerator *= 1.5f;
242  numerator -= (FloatType)vtkm::Dot(TetEdges[0], TetEdges[1]);
243  numerator -= (FloatType)vtkm::Dot(TetEdges[0], TetEdges[2]);
244  numerator -= (FloatType)vtkm::Dot(TetEdges[1], TetEdges[2]);
245 
246  OutType aspect_frobenius = (OutType)(numerator / denominator);
247 
248  if (aspect_frobenius > 0.0)
249  return vtkm::Min(aspect_frobenius, vtkm::Infinity<OutType>());
250 
251  return vtkm::Max(aspect_frobenius, vtkm::NegativeInfinity<OutType>());
252 }
253 } // namespace cellmetrics
254 } // namespace worklet
255 } // namespace vtkm
256 #endif // vtk_m_worklet_cellmetrics_CellAspectFrobeniusMetric_h
vtkm::Sqrt
VTKM_EXEC_CONT vtkm::Float32 Sqrt(vtkm::Float32 x)
Compute the square root of x.
Definition: Math.h:958
vtkm::ErrorCode
ErrorCode
Definition: ErrorCode.h:19
vtkm::MagnitudeSquared
VTKM_EXEC_CONT detail::FloatingPointReturnType< T >::Type MagnitudeSquared(const T &x)
Returns the square of the magnitude of a vector.
Definition: VectorAnalysis.h:64
VTKM_EXEC
#define VTKM_EXEC
Definition: ExportMacros.h:51
vtkm
Groups connected points that have the same field value.
Definition: Atomic.h:19
UNUSED
#define UNUSED(expr)
Definition: CellAspectFrobeniusMetric.h:42
vtkm::IdComponent
vtkm::Int32 IdComponent
Represents a component ID (index of component in a vector).
Definition: Types.h:168
vtkm::Magnitude
VTKM_EXEC_CONT detail::FloatingPointReturnType< T >::Type Magnitude(const T &x)
Returns the magnitude of a vector.
Definition: VectorAnalysis.h:100
vtkm::worklet::cellmetrics::FloatType
vtkm::FloatDefault FloatType
Definition: CellAspectFrobeniusMetric.h:50
CellShape.h
vtkm::worklet::cellmetrics::CellAspectFrobeniusMetric
VTKM_EXEC OutType CellAspectFrobeniusMetric(const vtkm::IdComponent &numPts, const PointCoordVecType &pts, CellShapeType shape, vtkm::ErrorCode &ec)
Definition: CellAspectFrobeniusMetric.h:56
ErrorCode.h
VectorAnalysis.h
vtkm::ErrorCode::InvalidCellMetric
@ InvalidCellMetric
vtkm::Cross
VTKM_EXEC_CONT vtkm::Vec< typename detail::FloatingPointReturnType< T >::Type, 3 > Cross(const vtkm::Vec< T, 3 > &x, const vtkm::Vec< T, 3 > &y)
Find the cross product of two vectors.
Definition: VectorAnalysis.h:177
vtkm::FloatDefault
vtkm::Float32 FloatDefault
The floating point type to use when no other precision is specified.
Definition: Types.h:198
CellTraits.h
VecTraits.h
vtkm::ErrorCode::InvalidNumberOfPoints
@ InvalidNumberOfPoints