VTK-m  2.0
exec/CellMeasure.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 //
6 // This software is distributed WITHOUT ANY WARRANTY; without even
7 // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
8 // PURPOSE. See the above copyright notice for more information.
9 //============================================================================
10 #ifndef vtk_m_exec_CellMeasure_h
11 #define vtk_m_exec_CellMeasure_h
12 
16 #include <vtkm/CellShape.h>
17 #include <vtkm/CellTraits.h>
18 #include <vtkm/ErrorCode.h>
19 #include <vtkm/VecTraits.h>
20 #include <vtkm/VectorAnalysis.h>
21 #include <vtkm/exec/FunctorBase.h>
22 
23 namespace vtkm
24 {
25 namespace exec
26 {
27 
29 template <typename OutType, typename PointCoordVecType, typename CellShapeType>
30 VTKM_EXEC OutType CellMeasure(const vtkm::IdComponent& numPts,
31  const PointCoordVecType& pts,
32  CellShapeType shape,
33  vtkm::ErrorCode&)
34 {
35  (void)numPts;
36  (void)pts;
37  (void)shape;
38  return OutType(0.0);
39 }
40 
41 // ========================= Arc length cells ==================================
43 template <typename OutType, typename PointCoordVecType>
44 VTKM_EXEC OutType CellMeasure(const vtkm::IdComponent& numPts,
45  const PointCoordVecType& pts,
46  vtkm::CellShapeTagLine,
47  vtkm::ErrorCode& ec)
48 {
49  OutType arcLength(0.0);
50  if (numPts < 2)
51  {
52  ec = vtkm::ErrorCode::InvalidCellMetric;
53  }
54  else
55  {
56  arcLength = static_cast<OutType>(Magnitude(pts[1] - pts[0]));
57  for (int ii = 2; ii < numPts; ++ii)
58  {
59  arcLength += static_cast<OutType>(Magnitude(pts[ii] - pts[ii - 1]));
60  }
61  }
62  return arcLength;
63 }
64 
65 // =============================== Area cells ==================================
67 template <typename OutType, typename PointCoordVecType>
68 VTKM_EXEC OutType CellMeasure(const vtkm::IdComponent& numPts,
69  const PointCoordVecType& pts,
70  vtkm::CellShapeTagTriangle,
71  vtkm::ErrorCode& ec)
72 {
73  if (numPts != 3)
74  {
75  ec = vtkm::ErrorCode::InvalidNumberOfPoints;
76  return OutType(0.0);
77  }
78  typename PointCoordVecType::ComponentType v1 = pts[1] - pts[0];
79  typename PointCoordVecType::ComponentType v2 = pts[2] - pts[0];
80  OutType area = OutType(0.5) * static_cast<OutType>(Magnitude(Cross(v1, v2)));
81  return area;
82 }
83 
85 template <typename OutType, typename PointCoordVecType>
86 VTKM_EXEC OutType CellMeasure(const vtkm::IdComponent& numPts,
87  const PointCoordVecType& pts,
88  vtkm::CellShapeTagQuad,
89  vtkm::ErrorCode& ec)
90 {
91  if (numPts != 4)
92  {
93  ec = vtkm::ErrorCode::InvalidNumberOfPoints;
94  return OutType(0.0);
95  }
96 
97  typename PointCoordVecType::ComponentType edges[4] = {
98  pts[1] - pts[0],
99  pts[2] - pts[1],
100  pts[3] - pts[2],
101  pts[0] - pts[3],
102  };
103 
104  if (vtkm::MagnitudeSquared(edges[0]) == OutType(0.0) ||
105  vtkm::MagnitudeSquared(edges[1]) == OutType(0.0) ||
106  vtkm::MagnitudeSquared(edges[2]) == OutType(0.0) ||
107  vtkm::MagnitudeSquared(edges[3]) == OutType(0.0))
108  return OutType(0.0);
109 
110  typename PointCoordVecType::ComponentType cornerNormals[4] = {
111  Cross(edges[3], edges[0]),
112  Cross(edges[0], edges[1]),
113  Cross(edges[1], edges[2]),
114  Cross(edges[2], edges[3]),
115  };
116 
117  // principal axes
118  typename PointCoordVecType::ComponentType principalAxes[2] = {
119  edges[0] - edges[2],
120  edges[1] - edges[3],
121  };
122 
123  // Unit normal at the quadrilateral center
124  typename PointCoordVecType::ComponentType unitCenterNormal =
125  Cross(principalAxes[0], principalAxes[1]);
126  Normalize(unitCenterNormal);
127 
128  OutType area =
129  static_cast<OutType>(
130  (Dot(unitCenterNormal, cornerNormals[0]) + Dot(unitCenterNormal, cornerNormals[1]) +
131  Dot(unitCenterNormal, cornerNormals[2]) + Dot(unitCenterNormal, cornerNormals[3]))) *
132  OutType(0.25);
133  return area;
134 }
135 
136 template <typename OutType, typename PointCoordVecType>
137 VTKM_EXEC OutType ComputeMeasure(const vtkm::IdComponent&,
138  const PointCoordVecType&,
139  vtkm::CellShapeTagPolygon,
140  vtkm::ErrorCode& ec)
141 {
142  ec = vtkm::ErrorCode::InvalidCellMetric;
143  return OutType(0.0);
144 }
145 
146 
147 // ============================= Volume cells ==================================
149 template <typename OutType, typename PointCoordVecType>
150 VTKM_EXEC OutType CellMeasure(const vtkm::IdComponent& numPts,
151  const PointCoordVecType& pts,
152  vtkm::CellShapeTagTetra,
153  vtkm::ErrorCode& ec)
154 {
155  if (numPts != 4)
156  {
157  ec = vtkm::ErrorCode::InvalidNumberOfPoints;
158  return OutType(0.0);
159  }
160 
161  typename PointCoordVecType::ComponentType v1 = pts[1] - pts[0];
162  typename PointCoordVecType::ComponentType v2 = pts[2] - pts[0];
163  typename PointCoordVecType::ComponentType v3 = pts[3] - pts[0];
164  OutType volume = static_cast<OutType>(Dot(Cross(v1, v2), v3)) / OutType(6.0);
165  return volume;
166 }
167 
169 template <typename OutType, typename PointCoordVecType>
170 VTKM_EXEC OutType CellMeasure(const vtkm::IdComponent& numPts,
171  const PointCoordVecType& pts,
172  vtkm::CellShapeTagHexahedron,
173  vtkm::ErrorCode& ec)
174 {
175  if (numPts != 8)
176  {
177  ec = vtkm::ErrorCode::InvalidNumberOfPoints;
178  return OutType(0.0);
179  }
180 
181  auto efg1 = pts[1];
182  efg1 += pts[2];
183  efg1 += pts[5];
184  efg1 += pts[6];
185  efg1 -= pts[0];
186  efg1 -= pts[3];
187  efg1 -= pts[4];
188  efg1 -= pts[7];
189 
190  auto efg2 = pts[2];
191  efg2 += pts[3];
192  efg2 += pts[6];
193  efg2 += pts[7];
194  efg2 -= pts[0];
195  efg2 -= pts[1];
196  efg2 -= pts[4];
197  efg2 -= pts[5];
198 
199  auto efg3 = pts[4];
200  efg3 += pts[5];
201  efg3 += pts[6];
202  efg3 += pts[7];
203  efg3 -= pts[0];
204  efg3 -= pts[1];
205  efg3 -= pts[2];
206  efg3 -= pts[3];
207 
208  OutType volume = static_cast<OutType>(Dot(Cross(efg2, efg3), efg1)) / OutType(64.0);
209  return volume;
210 }
211 
213 template <typename OutType, typename PointCoordVecType>
214 VTKM_EXEC OutType CellMeasure(const vtkm::IdComponent& numPts,
215  const PointCoordVecType& pts,
216  vtkm::CellShapeTagWedge,
217  vtkm::ErrorCode& ec)
218 {
219  if (numPts != 6)
220  {
221  ec = vtkm::ErrorCode::InvalidNumberOfPoints;
222  return OutType(0.0);
223  }
224 
225  typename PointCoordVecType::ComponentType v0 = pts[1] - pts[0];
226  typename PointCoordVecType::ComponentType v1 = pts[2] - pts[0];
227  typename PointCoordVecType::ComponentType v2 = pts[3] - pts[0];
228  OutType volume = static_cast<OutType>(Dot(Cross(v0, v1), v2)) / OutType(6.0);
229 
230  typename PointCoordVecType::ComponentType v3 = pts[4] - pts[1];
231  typename PointCoordVecType::ComponentType v4 = pts[5] - pts[1];
232  typename PointCoordVecType::ComponentType v5 = pts[3] - pts[1];
233  volume += static_cast<OutType>(Dot(Cross(v3, v4), v5)) / OutType(6.0);
234 
235  typename PointCoordVecType::ComponentType v6 = pts[5] - pts[1];
236  typename PointCoordVecType::ComponentType v7 = pts[2] - pts[1];
237  typename PointCoordVecType::ComponentType v8 = pts[3] - pts[1];
238  volume += static_cast<OutType>(Dot(Cross(v6, v7), v8)) / OutType(6.0);
239 
240  return volume;
241 }
242 
244 template <typename OutType, typename PointCoordVecType>
245 VTKM_EXEC OutType CellMeasure(const vtkm::IdComponent& numPts,
246  const PointCoordVecType& pts,
247  vtkm::CellShapeTagPyramid,
248  vtkm::ErrorCode& ec)
249 {
250  if (numPts != 5)
251  {
252  ec = vtkm::ErrorCode::InvalidNumberOfPoints;
253  return OutType(0.0);
254  }
255 
256  typename PointCoordVecType::ComponentType v1 = pts[1] - pts[0];
257  typename PointCoordVecType::ComponentType v2 = pts[3] - pts[0];
258  typename PointCoordVecType::ComponentType v3 = pts[4] - pts[0];
259  OutType volume = static_cast<OutType>(Dot(Cross(v1, v2), v3)) / OutType(6.0);
260 
261  typename PointCoordVecType::ComponentType v4 = pts[1] - pts[2];
262  typename PointCoordVecType::ComponentType v5 = pts[3] - pts[2];
263  typename PointCoordVecType::ComponentType v6 = pts[4] - pts[2];
264  volume += static_cast<OutType>(Dot(Cross(v5, v4), v6)) / OutType(6.0);
265 
266  return volume;
267 }
268 }
269 }
270 
271 #endif // vtk_m_exec_CellMeasure_h
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::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
CellShape.h
ErrorCode.h
VectorAnalysis.h
vtkm::exec::CellMeasure
VTKM_EXEC OutType CellMeasure(const vtkm::IdComponent &numPts, const PointCoordVecType &pts, CellShapeType shape, vtkm::ErrorCode &)
By default, cells have zero measure unless this template is specialized below.
Definition: exec/CellMeasure.h:30
vtkm::ErrorCode::InvalidCellMetric
@ InvalidCellMetric
FunctorBase.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
vtkm::exec::ComputeMeasure
VTKM_EXEC OutType ComputeMeasure(const vtkm::IdComponent &, const PointCoordVecType &, vtkm::CellShapeTagPolygon, vtkm::ErrorCode &ec)
Definition: exec/CellMeasure.h:137
CellTraits.h
VecTraits.h
vtkm::ErrorCode::InvalidNumberOfPoints
@ InvalidNumberOfPoints
Generated by   doxygen 1.8.17