VTK-m  2.0
CellDimensionMetric.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 2018 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
10 // Copyright 2018 UT-Battelle, LLC.
11 // Copyright 2018 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_CellDimensionMetric_h
21 #define vtk_m_worklet_cellmetrics_CellDimensionMetric_h
22 
23 /*
24  * Mesh quality metric functions that compute the dimension of mesh cells. The
25  * Dimension metric is defined as the volume divided by two times the gradient
26  * of the volume.
27  *
28  * This metric was designed in context of Sandia's Pronto code for stable time
29  * step calculation.
30  *
31  * These metric computations are adapted from the VTK implementation of the
32  * Verdict library, which provides a set of mesh/cell metrics for evaluating the
33  * geometric qualities of regions of mesh spaces.
34  *
35  * See: The Verdict Library Reference Manual (for per-cell-type metric formulae)
36  * See: vtk/ThirdParty/verdict/vtkverdict (for VTK code implementation of this
37  * metric)
38  */
39 
40 #include <vtkm/CellShape.h>
41 #include <vtkm/CellTraits.h>
42 #include <vtkm/ErrorCode.h>
43 #include <vtkm/VecTraits.h>
44 #include <vtkm/VectorAnalysis.h>
45 
46 #define UNUSED(expr) (void)(expr);
47 
48 namespace vtkm
49 {
50 namespace worklet
51 {
52 namespace cellmetrics
53 {
54 
55 using FloatType = vtkm::FloatDefault;
56 
57 // ========================= Unsupported cells ==================================
58 
59 // Dimension is only defined for Hexahedral cell typesk
60 template <typename OutType, typename PointCoordVecType, typename CellShapeType>
61 VTKM_EXEC OutType CellDimensionMetric(const vtkm::IdComponent& numPts,
62  const PointCoordVecType& pts,
63  CellShapeType shape,
64  vtkm::ErrorCode&)
65 {
66  UNUSED(numPts);
67  UNUSED(pts);
68  UNUSED(shape);
69  return OutType(-1.);
70 }
71 
72 template <typename OutType, typename PointCoordVecType>
73 VTKM_EXEC OutType CellDimensionMetric(const vtkm::IdComponent& numPts,
74  const PointCoordVecType& pts,
75  vtkm::CellShapeTagHexahedron,
76  vtkm::ErrorCode& ec)
77 {
78  UNUSED(numPts);
79  UNUSED(ec);
80 
81  OutType gradop[8][3];
82  OutType x1 = static_cast<OutType>(pts[0][0]);
83  OutType x2 = static_cast<OutType>(pts[1][0]);
84  OutType x3 = static_cast<OutType>(pts[2][0]);
85  OutType x4 = static_cast<OutType>(pts[3][0]);
86  OutType x5 = static_cast<OutType>(pts[4][0]);
87  OutType x6 = static_cast<OutType>(pts[5][0]);
88  OutType x7 = static_cast<OutType>(pts[6][0]);
89  OutType x8 = static_cast<OutType>(pts[7][0]);
90 
91  OutType y1 = static_cast<OutType>(pts[0][1]);
92  OutType y2 = static_cast<OutType>(pts[1][1]);
93  OutType y3 = static_cast<OutType>(pts[2][1]);
94  OutType y4 = static_cast<OutType>(pts[3][1]);
95  OutType y5 = static_cast<OutType>(pts[4][1]);
96  OutType y6 = static_cast<OutType>(pts[5][1]);
97  OutType y7 = static_cast<OutType>(pts[6][1]);
98  OutType y8 = static_cast<OutType>(pts[7][1]);
99 
100  OutType z1 = static_cast<OutType>(pts[0][2]);
101  OutType z2 = static_cast<OutType>(pts[1][2]);
102  OutType z3 = static_cast<OutType>(pts[2][2]);
103  OutType z4 = static_cast<OutType>(pts[3][2]);
104  OutType z5 = static_cast<OutType>(pts[4][2]);
105  OutType z6 = static_cast<OutType>(pts[5][2]);
106  OutType z7 = static_cast<OutType>(pts[6][2]);
107  OutType z8 = static_cast<OutType>(pts[7][2]);
108 
109  OutType z24 = z2 - z4;
110  OutType z52 = z5 - z2;
111  OutType z45 = z4 - z5;
112  gradop[0][0] = (y2 * (z6 - z3 - z45) + y3 * z24 + y4 * (z3 - z8 - z52) + y5 * (z8 - z6 - z24) +
113  y6 * z52 + y8 * z45) /
114  (OutType)12.0;
115 
116  OutType z31 = z3 - z1;
117  OutType z63 = z6 - z3;
118  OutType z16 = z1 - z6;
119  gradop[1][0] = (y3 * (z7 - z4 - z16) + y4 * z31 + y1 * (z4 - z5 - z63) + y6 * (z5 - z7 - z31) +
120  y7 * z63 + y5 * z16) /
121  (OutType)12.0;
122 
123  OutType z42 = z4 - z2;
124  OutType z74 = z7 - z4;
125  OutType z27 = z2 - z7;
126  gradop[2][0] = (y4 * (z8 - z1 - z27) + y1 * z42 + y2 * (z1 - z6 - z74) + y7 * (z6 - z8 - z42) +
127  y8 * z74 + y6 * z27) /
128  (OutType)12.0;
129 
130  OutType z13 = z1 - z3;
131  OutType z81 = z8 - z1;
132  OutType z38 = z3 - z8;
133  gradop[3][0] = (y1 * (z5 - z2 - z38) + y2 * z13 + y3 * (z2 - z7 - z81) + y8 * (z7 - z5 - z13) +
134  y5 * z81 + y7 * z38) /
135  (OutType)12.0;
136 
137  OutType z86 = z8 - z6;
138  OutType z18 = z1 - z8;
139  OutType z61 = z6 - z1;
140  gradop[4][0] = (y8 * (z4 - z7 - z61) + y7 * z86 + y6 * (z7 - z2 - z18) + y1 * (z2 - z4 - z86) +
141  y4 * z18 + y2 * z61) /
142  (OutType)12.0;
143 
144  OutType z57 = z5 - z7;
145  OutType z25 = z2 - z5;
146  OutType z72 = z7 - z2;
147  gradop[5][0] = (y5 * (z1 - z8 - z72) + y8 * z57 + y7 * (z8 - z3 - z25) + y2 * (z3 - z1 - z57) +
148  y1 * z25 + y3 * z72) /
149  (OutType)12.0;
150 
151  OutType z68 = z6 - z8;
152  OutType z36 = z3 - z6;
153  OutType z83 = z8 - z3;
154  gradop[6][0] = (y6 * (z2 - z5 - z83) + y5 * z68 + y8 * (z5 - z4 - z36) + y3 * (z4 - z2 - z68) +
155  y2 * z36 + y4 * z83) /
156  (OutType)12.0;
157 
158  OutType z75 = z7 - z5;
159  OutType z47 = z4 - z7;
160  OutType z54 = z5 - z4;
161  gradop[7][0] = (y7 * (z3 - z6 - z54) + y6 * z75 + y5 * (z6 - z1 - z47) + y4 * (z1 - z3 - z75) +
162  y3 * z47 + y1 * z54) /
163  (OutType)12.0;
164 
165  OutType x24 = x2 - x4;
166  OutType x52 = x5 - x2;
167  OutType x45 = x4 - x5;
168  gradop[0][1] = (z2 * (x6 - x3 - x45) + z3 * x24 + z4 * (x3 - x8 - x52) + z5 * (x8 - x6 - x24) +
169  z6 * x52 + z8 * x45) /
170  (OutType)12.0;
171 
172  OutType x31 = x3 - x1;
173  OutType x63 = x6 - x3;
174  OutType x16 = x1 - x6;
175  gradop[1][1] = (z3 * (x7 - x4 - x16) + z4 * x31 + z1 * (x4 - x5 - x63) + z6 * (x5 - x7 - x31) +
176  z7 * x63 + z5 * x16) /
177  (OutType)12.0;
178 
179  OutType x42 = x4 - x2;
180  OutType x74 = x7 - x4;
181  OutType x27 = x2 - x7;
182  gradop[2][1] = (z4 * (x8 - x1 - x27) + z1 * x42 + z2 * (x1 - x6 - x74) + z7 * (x6 - x8 - x42) +
183  z8 * x74 + z6 * x27) /
184  (OutType)12.0;
185 
186  OutType x13 = x1 - x3;
187  OutType x81 = x8 - x1;
188  OutType x38 = x3 - x8;
189  gradop[3][1] = (z1 * (x5 - x2 - x38) + z2 * x13 + z3 * (x2 - x7 - x81) + z8 * (x7 - x5 - x13) +
190  z5 * x81 + z7 * x38) /
191  (OutType)12.0;
192 
193  OutType x86 = x8 - x6;
194  OutType x18 = x1 - x8;
195  OutType x61 = x6 - x1;
196  gradop[4][1] = (z8 * (x4 - x7 - x61) + z7 * x86 + z6 * (x7 - x2 - x18) + z1 * (x2 - x4 - x86) +
197  z4 * x18 + z2 * x61) /
198  (OutType)12.0;
199 
200  OutType x57 = x5 - x7;
201  OutType x25 = x2 - x5;
202  OutType x72 = x7 - x2;
203  gradop[5][1] = (z5 * (x1 - x8 - x72) + z8 * x57 + z7 * (x8 - x3 - x25) + z2 * (x3 - x1 - x57) +
204  z1 * x25 + z3 * x72) /
205  (OutType)12.0;
206 
207  OutType x68 = x6 - x8;
208  OutType x36 = x3 - x6;
209  OutType x83 = x8 - x3;
210  gradop[6][1] = (z6 * (x2 - x5 - x83) + z5 * x68 + z8 * (x5 - x4 - x36) + z3 * (x4 - x2 - x68) +
211  z2 * x36 + z4 * x83) /
212  (OutType)12.0;
213 
214  OutType x75 = x7 - x5;
215  OutType x47 = x4 - x7;
216  OutType x54 = x5 - x4;
217  gradop[7][1] = (z7 * (x3 - x6 - x54) + z6 * x75 + z5 * (x6 - x1 - x47) + z4 * (x1 - x3 - x75) +
218  z3 * x47 + z1 * x54) /
219  (OutType)12.0;
220 
221  OutType y24 = y2 - y4;
222  OutType y52 = y5 - y2;
223  OutType y45 = y4 - y5;
224  gradop[0][2] = (x2 * (y6 - y3 - y45) + x3 * y24 + x4 * (y3 - y8 - y52) + x5 * (y8 - y6 - y24) +
225  x6 * y52 + x8 * y45) /
226  (OutType)12.0;
227 
228  OutType y31 = y3 - y1;
229  OutType y63 = y6 - y3;
230  OutType y16 = y1 - y6;
231  gradop[1][2] = (x3 * (y7 - y4 - y16) + x4 * y31 + x1 * (y4 - y5 - y63) + x6 * (y5 - y7 - y31) +
232  x7 * y63 + x5 * y16) /
233  (OutType)12.0;
234 
235  OutType y42 = y4 - y2;
236  OutType y74 = y7 - y4;
237  OutType y27 = y2 - y7;
238  gradop[2][2] = (x4 * (y8 - y1 - y27) + x1 * y42 + x2 * (y1 - y6 - y74) + x7 * (y6 - y8 - y42) +
239  x8 * y74 + x6 * y27) /
240  (OutType)12.0;
241 
242  OutType y13 = y1 - y3;
243  OutType y81 = y8 - y1;
244  OutType y38 = y3 - y8;
245  gradop[3][2] = (x1 * (y5 - y2 - y38) + x2 * y13 + x3 * (y2 - y7 - y81) + x8 * (y7 - y5 - y13) +
246  x5 * y81 + x7 * y38) /
247  (OutType)12.0;
248 
249  OutType y86 = y8 - y6;
250  OutType y18 = y1 - y8;
251  OutType y61 = y6 - y1;
252  gradop[4][2] = (x8 * (y4 - y7 - y61) + x7 * y86 + x6 * (y7 - y2 - y18) + x1 * (y2 - y4 - y86) +
253  x4 * y18 + x2 * y61) /
254  (OutType)12.0;
255 
256  OutType y57 = y5 - y7;
257  OutType y25 = y2 - y5;
258  OutType y72 = y7 - y2;
259  gradop[5][2] = (x5 * (y1 - y8 - y72) + x8 * y57 + x7 * (y8 - y3 - y25) + x2 * (y3 - y1 - y57) +
260  x1 * y25 + x3 * y72) /
261  (OutType)12.0;
262 
263  OutType y68 = y6 - y8;
264  OutType y36 = y3 - y6;
265  OutType y83 = y8 - y3;
266  gradop[6][2] = (x6 * (y2 - y5 - y83) + x5 * y68 + x8 * (y5 - y4 - y36) + x3 * (y4 - y2 - y68) +
267  x2 * y36 + x4 * y83) /
268  (OutType)12.0;
269 
270  OutType y75 = y7 - y5;
271  OutType y47 = y4 - y7;
272  OutType y54 = y5 - y4;
273  gradop[7][2] = (x7 * (y3 - y6 - y54) + x6 * y75 + x5 * (y6 - y1 - y47) + x4 * (y1 - y3 - y75) +
274  x3 * y47 + x1 * y54) /
275  (OutType)12.0;
276  OutType volume = OutType(pts[0][0]) * gradop[0][0] + OutType(pts[1][0]) * gradop[1][0] +
277  OutType(pts[2][0]) * gradop[2][0] + OutType(pts[3][0]) * gradop[3][0] +
278  OutType(pts[4][0]) * gradop[4][0] + OutType(pts[5][0]) * gradop[5][0] +
279  OutType(pts[6][0]) * gradop[6][0] + OutType(pts[7][0]) * gradop[7][0];
280  OutType two = (OutType)2.;
281  OutType aspect = (OutType).5 * vtkm::Pow(volume, two) /
282  (vtkm::Pow(gradop[0][0], two) + vtkm::Pow(gradop[1][0], two) + vtkm::Pow(gradop[2][0], two) +
283  vtkm::Pow(gradop[3][0], two) + vtkm::Pow(gradop[4][0], two) + vtkm::Pow(gradop[5][0], two) +
284  vtkm::Pow(gradop[6][0], two) + vtkm::Pow(gradop[7][0], two) + vtkm::Pow(gradop[0][1], two) +
285  vtkm::Pow(gradop[1][1], two) + vtkm::Pow(gradop[2][1], two) + vtkm::Pow(gradop[3][1], two) +
286  vtkm::Pow(gradop[4][1], two) + vtkm::Pow(gradop[5][1], two) + vtkm::Pow(gradop[6][1], two) +
287  vtkm::Pow(gradop[7][1], two) + vtkm::Pow(gradop[0][2], two) + vtkm::Pow(gradop[1][2], two) +
288  vtkm::Pow(gradop[2][2], two) + vtkm::Pow(gradop[3][2], two) + vtkm::Pow(gradop[4][2], two) +
289  vtkm::Pow(gradop[5][2], two) + vtkm::Pow(gradop[6][2], two) + vtkm::Pow(gradop[7][2], two));
290 
291  return vtkm::Sqrt(aspect);
292 }
293 
294 } // cell metrics
295 } // worklet
296 } // vtkm
297 
298 #endif // vtk_m_worklet_cellmetrics_CellDimensionMetric_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: CellDimensionMetric.h:46
vtkm::IdComponent
vtkm::Int32 IdComponent
Represents a component ID (index of component in a vector).
Definition: Types.h:168
vtkm::worklet::cellmetrics::FloatType
vtkm::FloatDefault FloatType
Definition: CellAspectFrobeniusMetric.h:50
CellShape.h
ErrorCode.h
VectorAnalysis.h
vtkm::FloatDefault
vtkm::Float32 FloatDefault
The floating point type to use when no other precision is specified.
Definition: Types.h:198
vtkm::worklet::cellmetrics::CellDimensionMetric
VTKM_EXEC OutType CellDimensionMetric(const vtkm::IdComponent &numPts, const PointCoordVecType &pts, CellShapeType shape, vtkm::ErrorCode &)
Definition: CellDimensionMetric.h:61
CellTraits.h
VecTraits.h
Generated by   doxygen 1.8.17