VTK-m  2.0
CellMaxAspectFrobeniusMetric.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_CellMaxAspectFrobeniusMetric_h
21 #define vtk_m_worklet_cellmetrics_CellMaxAspectFrobeniusMetric_h
22 
23 /*
24 * Mesh quality metric functions that compute the maximum aspect frobenius of certain mesh cells,
25 * each of which are composed of two or more triangles or tetrahedrons. The output aspect metric
26 * value is the maximum among all triangles or tetrahedrons.
27 * The aspect frobenius metric generally measures the degree of regularity of a cell, with
28 * a value of 1 representing a regular cell.
29 ** These metric computations are adapted from the VTK implementation of the Verdict library,
30 * which provides a set of mesh/cell metrics for evaluating the geometric qualities of regions
31 * of mesh spaces.
32 ** The maximum aspect frobenious computations for pyramid and wedge cell types are not defined in the
33 * VTK implementation, but are provided here.
34 ** See: The Verdict Library Reference Manual (for per-cell-type metric formulae)
35 * See: vtk/ThirdParty/verdict/vtkverdict (for VTK code implementation of this metric)
36 */
37 
39 #include <vtkm/CellShape.h>
40 #include <vtkm/CellTraits.h>
41 #include <vtkm/ErrorCode.h>
42 #include <vtkm/VecTraits.h>
43 #include <vtkm/VectorAnalysis.h>
44 
45 #define UNUSED(expr) (void)(expr);
46 
47 namespace vtkm
48 {
49 namespace worklet
50 {
51 namespace cellmetrics
52 {
53 
55 
56 //This approximates the aspect frobenius of a tetrahedron, except for slight
57 //mathematical differences. In the standard aspect frobenius metric, a tetrahedron
58 //is compared to a reference right equilateral tetrahedron. However, in the max
59 //aspect frobenius metric of hexahedrons, the component tetrahedrons are compared
60 //to reference right isoceles tetrahedrons. Thus, some of the calculations differ
61 //to account for the change in reference tetrahedron. This condition computation
62 //is not to be confused with the separate CellConditionMetric metric, but is similar
63 //in computation.
64 template <typename OutType, typename VecType>
65 VTKM_EXEC OutType ComputeTetCondition(const VecType edges[])
66 {
67  //Compute the determinant/volume of the reference tet.
68  //(right isosceles tet for max aspect frobenius of hexs, pyramids, and wedges)
69  OutType det = (OutType)vtkm::Dot(edges[0], vtkm::Cross(edges[1], edges[2]));
70 
71  if (det <= vtkm::NegativeInfinity<OutType>())
72  return vtkm::Infinity<OutType>();
73 
74  OutType term1 = static_cast<OutType>(
75  vtkm::Dot(edges[0], edges[0]) + vtkm::Dot(edges[1], edges[1]) + vtkm::Dot(edges[2], edges[2]));
76 
77  VecType crosses[3] = { vtkm::Cross(edges[0], edges[1]),
78  vtkm::Cross(edges[1], edges[2]),
79  vtkm::Cross(edges[2], edges[0]) };
80 
81  OutType term2 =
82  static_cast<OutType>(vtkm::Dot(crosses[0], crosses[0]) + vtkm::Dot(crosses[1], crosses[1]) +
83  vtkm::Dot(crosses[2], crosses[2]));
84 
85  return vtkm::Sqrt(term1 * term2) / det;
86 }
87 
88 // ========================= Unsupported cells ==================================
89 
90 // By default, cells have undefined aspect frobenius unless the shape type template is specialized below.
91 template <typename OutType, typename PointCoordVecType, typename CellShapeType>
93  const PointCoordVecType&,
94  CellShapeType,
95  vtkm::ErrorCode& ec)
96 {
98  return OutType(0.0);
99 }
100 
101 //If the polygon has 3 vertices or 4 vertices, then just call
102 //the functions for Triangle and Quad cell types. Otherwise,
103 //this metric is not supported for (n>4)-vertex polygons, such
104 //as pentagons or hexagons, or (n<3)-vertex polygons, such as lines or points.
105 template <typename OutType, typename PointCoordVecType>
107  const PointCoordVecType& pts,
108  vtkm::CellShapeTagPolygon,
109  vtkm::ErrorCode& ec)
110 {
111  if (numPts == 3)
112  return vtkm::worklet::cellmetrics::CellAspectFrobeniusMetric<OutType>(
113  numPts, pts, vtkm::CellShapeTagTriangle(), ec);
114  if (numPts == 4)
115  return CellMaxAspectFrobeniusMetric<OutType>(numPts, pts, vtkm::CellShapeTagQuad(), ec);
116  else
117  {
119  return OutType(0.0);
120  }
121 }
122 
123 //The max aspect frobenius metric is not supported for lines/edges.
124 template <typename OutType, typename PointCoordVecType>
126  const PointCoordVecType&,
127  vtkm::CellShapeTagLine,
128  vtkm::ErrorCode& ec)
129 {
131  return OutType(0.0);
132 }
133 
134 //The max aspect frobenius metric is not uniquely defined for triangles,
135 //since the standard aspect frobenius metric is used for triangles.
136 //Thus, this implementation simply calls the aspect frobenius metric.
137 template <typename OutType, typename PointCoordVecType>
139  const PointCoordVecType& pts,
140  vtkm::CellShapeTagTriangle,
141  vtkm::ErrorCode& ec)
142 {
143  if (numPts != 3)
144  {
146  return OutType(0.0);
147  }
148  return vtkm::worklet::cellmetrics::CellAspectFrobeniusMetric<OutType>(
149  numPts, pts, vtkm::CellShapeTagTriangle(), ec);
150 }
151 
152 //The max aspect frobenius metric is not supported for pyramids.
153 template <typename OutType, typename PointCoordVecType>
155  const PointCoordVecType&,
156  vtkm::CellShapeTagPyramid,
157  vtkm::ErrorCode& ec)
158 {
160  return OutType(0.0);
161 }
162 
163 // ========================= 2D cells ==================================
164 
165 // Computes the max aspect frobenius of a quad.
166 // Formula: The maximum aspect frobenius metric among the four triangles formed
167 // at the four corner points of the quad. Given a corner point, two other points are
168 // chosen in a uniform, counter-clockwise manner to form a triangle. The aspect frobenius metric
169 // is computed on this triangle. The maximum among this four computed triangle metrics
170 // is returned as output.
171 // Equals 1 for a unit square.
172 // Acceptable range: [1,1.3]
173 // Full range: [1,FLOAT_MAX]
174 template <typename OutType, typename PointCoordVecType>
176  const PointCoordVecType& pts,
177  vtkm::CellShapeTagQuad,
178  vtkm::ErrorCode& ec)
179 {
180  if (numPts != 4)
181  {
183  return OutType(0.0);
184  }
185  //The 4 edges of a quad
186  using Edge = typename PointCoordVecType::ComponentType;
187  const Edge QuadEdges[4] = { pts[1] - pts[0], pts[2] - pts[1], pts[3] - pts[2], pts[0] - pts[3] };
188 
189  FloatType a2 = vtkm::MagnitudeSquared(QuadEdges[0]);
190  FloatType b2 = vtkm::MagnitudeSquared(QuadEdges[1]);
191  FloatType c2 = vtkm::MagnitudeSquared(QuadEdges[2]);
192  FloatType d2 = vtkm::MagnitudeSquared(QuadEdges[3]);
193 
194  //Compute the length of the cross product for each of the 4 reference triangles.
195  //The result is twice the area of the triangle.
196  FloatType ab = vtkm::Magnitude(vtkm::Cross(QuadEdges[0], QuadEdges[1]));
197  FloatType bc = vtkm::Magnitude(vtkm::Cross(QuadEdges[1], QuadEdges[2]));
198  FloatType cd = vtkm::Magnitude(vtkm::Cross(QuadEdges[2], QuadEdges[3]));
199  FloatType da = vtkm::Magnitude(vtkm::Cross(QuadEdges[3], QuadEdges[0]));
200 
201  if (ab < vtkm::NegativeInfinity<OutType>() || bc < vtkm::NegativeInfinity<OutType>() ||
202  cd < vtkm::NegativeInfinity<OutType>() || da < vtkm::NegativeInfinity<OutType>())
203  return vtkm::Infinity<OutType>();
204 
205  FloatType qmax = (a2 + b2) / ab; //Initial max aspect frobenius (triangle 0)
206 
207  FloatType qcur = (b2 + c2) / bc; //Keep checking/updating max (triangles 1 - 3)
208  qmax = qmax > qcur ? qmax : qcur;
209 
210  qcur = (c2 + d2) / cd;
211  qmax = qmax > qcur ? qmax : qcur;
212 
213  qcur = (d2 + a2) / da;
214  qmax = qmax > qcur ? qmax : qcur;
215 
216  OutType max_aspect_frobenius = 0.5f * (OutType)qmax;
217 
218  if (max_aspect_frobenius > 0)
219  vtkm::Min(max_aspect_frobenius, vtkm::Infinity<OutType>());
220 
221  return vtkm::Max(max_aspect_frobenius, vtkm::NegativeInfinity<OutType>());
222 }
223 
224 // ============================= 3D Volume cells ==================================i
225 
226 // Computes the aspect frobenius of a tetrahedron.
227 // Formula: Sum of lengths of 3 edges, divided by a multiple of the triangle area.
228 // Equals 1 for a right regular tetrahedron (4 equilateral triangles).
229 // Acceptable range: [1,1.3]
230 // Full range: [1,FLOAT_MAX]
231 template <typename OutType, typename PointCoordVecType>
233  const PointCoordVecType& pts,
234  vtkm::CellShapeTagTetra,
235  vtkm::ErrorCode& ec)
236 {
237  if (numPts != 4)
238  {
240  return OutType(0.0);
241  }
242 
243  return vtkm::worklet::cellmetrics::CellAspectFrobeniusMetric<OutType, PointCoordVecType>(
244  numPts, pts, vtkm::CellShapeTagTetra(), ec);
245 }
246 
247 // Computes the maximum aspect frobenius of a hexahedron.
248 // Formula: The maximum aspect frobenius metric among the eight tetrahedrons formed
249 // at the eight corner points of the hex. Given a corner point, three other points are
250 // chosen in a uniform, counter-clockwise manner to form a tetrahedron. The aspect frobenius metric
251 // is computed on this tet, with respect to a reference right isosceles tet. The maximum among
252 // these eight computed tet metrics is returned as output.
253 // Equals 1 for a unit cube (right isosceles tet formed at all 8 corner points).
254 // Acceptable range: [1,1.3]
255 // Full range: [1,FLOAT_MAX]
256 template <typename OutType, typename PointCoordVecType>
258  const PointCoordVecType& pts,
259  vtkm::CellShapeTagHexahedron,
260  vtkm::ErrorCode& ec)
261 {
262  if (numPts != 8)
263  {
265  return OutType(0.0);
266  }
267 
268  using Edge = typename PointCoordVecType::ComponentType;
269 
270  //8 tets: one constructed at each different corner of the hex
271  //For each tet: Two base edges and one vertical edge, used to compute the tet volume
272  const Edge TetEdges[8][3] = {
273  {
274  pts[1] - pts[0], //Base edge 1
275  pts[3] - pts[0], //Base edge 2
276  pts[4] - pts[0] //Vertical edge 3
277  }, //tet 0
278 
279  { pts[2] - pts[1], pts[0] - pts[1], pts[5] - pts[1] }, //tet 1
280 
281  { pts[3] - pts[2], pts[1] - pts[2], pts[6] - pts[2] }, //tet 2
282 
283  { pts[0] - pts[3], pts[2] - pts[3], pts[7] - pts[3] }, //tet 3
284 
285  { pts[7] - pts[4], pts[5] - pts[4], pts[0] - pts[4] }, //tet 4
286 
287  { pts[4] - pts[5], pts[6] - pts[5], pts[1] - pts[5] }, //tet 5
288 
289  { pts[5] - pts[6], pts[7] - pts[6], pts[2] - pts[6] }, //tet 6
290 
291  { pts[6] - pts[7], pts[4] - pts[7], pts[3] - pts[7] } //tet 7
292  };
293 
294  //For each tet, compute the condition metric, which approximates the deviation of the
295  //tet's volume to that of a right isoceles tetrahedron. Return the maximum metric value
296  //among all 8 tets as the maximum aspect frobenius.
297  OutType max_aspect_frobenius = ComputeTetCondition<OutType, Edge>(TetEdges[0]);
298  OutType curr;
299  for (vtkm::Id i = 1; i < 8; i++)
300  {
301  curr = ComputeTetCondition<OutType, Edge>(TetEdges[i]);
302  if (curr <= 0.0)
303  return vtkm::Infinity<OutType>();
304  if (curr > max_aspect_frobenius)
305  max_aspect_frobenius = curr;
306  }
307 
308  max_aspect_frobenius *= (OutType)0.3333333;
309 
310  if (max_aspect_frobenius > 0)
311  return vtkm::Min(max_aspect_frobenius, vtkm::Infinity<OutType>());
312 
313  return vtkm::Max(max_aspect_frobenius, OutType(0.0));
314 }
315 
316 // Computes the maximum aspect frobenius of a wedge.
317 // Formula: The maximum aspect frobenius metric among the six tetrahedrons formed
318 // from the six corner points of the two triangular faces. Given a corner point, three other points are
319 // chosen in a uniform, counter-clockwise manner to form a tetrahedron. The aspect frobenius metric
320 // is computed on this tet, with respect to an equilateral tet. The maximum among
321 // these six computed tet metrics is returned as output.
322 // Equals 1 for a unit wedge (two equilateral triangles of unit edge length and 3 unit squares).
323 // Acceptable range: [1,1.3]
324 // Full range: [1,FLOAT_MAX]
325 template <typename OutType, typename PointCoordVecType>
327  const PointCoordVecType& pts,
328  vtkm::CellShapeTagWedge,
329  vtkm::ErrorCode& ec)
330 {
331  if (numPts != 6)
332  {
334  return OutType(0.0);
335  }
336 
337  using Vector = typename PointCoordVecType::ComponentType;
338 
339  //Four positively-oriented points of each tet
340  const vtkm::Vec<Vector, 4> tetras[6] = {
341  vtkm::Vec<Vector, 4>(pts[0], pts[1], pts[2], pts[3]), //tet 0
342  vtkm::Vec<Vector, 4>(pts[1], pts[2], pts[0], pts[4]), //tet 1
343  vtkm::Vec<Vector, 4>(pts[2], pts[0], pts[1], pts[5]), //tet 2
344  vtkm::Vec<Vector, 4>(pts[3], pts[5], pts[4], pts[0]), //tet 3
345  vtkm::Vec<Vector, 4>(pts[4], pts[3], pts[5], pts[1]), //tet 4
346  vtkm::Vec<Vector, 4>(pts[5], pts[4], pts[3], pts[2]) //tet 5
347  };
348 
349  //For each tet, call the aspect frobenius metric.
350  //Return the maximum metric value among all 6 tets as the maximum aspect frobenius.
351  const vtkm::IdComponent tetPts = 4;
352  OutType max_aspect_frobenius = vtkm::worklet::cellmetrics::CellAspectFrobeniusMetric<OutType>(
353  tetPts, tetras[0], vtkm::CellShapeTagTetra(), ec);
354  OutType curr;
355  for (vtkm::Id i = 1; i < 6; i++)
356  {
357  curr = vtkm::worklet::cellmetrics::CellAspectFrobeniusMetric<OutType>(
358  tetPts, tetras[i], vtkm::CellShapeTagTetra(), ec);
359  if (curr > max_aspect_frobenius)
360  max_aspect_frobenius = curr;
361  }
362 
363  //Divide by metric value of unit wedge (normalization)
364  max_aspect_frobenius /= 1.16477;
365 
366  if (max_aspect_frobenius > 0)
367  {
368  return vtkm::Min(max_aspect_frobenius, vtkm::Infinity<OutType>());
369  }
370  else
371  {
372  return vtkm::Max(max_aspect_frobenius, vtkm::NegativeInfinity<OutType>());
373  }
374 }
375 } // namespace cellmetrics
376 } // namespace worklet
377 } // namespace vtkm
378 #endif // vtk_m_worklet_cellmetrics_CellMaxAspectFrobeniusMetric_h
vtkm::worklet::cellmetrics::CellMaxAspectFrobeniusMetric
VTKM_EXEC OutType CellMaxAspectFrobeniusMetric(const vtkm::IdComponent &, const PointCoordVecType &, CellShapeType, vtkm::ErrorCode &ec)
Definition: CellMaxAspectFrobeniusMetric.h:92
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
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::Id
vtkm::Int32 Id
Represents an ID (index into arrays).
Definition: Types.h:191
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::worklet::cellmetrics::ComputeTetCondition
VTKM_EXEC OutType ComputeTetCondition(const VecType edges[])
Definition: CellMaxAspectFrobeniusMetric.h:65
CellAspectFrobeniusMetric.h
vtkm::Vec
A short fixed-length array.
Definition: Types.h:767
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