VTK-m  2.0
CellShapeMetric.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_CellShapeMetric_h
21 #define vtk_m_worklet_CellShapeMetric_h
22 /*
23 * Mesh quality metric functions that compute the shape, or weighted Jacobian, of mesh cells.
24 * The Jacobian of a cell is weighted by the condition metric value of the cell.
25 ** These metric computations are adapted from the VTK implementation of the Verdict library,
26 * which provides a set of cell metrics for evaluating the geometric qualities of regions of mesh spaces.
27 ** See: The Verdict Library Reference Manual (for per-cell-type metric formulae)
28 * See: vtk/ThirdParty/verdict/vtkverdict (for VTK code implementation of this metric)
29 */
30 
31 #include "CellConditionMetric.h"
32 #include "CellJacobianMetric.h"
33 #include "TypeOfCellHexahedral.h"
34 #include "TypeOfCellQuadrilateral.h"
35 #include "TypeOfCellTetrahedral.h"
36 #include "TypeOfCellTriangle.h"
37 #include <vtkm/CellShape.h>
38 #include <vtkm/CellTraits.h>
39 #include <vtkm/ErrorCode.h>
40 #include <vtkm/VecTraits.h>
41 #include <vtkm/VectorAnalysis.h>
42 
43 namespace vtkm
44 {
45 namespace worklet
46 {
47 namespace cellmetrics
48 {
49 
50 // By default, cells have no shape unless the shape type template is specialized below.
51 template <typename OutType, typename PointCoordVecType, typename CellShapeType>
52 VTKM_EXEC OutType CellShapeMetric(const vtkm::IdComponent& numPts,
53  const PointCoordVecType& pts,
54  CellShapeType shape,
55  vtkm::ErrorCode&)
56 {
57  UNUSED(numPts);
58  UNUSED(pts);
59  UNUSED(shape);
60  return OutType(-1.0);
61 }
62 
63 // =============================== Shape metric cells ==================================
64 
65 // Compute the shape quality metric of a triangular cell.
66 template <typename OutType, typename PointCoordVecType>
67 VTKM_EXEC OutType CellShapeMetric(const vtkm::IdComponent& numPts,
68  const PointCoordVecType& pts,
69  vtkm::CellShapeTagTriangle,
70  vtkm::ErrorCode& ec)
71 {
72  if (numPts != 3)
73  {
74  ec = vtkm::ErrorCode::InvalidNumberOfPoints;
75  return OutType(0.0);
76  }
77  using Scalar = OutType;
78  using CollectionOfPoints = PointCoordVecType;
79 
80  const Scalar condition =
81  vtkm::worklet::cellmetrics::CellConditionMetric<Scalar, CollectionOfPoints>(
82  numPts, pts, vtkm::CellShapeTagTriangle(), ec);
83  const Scalar q(1 / condition);
84 
85  return q;
86 }
87 
89 template <typename OutType, typename PointCoordVecType>
90 VTKM_EXEC OutType CellShapeMetric(const vtkm::IdComponent& numPts,
91  const PointCoordVecType& pts,
92  vtkm::CellShapeTagQuad,
93  vtkm::ErrorCode& ec)
94 {
95  if (numPts != 4)
96  {
97  ec = vtkm::ErrorCode::InvalidNumberOfPoints;
98  return OutType(0.0);
99  }
100 
101  using Scalar = OutType;
102  using CollectionOfPoints = PointCoordVecType;
103  using Vector = typename PointCoordVecType::ComponentType;
104  const Scalar two(2.0);
105  const Scalar alpha0 = GetQuadAlpha0<Scalar, Vector, CollectionOfPoints>(pts);
106  const Scalar alpha1 = GetQuadAlpha1<Scalar, Vector, CollectionOfPoints>(pts);
107  const Scalar alpha2 = GetQuadAlpha2<Scalar, Vector, CollectionOfPoints>(pts);
108  const Scalar alpha3 = GetQuadAlpha3<Scalar, Vector, CollectionOfPoints>(pts);
109  const Scalar l0Squared =
110  vtkm::Pow(GetQuadL0Magnitude<Scalar, Vector, CollectionOfPoints>(pts), 2);
111  const Scalar l1Squared =
112  vtkm::Pow(GetQuadL1Magnitude<Scalar, Vector, CollectionOfPoints>(pts), 2);
113  const Scalar l2Squared =
114  vtkm::Pow(GetQuadL2Magnitude<Scalar, Vector, CollectionOfPoints>(pts), 2);
115  const Scalar l3Squared =
116  vtkm::Pow(GetQuadL3Magnitude<Scalar, Vector, CollectionOfPoints>(pts), 2);
117 
118  const Scalar min = vtkm::Min(
119  (alpha0 / (l0Squared + l3Squared)),
120  vtkm::Min((alpha1 / (l1Squared + l0Squared)),
121  vtkm::Min((alpha2 / (l2Squared + l1Squared)), (alpha3 / (l3Squared + l2Squared)))));
122  const Scalar q(two * min);
123 
124  return q;
125 }
126 
127 // ============================= 3D cells ==================================
129 template <typename OutType, typename PointCoordVecType>
130 VTKM_EXEC OutType CellShapeMetric(const vtkm::IdComponent& numPts,
131  const PointCoordVecType& pts,
132  vtkm::CellShapeTagTetra,
133  vtkm::ErrorCode& ec)
134 {
135  if (numPts != 4)
136  {
137  ec = vtkm::ErrorCode::InvalidNumberOfPoints;
138  return OutType(0.0);
139  }
140  using Scalar = OutType;
141  using CollectionOfPoints = PointCoordVecType;
142  using Vector = typename PointCoordVecType::ComponentType;
143 
144  const Scalar zero(0.0);
145  const Scalar twoThirds = (Scalar)(2.0 / 3.0);
146  const Scalar threeHalves(1.5);
147  const Scalar rtTwo = (Scalar)(vtkm::Sqrt(2.0));
148  const Scalar three(3.0);
149  const Scalar jacobian =
150  vtkm::worklet::cellmetrics::CellJacobianMetric<Scalar, CollectionOfPoints>(
151  numPts, pts, vtkm::CellShapeTagTetra(), ec);
152 
153  if (jacobian <= zero)
154  {
155  return zero;
156  }
157 
158  const Vector l0 = GetTetraL0<Scalar, Vector, CollectionOfPoints>(pts);
159  const Vector l2 = GetTetraL2<Scalar, Vector, CollectionOfPoints>(pts);
160  const Vector l3 = GetTetraL3<Scalar, Vector, CollectionOfPoints>(pts);
161  const Vector negl2 = -1 * l2;
162 
163  const Scalar l0l0 = static_cast<Scalar>(vtkm::Dot(l0, l0));
164  const Scalar l2l2 = static_cast<Scalar>(vtkm::Dot(l2, l2));
165  const Scalar l3l3 = static_cast<Scalar>(vtkm::Dot(l3, l3));
166  const Scalar l0negl2 = static_cast<Scalar>(vtkm::Dot(l0, negl2));
167  const Scalar l0l3 = static_cast<Scalar>(vtkm::Dot(l0, l3));
168  const Scalar negl2l3 = static_cast<Scalar>(vtkm::Dot(negl2, l3));
169 
170  const Scalar numerator = three * vtkm::Pow(jacobian * rtTwo, twoThirds);
171  Scalar denominator = (threeHalves * (l0l0 + l2l2 + l3l3)) - (l0negl2 + l0l3 + negl2l3);
172  if (denominator <= zero)
173  {
174  return zero;
175  }
176  Scalar q(numerator / denominator);
177  return q;
178 }
179 
181 template <typename OutType, typename PointCoordVecType>
182 VTKM_EXEC OutType CellShapeMetric(const vtkm::IdComponent& numPts,
183  const PointCoordVecType& pts,
184  vtkm::CellShapeTagHexahedron,
185  vtkm::ErrorCode& ec)
186 {
187  if (numPts != 8)
188  {
189  ec = vtkm::ErrorCode::InvalidNumberOfPoints;
190  return OutType(0.0);
191  }
192  using Scalar = OutType;
193  using CollectionOfPoints = PointCoordVecType;
194  using Vector = typename PointCoordVecType::ComponentType;
195 
196  const Scalar zero(0.0);
197  const Scalar twoThirds = (Scalar)(2.0 / 3.0);
198 
199  const Scalar alpha0 = GetHexAlphai<Scalar, Vector, CollectionOfPoints>(pts, vtkm::Id(0));
200  const Scalar alpha1 = GetHexAlphai<Scalar, Vector, CollectionOfPoints>(pts, vtkm::Id(1));
201  const Scalar alpha2 = GetHexAlphai<Scalar, Vector, CollectionOfPoints>(pts, vtkm::Id(2));
202  const Scalar alpha3 = GetHexAlphai<Scalar, Vector, CollectionOfPoints>(pts, vtkm::Id(3));
203  const Scalar alpha4 = GetHexAlphai<Scalar, Vector, CollectionOfPoints>(pts, vtkm::Id(4));
204  const Scalar alpha5 = GetHexAlphai<Scalar, Vector, CollectionOfPoints>(pts, vtkm::Id(5));
205  const Scalar alpha6 = GetHexAlphai<Scalar, Vector, CollectionOfPoints>(pts, vtkm::Id(6));
206  const Scalar alpha7 = GetHexAlphai<Scalar, Vector, CollectionOfPoints>(pts, vtkm::Id(7));
207  const Scalar alpha8 = GetHexAlphai<Scalar, Vector, CollectionOfPoints>(pts, vtkm::Id(8));
208 
209  const Scalar A0squared =
210  GetHexAiNormSquared<Scalar, Vector, CollectionOfPoints>(pts, vtkm::Id(0));
211  const Scalar A1squared =
212  GetHexAiNormSquared<Scalar, Vector, CollectionOfPoints>(pts, vtkm::Id(1));
213  const Scalar A2squared =
214  GetHexAiNormSquared<Scalar, Vector, CollectionOfPoints>(pts, vtkm::Id(2));
215  const Scalar A3squared =
216  GetHexAiNormSquared<Scalar, Vector, CollectionOfPoints>(pts, vtkm::Id(3));
217  const Scalar A4squared =
218  GetHexAiNormSquared<Scalar, Vector, CollectionOfPoints>(pts, vtkm::Id(4));
219  const Scalar A5squared =
220  GetHexAiNormSquared<Scalar, Vector, CollectionOfPoints>(pts, vtkm::Id(5));
221  const Scalar A6squared =
222  GetHexAiNormSquared<Scalar, Vector, CollectionOfPoints>(pts, vtkm::Id(6));
223  const Scalar A7squared =
224  GetHexAiNormSquared<Scalar, Vector, CollectionOfPoints>(pts, vtkm::Id(7));
225  const Scalar A8squared =
226  GetHexAiNormSquared<Scalar, Vector, CollectionOfPoints>(pts, vtkm::Id(8));
227 
228  if (alpha0 <= zero || alpha1 <= zero || alpha2 <= zero || alpha3 <= zero || alpha4 <= zero ||
229  alpha5 <= zero || alpha6 <= zero || alpha7 <= zero || alpha8 <= zero || A0squared <= zero ||
230  A1squared <= zero || A2squared <= zero || A3squared <= zero || A4squared <= zero ||
231  A5squared <= zero || A6squared <= zero || A7squared <= zero || A8squared <= zero)
232  {
233  return zero;
234  }
235  // find min according to verdict manual
236  Scalar min = vtkm::Pow(alpha0, twoThirds) / A0squared;
237  Scalar temp = vtkm::Pow(alpha1, twoThirds) / A1squared;
238  min = temp < min ? temp : min;
239 
240  temp = vtkm::Pow(alpha2, twoThirds) / A2squared;
241  min = temp < min ? temp : min;
242 
243  temp = vtkm::Pow(alpha3, twoThirds) / A3squared;
244  min = temp < min ? temp : min;
245 
246  temp = vtkm::Pow(alpha4, twoThirds) / A4squared;
247  min = temp < min ? temp : min;
248 
249  temp = vtkm::Pow(alpha5, twoThirds) / A5squared;
250  min = temp < min ? temp : min;
251 
252  temp = vtkm::Pow(alpha6, twoThirds) / A6squared;
253  min = temp < min ? temp : min;
254 
255  temp = vtkm::Pow(alpha7, twoThirds) / A7squared;
256  min = temp < min ? temp : min;
257 
258  temp = vtkm::Pow(alpha8, twoThirds) / A8squared;
259  min = temp < min ? temp : min;
260  // determine the shape
261  Scalar q(3 * min);
262  return q;
263 }
264 } // cellmetrics
265 } // worklet
266 } // vtkm
267 #endif // vtk_m_worklet_CellShapeMetric_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_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
CellShape.h
vtkm::Id
vtkm::Int32 Id
Represents an ID (index into arrays).
Definition: Types.h:191
ErrorCode.h
VectorAnalysis.h
vtkm::worklet::cellmetrics::CellShapeMetric
VTKM_EXEC OutType CellShapeMetric(const vtkm::IdComponent &numPts, const PointCoordVecType &pts, CellShapeType shape, vtkm::ErrorCode &)
Definition: CellShapeMetric.h:52
CellJacobianMetric.h
TypeOfCellTetrahedral.h
TypeOfCellTriangle.h
TypeOfCellHexahedral.h
TypeOfCellQuadrilateral.h
CellTraits.h
CellConditionMetric.h
VecTraits.h
vtkm::ErrorCode::InvalidNumberOfPoints
@ InvalidNumberOfPoints
Generated by   doxygen 1.8.17