VTK-m  2.0
CellScaledJacobianMetric.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_ScaledJacobian_h
21 #define vtk_m_worklet_cellmetrics_ScaledJacobian_h
22 
23 /*
24 * Mesh quality metric functions that computes the scaled jacobian of mesh cells.
25 * The jacobian of a cell is defined as the determinant of the Jociabian matrix
26 ** These metric computations are adapted from the VTK implementation of the Verdict library,
27 * which provides a set of mesh/cell metrics for evaluating the geometric qualities of regions
28 * of mesh spaces.
29 ** See: The Verdict Library Reference Manual (for per-cell-type metric formulae)
30 * See: vtk/ThirdParty/verdict/vtkverdict (for VTK code implementation of this metric)
31 */
32 
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 #define UNUSED(expr) (void)(expr);
44 
45 namespace vtkm
46 {
47 namespace worklet
48 {
49 namespace cellmetrics
50 {
51 
52 using FloatType = vtkm::FloatDefault;
53 
54 // ========================= Unsupported cells ==================================
55 
56 // By default, cells have zero shape unless the shape type template is specialized below.
57 template <typename OutType, typename PointCoordVecType, typename CellShapeType>
58 VTKM_EXEC OutType CellScaledJacobianMetric(const vtkm::IdComponent& numPts,
59  const PointCoordVecType& pts,
60  CellShapeType shape,
61  vtkm::ErrorCode&)
62 {
63  UNUSED(numPts);
64  UNUSED(pts);
65  UNUSED(shape);
66  return OutType(-2.0);
67 }
68 
69 // ========================= 2D cells ==================================
70 //Compute the scaled jacobian of a triangle
71 //Formula: q = ((2*sqrt(3))/3) * (J/Lmax)
72 // - J -> jacobian, if surface normal N is center of triangle and
73 // and N*L2*L1 < 0, then -jacobian
74 // - Lmax -> max{ |L0| * |L1|, |L0| * |L2|, |L1| * |L2| }
75 //Equals 1 for equilateral unit triangle
76 //Acceptable Range: [0.5, 2*sqrt(3)/3]
77 //Normal Range : [ -(2*sqrt(3)/3) , 2*sqrt(3)/3]
78 //Full Range : [-FLOAT_MAX, FLOAT_MAX]
79 template <typename OutType, typename PointCoordVecType>
80 VTKM_EXEC OutType CellScaledJacobianMetric(const vtkm::IdComponent& numPts,
81  const PointCoordVecType& pts,
82  vtkm::CellShapeTagTriangle,
83  vtkm::ErrorCode& ec)
84 {
85  if (numPts != 3)
86  {
87  ec = vtkm::ErrorCode::InvalidNumberOfPoints;
88  return OutType(0.0);
89  }
90  using Scalar = OutType;
91  using CollectionOfPoints = PointCoordVecType;
92  using Vector = typename PointCoordVecType::ComponentType;
93 
94  const Vector l0 = GetTriangleL0<Scalar, Vector, CollectionOfPoints>(pts);
95  const Vector l1 = GetTriangleL1<Scalar, Vector, CollectionOfPoints>(pts);
96  const Vector l2 = GetTriangleL2<Scalar, Vector, CollectionOfPoints>(pts);
97 
98  const Scalar modifier = (Scalar)(2 * vtkm::Sqrt(3) / 3);
99  const Scalar l0_magnitude = GetTriangleL0Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
100  const Scalar l1_magnitude = GetTriangleL1Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
101  const Scalar l2_magnitude = GetTriangleL2Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
102  const Scalar l0l1_product = l0_magnitude * l1_magnitude;
103  const Scalar l0l2_product = l0_magnitude * l2_magnitude;
104  const Scalar l1l2_product = l1_magnitude * l2_magnitude;
105  const Scalar productMax = vtkm::Max(l0l1_product, vtkm::Max(l0l2_product, l1l2_product));
106  if (productMax <= Scalar(0.0))
107  {
108  return Scalar(0.0);
109  }
110  // compute jacobian of triangle
111  Vector TriCross = vtkm::Cross(l2, l1);
112  Scalar scaledJacobian = static_cast<OutType>(vtkm::Magnitude(TriCross));
113 
114  //add all pieces together
115  //TODO change
116  Vector normalVector = (Scalar(1.0) / Scalar(3.0)) * (l0 + l1 + l2);
117  Vector surfaceNormalAtCenter = vtkm::TriangleNormal(normalVector, normalVector, normalVector);
118  if (vtkm::Dot(surfaceNormalAtCenter, TriCross) < 0)
119  {
120  scaledJacobian *= -1;
121  }
122  scaledJacobian *= modifier;
123  const Scalar q = scaledJacobian / productMax;
124 
125  return q;
126 }
127 
128 // Compute the scaled jacobian of a quadrilateral.
129 // Formula: min{J0/(L0*L3), J1/(L1*L0), J2/(L2*L1), J3/(L3*L2)}
130 // -Ji -> Jacobian at corner i, the intersection of the edge vectors
131 // it is divided by
132 // Equals 1 for a unit square
133 //Acceptable Range: [0.3, 1]
134 //Normal Range : [-1, 1]
135 // Full range : [-1, 1]
136 template <typename OutType, typename PointCoordVecType>
137 VTKM_EXEC OutType CellScaledJacobianMetric(const vtkm::IdComponent& numPts,
138  const PointCoordVecType& pts,
139  vtkm::CellShapeTagQuad,
140  vtkm::ErrorCode& ec)
141 {
142  if (numPts != 4)
143  {
144  ec = vtkm::ErrorCode::InvalidNumberOfPoints;
145  return OutType(0.0);
146  }
147  using Scalar = OutType;
148  using CollectionOfPoints = PointCoordVecType;
149  using Vector = typename PointCoordVecType::ComponentType;
150  //The 4 edges of a quadrilateral
151  const Scalar l0_magnitude = GetQuadL0Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
152  const Scalar l1_magnitude = GetQuadL1Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
153  const Scalar l2_magnitude = GetQuadL2Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
154  const Scalar l3_magnitude = GetQuadL3Magnitude<Scalar, Vector, CollectionOfPoints>(pts);
155  const Scalar negativeInfinity = vtkm::NegativeInfinity<Scalar>();
156 
157  if (l0_magnitude < negativeInfinity || l1_magnitude < negativeInfinity ||
158  l2_magnitude < negativeInfinity || l3_magnitude < negativeInfinity)
159  {
160  return Scalar(0.0);
161  }
162  const Scalar l0l3_product = l0_magnitude * l3_magnitude;
163  const Scalar l1l0_product = l1_magnitude * l0_magnitude;
164  const Scalar l2l1_product = l2_magnitude * l1_magnitude;
165  const Scalar l3l2_product = l3_magnitude * l2_magnitude;
166 
167  /*
168 3 * 0
169 0 * 1
170 1 * 2
171 2 * 3
172 */
173 
174  Scalar alpha0Scaled = GetQuadAlpha0<Scalar, Vector, CollectionOfPoints>(pts);
175  Scalar alpha1Scaled = GetQuadAlpha1<Scalar, Vector, CollectionOfPoints>(pts);
176  Scalar alpha2Scaled = GetQuadAlpha2<Scalar, Vector, CollectionOfPoints>(pts);
177  Scalar alpha3Scaled = GetQuadAlpha3<Scalar, Vector, CollectionOfPoints>(pts);
178 
179  alpha0Scaled /= l0l3_product;
180  alpha1Scaled /= l1l0_product;
181  alpha2Scaled /= l2l1_product;
182  alpha3Scaled /= l3l2_product;
183  const Scalar q =
184  vtkm::Min(alpha0Scaled, vtkm::Min(alpha1Scaled, vtkm::Min(alpha2Scaled, alpha3Scaled)));
185  return q;
186 }
187 
188 // ============================= 3D Volume cells ==================================
189 // Compute the scaled jacobian of a hexahedron.
190 // Formula: q = min{Ai}
191 // Ai -> for i 1...8 (Jacobian determinant at respective corner, divided by corresponding edge lengths
192 // Equals 1 for a unit cube
193 // Acceptable Range: [0.5, 1]
194 // Normal Range : [-1, 1]
195 // Full range : [1,FLOAT_MAX]
196 template <typename OutType, typename PointCoordVecType>
197 VTKM_EXEC OutType CellScaledJacobianMetric(const vtkm::IdComponent& numPts,
198  const PointCoordVecType& pts,
199  vtkm::CellShapeTagHexahedron,
200  vtkm::ErrorCode& ec)
201 {
202  if (numPts != 8)
203  {
204  ec = vtkm::ErrorCode::InvalidNumberOfPoints;
205  return OutType(0.0);
206  }
207 
208  //The 12 edges of a hexahedron
209  using Edge = typename PointCoordVecType::ComponentType;
210  Edge HexEdges[12] = { pts[1] - pts[0], pts[2] - pts[1], pts[3] - pts[2], pts[3] - pts[0],
211  pts[4] - pts[0], pts[5] - pts[1], pts[6] - pts[2], pts[7] - pts[3],
212  pts[5] - pts[4], pts[6] - pts[5], pts[7] - pts[6], pts[7] - pts[4] };
213  Edge principleXAxis = HexEdges[0] + (pts[2] - pts[3]) + HexEdges[8] + (pts[6] - pts[7]);
214  Edge principleYAxis = (pts[3] - pts[0]) + HexEdges[1] + (pts[7] - pts[4]) + HexEdges[9];
215  Edge principleZAxis = HexEdges[4] + HexEdges[5] + HexEdges[6] + HexEdges[7];
216  Edge hexMatrices[9][3] = { { HexEdges[0], HexEdges[3], HexEdges[4] },
217  { HexEdges[1], (-1 * HexEdges[0]), HexEdges[5] },
218  { HexEdges[2], (-1 * HexEdges[1]), HexEdges[6] },
219  { (-1 * HexEdges[3]), (-1 * HexEdges[2]), HexEdges[7] },
220  { HexEdges[11], HexEdges[8], (-1 * HexEdges[4]) },
221  { (-1 * HexEdges[8]), HexEdges[9], (-1 * HexEdges[5]) },
222  { (-1 * HexEdges[9]), HexEdges[10], (-1 * HexEdges[6]) },
223  { (-1 * HexEdges[10]), (-1 * HexEdges[11]), (-1 * HexEdges[7]) },
224  { principleXAxis, principleYAxis, principleZAxis } };
225  OutType currDeterminant, minDeterminant = vtkm::Infinity<OutType>();
226  FloatType lenSquared1, lenSquared2, lenSquared3, minLengthSquared = vtkm::Infinity<FloatType>();
227  vtkm::IdComponent matrixIndex;
228  for (matrixIndex = 0; matrixIndex < 9; matrixIndex++)
229  {
230  lenSquared1 = (FloatType)vtkm::MagnitudeSquared(hexMatrices[matrixIndex][0]);
231  minLengthSquared = lenSquared1 < minLengthSquared ? lenSquared1 : minLengthSquared;
232  lenSquared2 = (FloatType)vtkm::MagnitudeSquared(hexMatrices[matrixIndex][1]);
233  minLengthSquared = lenSquared2 < minLengthSquared ? lenSquared2 : minLengthSquared;
234  lenSquared3 = (FloatType)vtkm::MagnitudeSquared(hexMatrices[matrixIndex][2]);
235  minLengthSquared = lenSquared3 < minLengthSquared ? lenSquared3 : minLengthSquared;
236 
237  vtkm::Normalize(hexMatrices[matrixIndex][0]);
238  vtkm::Normalize(hexMatrices[matrixIndex][1]);
239  vtkm::Normalize(hexMatrices[matrixIndex][2]);
240  currDeterminant =
241  (OutType)vtkm::Dot(hexMatrices[matrixIndex][0],
242  vtkm::Cross(hexMatrices[matrixIndex][1], hexMatrices[matrixIndex][2]));
243  if (currDeterminant < minDeterminant)
244  {
245  minDeterminant = currDeterminant;
246  }
247  }
248  if (minLengthSquared < vtkm::NegativeInfinity<FloatType>())
249  {
250  return vtkm::Infinity<OutType>();
251  }
252  OutType toReturn = minDeterminant;
253  if (toReturn > 0)
254  return vtkm::Min(toReturn, vtkm::Infinity<OutType>()); //normal case
255 
256  return vtkm::Max(toReturn, vtkm::NegativeInfinity<OutType>());
257 }
258 
259 // Compute the scaled jacobian of a tetrahedron.
260 // Formula: q = J*sqrt(2)/Lamda_max
261 // J -> jacobian,(L2 * L0) * L3
262 // Lamda_max -> max{ L0*L2*L3, L0*L1*L4, L1*L2*L5, L3*L4*L5}
263 // Equals Sqrt(2) / 2 for unit equilateral tetrahedron
264 // Acceptable Range: [0, FLOAT_MAX]
265 // Normal Range: [0, FLOAT_MAX]
266 // Full range: [FLOAT_MIN,FLOAT_MAX]
267 template <typename OutType, typename PointCoordVecType>
268 VTKM_EXEC OutType CellScaledJacobianMetric(const vtkm::IdComponent& numPts,
269  const PointCoordVecType& pts,
270  vtkm::CellShapeTagTetra,
271  vtkm::ErrorCode& ec)
272 {
273  if (numPts != 4)
274  {
275  ec = vtkm::ErrorCode::InvalidNumberOfPoints;
276  return OutType(0.0);
277  }
278 
279  //the edge and side sets
280  using Edge = typename PointCoordVecType::ComponentType;
281  Edge Edges[6] = { pts[1] - pts[0], pts[2] - pts[1], pts[0] - pts[2],
282  pts[3] - pts[0], pts[3] - pts[1], pts[3] - pts[2] };
283  OutType EdgesSquared[6];
284  OutType jacobian = static_cast<OutType>(vtkm::Dot(vtkm::Cross(Edges[2], Edges[0]), Edges[3]));
285  // compute the scaled jacobian
286  OutType currSide, maxSide = vtkm::NegativeInfinity<OutType>();
287  vtkm::IdComponent edgeIndex, sideIndex;
288  for (edgeIndex = 0; edgeIndex < 6; edgeIndex++)
289  {
290  EdgesSquared[edgeIndex] = static_cast<OutType>(vtkm::MagnitudeSquared(Edges[edgeIndex]));
291  }
292  OutType Sides[4] = { EdgesSquared[0] * EdgesSquared[2] * EdgesSquared[3],
293  EdgesSquared[0] * EdgesSquared[1] * EdgesSquared[4],
294  EdgesSquared[1] * EdgesSquared[2] * EdgesSquared[5],
295  EdgesSquared[3] * EdgesSquared[4] * EdgesSquared[5] };
296  for (sideIndex = 0; sideIndex < 4; sideIndex++)
297  {
298  currSide = Sides[sideIndex];
299  maxSide = currSide > maxSide ? currSide : maxSide;
300  }
301  maxSide = vtkm::Sqrt(maxSide);
302  OutType toUseInCalculation = jacobian > maxSide ? jacobian : maxSide;
303  if (toUseInCalculation < vtkm::NegativeInfinity<OutType>())
304  {
305  return vtkm::Infinity<OutType>();
306  }
307  return (vtkm::Sqrt<OutType>(2) * jacobian) / toUseInCalculation;
308 }
309 } // namespace cellmetrics
310 } // namespace worklet
311 } // namespace vtkm
312 #endif // vtk_m_worklet_cellmetrics_ScaledJacobian_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
vtkm::worklet::cellmetrics::CellScaledJacobianMetric
VTKM_EXEC OutType CellScaledJacobianMetric(const vtkm::IdComponent &numPts, const PointCoordVecType &pts, CellShapeType shape, vtkm::ErrorCode &)
Definition: CellScaledJacobianMetric.h:58
vtkm::Normalize
VTKM_EXEC_CONT void Normalize(T &x)
Changes a vector to be normal.
Definition: VectorAnalysis.h:168
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
ErrorCode.h
VectorAnalysis.h
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
TypeOfCellTetrahedral.h
TypeOfCellTriangle.h
vtkm::TriangleNormal
VTKM_EXEC_CONT vtkm::Vec< typename detail::FloatingPointReturnType< T >::Type, 3 > TriangleNormal(const vtkm::Vec< T, 3 > &a, const vtkm::Vec< T, 3 > &b, const vtkm::Vec< T, 3 > &c)
Find the normal of a triangle.
Definition: VectorAnalysis.h:200
TypeOfCellHexahedral.h
UNUSED
#define UNUSED(expr)
Definition: CellScaledJacobianMetric.h:43
vtkm::FloatDefault
vtkm::Float32 FloatDefault
The floating point type to use when no other precision is specified.
Definition: Types.h:198
TypeOfCellQuadrilateral.h
CellTraits.h
VecTraits.h
vtkm::ErrorCode::InvalidNumberOfPoints
@ InvalidNumberOfPoints
Generated by   doxygen 1.8.17