VTK-m  2.0
CellRelativeSizeSquaredMetric.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_CellRelativeSizeSquaredMetric_h
21 #define vtk_m_worklet_cellmetrics_CellRelativeSizeSquaredMetric_h
22 
23 /*
24  * Mesh quality metric functions that compute the relative size squared of mesh
25  * cells. The RSS of a cell is defined as the square of the minimum of: the area
26  * divided by the average area of an ensemble of triangles or the inverse. For
27  * 3D cells we use the volumes instead of the areas.
28  *
29  * These metric computations are adapted from the VTK implementation of the
30  * Verdict library, which provides a set of mesh/cell metrics for evaluating the
31  * geometric qualities of regions of mesh spaces.
32  *
33  * See: The Verdict Library Reference Manual (for per-cell-type metric formulae)
34  * See: vtk/ThirdParty/verdict/vtkverdict (for VTK code implementation of this
35  * metric)
36  */
37 
38 #include <vtkm/CellShape.h>
39 #include <vtkm/CellTraits.h>
40 #include <vtkm/Matrix.h>
41 #include <vtkm/VecTraits.h>
42 #include <vtkm/VectorAnalysis.h>
43 #include <vtkm/exec/CellMeasure.h>
44 
45 #define UNUSED(expr) (void)(expr);
46 
47 namespace vtkm
48 {
49 namespace worklet
50 {
51 namespace cellmetrics
52 {
53 
55 
56 // ========================= Unsupported cells ==================================
57 
58 // By default, cells have zero shape unless the shape type template is specialized below.
59 template <typename OutType, typename PointCoordVecType, typename CellShapeType>
61  const PointCoordVecType& pts,
62  const OutType& avgArea,
63  CellShapeType shape,
65 {
66  UNUSED(numPts);
67  UNUSED(pts);
68  UNUSED(avgArea);
69  UNUSED(shape);
70  return OutType(-1.);
71 }
72 
73 // ========================= 2D cells ==================================
74 
75 template <typename OutType, typename PointCoordVecType>
77  const PointCoordVecType& pts,
78  const OutType& avgArea,
79  vtkm::CellShapeTagTriangle tag,
80  vtkm::ErrorCode& ec)
81 {
82  UNUSED(ec);
83  if (numPts != 3)
84  {
86  return OutType(-1.);
87  }
88  OutType A = vtkm::exec::CellMeasure<OutType>(numPts, pts, tag, ec);
89  OutType R = A / avgArea;
90  if (R == OutType(0.))
91  return OutType(0.);
92  OutType q = vtkm::Pow(vtkm::Min(R, OutType(1.) / R), OutType(2.));
93  return OutType(q);
94 }
95 
96 template <typename OutType, typename PointCoordVecType>
98  const PointCoordVecType& pts,
99  const OutType& avgArea,
100  vtkm::CellShapeTagQuad tag,
101  vtkm::ErrorCode& ec)
102 {
103  UNUSED(ec);
104  if (numPts != 4)
105  {
107  return OutType(-1.);
108  }
109  OutType A = vtkm::exec::CellMeasure<OutType>(numPts, pts, tag, ec);
110  OutType R = A / avgArea;
111  if (R == OutType(0.))
112  return OutType(0.);
113  OutType q = vtkm::Pow(vtkm::Min(R, OutType(1.) / R), OutType(2.));
114  return OutType(q);
115 }
116 
117 // ========================= 3D cells ==================================
118 
119 template <typename OutType, typename PointCoordVecType>
121  const PointCoordVecType& pts,
122  const OutType& avgVolume,
123  vtkm::CellShapeTagTetra tag,
124  vtkm::ErrorCode& ec)
125 {
126  UNUSED(ec);
127  if (numPts != 4)
128  {
130  return OutType(-1.);
131  }
132  OutType V = vtkm::exec::CellMeasure<OutType>(numPts, pts, tag, ec);
133  OutType R = V / avgVolume;
134  if (R == OutType(0.))
135  return OutType(0.);
136  OutType q = vtkm::Pow(vtkm::Min(R, OutType(1.) / R), OutType(2.));
137  return OutType(q);
138 }
139 
140 template <typename OutType, typename PointCoordVecType>
142  const PointCoordVecType& pts,
143  const OutType& avgVolume,
144  vtkm::CellShapeTagHexahedron tag,
145  vtkm::ErrorCode& ec)
146 {
147  UNUSED(tag);
148  UNUSED(ec);
149  if (numPts != 8)
150  {
152  return OutType(-1.);
153  }
154  OutType X1x = static_cast<OutType>((pts[1][0] - pts[0][0]) + (pts[2][0] - pts[3][0]) +
155  (pts[5][0] - pts[4][0]) + (pts[6][0] - pts[7][0]));
156  OutType X1y = static_cast<OutType>((pts[1][1] - pts[0][1]) + (pts[2][1] - pts[3][1]) +
157  (pts[5][1] - pts[4][1]) + (pts[6][1] - pts[7][1]));
158  OutType X1z = static_cast<OutType>((pts[1][2] - pts[0][2]) + (pts[2][2] - pts[3][2]) +
159  (pts[5][2] - pts[4][2]) + (pts[6][2] - pts[7][2]));
160 
161  OutType X2x = static_cast<OutType>((pts[2][0] - pts[0][0]) + (pts[2][0] - pts[1][0]) +
162  (pts[7][0] - pts[4][0]) + (pts[6][0] - pts[5][0]));
163  OutType X2y = static_cast<OutType>((pts[2][1] - pts[0][1]) + (pts[2][1] - pts[1][1]) +
164  (pts[7][1] - pts[4][1]) + (pts[6][1] - pts[5][1]));
165  OutType X2z = static_cast<OutType>((pts[2][2] - pts[0][2]) + (pts[2][2] - pts[1][2]) +
166  (pts[7][2] - pts[4][2]) + (pts[6][2] - pts[5][2]));
167 
168  OutType X3x = static_cast<OutType>((pts[4][0] - pts[0][0]) + (pts[5][0] - pts[1][0]) +
169  (pts[6][0] - pts[2][0]) + (pts[7][0] - pts[3][0]));
170  OutType X3y = static_cast<OutType>((pts[4][1] - pts[0][1]) + (pts[5][1] - pts[1][1]) +
171  (pts[6][1] - pts[2][1]) + (pts[7][1] - pts[3][1]));
172  OutType X3z = static_cast<OutType>((pts[4][2] - pts[0][2]) + (pts[5][2] - pts[1][2]) +
173  (pts[6][2] - pts[2][2]) + (pts[7][2] - pts[3][2]));
175  vtkm::MatrixSetRow(A8, 0, vtkm::Vec<OutType, 3>(X1x, X1y, X1z));
176  vtkm::MatrixSetRow(A8, 1, vtkm::Vec<OutType, 3>(X2x, X2y, X2z));
177  vtkm::MatrixSetRow(A8, 2, vtkm::Vec<OutType, 3>(X3x, X3y, X3z));
178  OutType D = vtkm::MatrixDeterminant(A8);
179  D = D / (OutType(64.) * avgVolume);
180  if (D == OutType(0.))
181  return OutType(0.);
182  OutType q = vtkm::Pow(vtkm::Min(D, OutType(1.) / D), OutType(2.));
183  return OutType(q);
184 }
185 
186 } // namespace cellmetrics
187 } // namespace worklet
188 } // namespace vtkm
189 
190 #endif // vtk_m_worklet_cellmetrics_CellRelativeSizeSquaredMetric_h
vtkm::ErrorCode
ErrorCode
Definition: ErrorCode.h:19
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: CellRelativeSizeSquaredMetric.h:45
vtkm::IdComponent
vtkm::Int32 IdComponent
Represents a component ID (index of component in a vector).
Definition: Types.h:168
vtkm::MatrixSetRow
VTKM_EXEC_CONT void MatrixSetRow(vtkm::Matrix< T, NumRow, NumCol > &matrix, vtkm::IdComponent rowIndex, const vtkm::Vec< T, NumCol > &rowValues)
Convenience function for setting a row of a matrix.
Definition: Matrix.h:130
Matrix.h
vtkm::MatrixDeterminant
VTKM_EXEC_CONT T MatrixDeterminant(const vtkm::Matrix< T, Size, Size > &A)
Compute the determinant of a matrix.
Definition: Matrix.h:491
vtkm::worklet::cellmetrics::FloatType
vtkm::FloatDefault FloatType
Definition: CellAspectFrobeniusMetric.h:50
CellShape.h
CellMeasure.h
VectorAnalysis.h
vtkm::Vec
A short fixed-length array.
Definition: Types.h:767
vtkm::Matrix
Basic Matrix type.
Definition: Matrix.h:33
vtkm::FloatDefault
vtkm::Float32 FloatDefault
The floating point type to use when no other precision is specified.
Definition: Types.h:198
CellTraits.h
vtkm::worklet::cellmetrics::CellRelativeSizeSquaredMetric
VTKM_EXEC OutType CellRelativeSizeSquaredMetric(const vtkm::IdComponent &numPts, const PointCoordVecType &pts, const OutType &avgArea, CellShapeType shape, vtkm::ErrorCode &)
Definition: CellRelativeSizeSquaredMetric.h:60
VecTraits.h
vtkm::ErrorCode::InvalidNumberOfPoints
@ InvalidNumberOfPoints