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 /*!\file
13 * Functions that provide integral measures over cells.
14 */
15
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
28 /// By default, cells have zero measure unless this template is specialized below.
29 template <typename OutType, typename PointCoordVecType, typename CellShapeType>
CellMeasure(const vtkm::IdComponent & numPts,const PointCoordVecType & pts,CellShapeType shape,vtkm::ErrorCode &)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 ==================================
42 /// Compute the arc length of a poly-line cell.
43 template <typename OutType, typename PointCoordVecType>
CellMeasure(const vtkm::IdComponent & numPts,const PointCoordVecType & pts,vtkm::CellShapeTagLine,vtkm::ErrorCode & ec)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 ==================================
66 /// Compute the area of a triangular cell.
67 template <typename OutType, typename PointCoordVecType>
CellMeasure(const vtkm::IdComponent & numPts,const PointCoordVecType & pts,vtkm::CellShapeTagTriangle,vtkm::ErrorCode & ec)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
84 /// Compute the area of a quadrilateral cell.
85 template <typename OutType, typename PointCoordVecType>
CellMeasure(const vtkm::IdComponent & numPts,const PointCoordVecType & pts,vtkm::CellShapeTagQuad,vtkm::ErrorCode & ec)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 typename PointCoordVecType::ComponentType cornerNormals[4] = {
105 Cross(edges[3], edges[0]),
106 Cross(edges[0], edges[1]),
107 Cross(edges[1], edges[2]),
108 Cross(edges[2], edges[3]),
109 };
110
111 // principal axes
112 typename PointCoordVecType::ComponentType principalAxes[2] = {
113 edges[0] - edges[2],
114 edges[1] - edges[3],
115 };
116
117 // Unit normal at the quadrilateral center
118 typename PointCoordVecType::ComponentType unitCenterNormal =
119 Cross(principalAxes[0], principalAxes[1]);
120 Normalize(unitCenterNormal);
121
122 OutType area =
123 (Dot(unitCenterNormal, cornerNormals[0]) + Dot(unitCenterNormal, cornerNormals[1]) +
124 Dot(unitCenterNormal, cornerNormals[2]) + Dot(unitCenterNormal, cornerNormals[3])) *
125 OutType(0.25);
126 return area;
127 }
128
129 template <typename OutType, typename PointCoordVecType>
ComputeMeasure(const vtkm::IdComponent &,const PointCoordVecType &,vtkm::CellShapeTagPolygon,vtkm::ErrorCode & ec)130 VTKM_EXEC OutType ComputeMeasure(const vtkm::IdComponent&,
131 const PointCoordVecType&,
132 vtkm::CellShapeTagPolygon,
133 vtkm::ErrorCode& ec)
134 {
135 ec = vtkm::ErrorCode::InvalidCellMetric;
136 return OutType(0.0);
137 }
138
139
140 // ============================= Volume cells ==================================
141 /// Compute the volume of a tetrahedron.
142 template <typename OutType, typename PointCoordVecType>
CellMeasure(const vtkm::IdComponent & numPts,const PointCoordVecType & pts,vtkm::CellShapeTagTetra,vtkm::ErrorCode & ec)143 VTKM_EXEC OutType CellMeasure(const vtkm::IdComponent& numPts,
144 const PointCoordVecType& pts,
145 vtkm::CellShapeTagTetra,
146 vtkm::ErrorCode& ec)
147 {
148 if (numPts != 4)
149 {
150 ec = vtkm::ErrorCode::InvalidNumberOfPoints;
151 return OutType(0.0);
152 }
153
154 typename PointCoordVecType::ComponentType v1 = pts[1] - pts[0];
155 typename PointCoordVecType::ComponentType v2 = pts[2] - pts[0];
156 typename PointCoordVecType::ComponentType v3 = pts[3] - pts[0];
157 OutType volume = Dot(Cross(v1, v2), v3) / OutType(6.0);
158 return volume;
159 }
160
161 /// Compute the volume of a hexahedral cell (approximated via triple product of average edge along each parametric axis).
162 template <typename OutType, typename PointCoordVecType>
CellMeasure(const vtkm::IdComponent & numPts,const PointCoordVecType & pts,vtkm::CellShapeTagHexahedron,vtkm::ErrorCode & ec)163 VTKM_EXEC OutType CellMeasure(const vtkm::IdComponent& numPts,
164 const PointCoordVecType& pts,
165 vtkm::CellShapeTagHexahedron,
166 vtkm::ErrorCode& ec)
167 {
168 if (numPts != 8)
169 {
170 ec = vtkm::ErrorCode::InvalidNumberOfPoints;
171 return OutType(0.0);
172 }
173
174 auto efg1 = pts[1];
175 efg1 += pts[2];
176 efg1 += pts[5];
177 efg1 += pts[6];
178 efg1 -= pts[0];
179 efg1 -= pts[3];
180 efg1 -= pts[4];
181 efg1 -= pts[7];
182
183 auto efg2 = pts[2];
184 efg2 += pts[3];
185 efg2 += pts[6];
186 efg2 += pts[7];
187 efg2 -= pts[0];
188 efg2 -= pts[1];
189 efg2 -= pts[4];
190 efg2 -= pts[5];
191
192 auto efg3 = pts[4];
193 efg3 += pts[5];
194 efg3 += pts[6];
195 efg3 += pts[7];
196 efg3 -= pts[0];
197 efg3 -= pts[1];
198 efg3 -= pts[2];
199 efg3 -= pts[3];
200
201 OutType volume = Dot(Cross(efg2, efg3), efg1) / OutType(64.0);
202 return volume;
203 }
204
205 /// Compute the volume of a wedge cell (approximated as 3 tetrahedra).
206 template <typename OutType, typename PointCoordVecType>
CellMeasure(const vtkm::IdComponent & numPts,const PointCoordVecType & pts,vtkm::CellShapeTagWedge,vtkm::ErrorCode & ec)207 VTKM_EXEC OutType CellMeasure(const vtkm::IdComponent& numPts,
208 const PointCoordVecType& pts,
209 vtkm::CellShapeTagWedge,
210 vtkm::ErrorCode& ec)
211 {
212 if (numPts != 6)
213 {
214 ec = vtkm::ErrorCode::InvalidNumberOfPoints;
215 return OutType(0.0);
216 }
217
218 typename PointCoordVecType::ComponentType v0 = pts[1] - pts[0];
219 typename PointCoordVecType::ComponentType v1 = pts[2] - pts[0];
220 typename PointCoordVecType::ComponentType v2 = pts[3] - pts[0];
221 OutType volume = Dot(Cross(v0, v1), v2) / OutType(6.0);
222
223 typename PointCoordVecType::ComponentType v3 = pts[4] - pts[1];
224 typename PointCoordVecType::ComponentType v4 = pts[5] - pts[1];
225 typename PointCoordVecType::ComponentType v5 = pts[3] - pts[1];
226 volume += Dot(Cross(v3, v4), v5) / OutType(6.0);
227
228 typename PointCoordVecType::ComponentType v6 = pts[5] - pts[1];
229 typename PointCoordVecType::ComponentType v7 = pts[2] - pts[1];
230 typename PointCoordVecType::ComponentType v8 = pts[3] - pts[1];
231 volume += Dot(Cross(v6, v7), v8) / OutType(6.0);
232
233 return volume;
234 }
235
236 /// Compute the volume of a pyramid (approximated as 2 tetrahedra)
237 template <typename OutType, typename PointCoordVecType>
CellMeasure(const vtkm::IdComponent & numPts,const PointCoordVecType & pts,vtkm::CellShapeTagPyramid,vtkm::ErrorCode & ec)238 VTKM_EXEC OutType CellMeasure(const vtkm::IdComponent& numPts,
239 const PointCoordVecType& pts,
240 vtkm::CellShapeTagPyramid,
241 vtkm::ErrorCode& ec)
242 {
243 if (numPts != 5)
244 {
245 ec = vtkm::ErrorCode::InvalidNumberOfPoints;
246 return OutType(0.0);
247 }
248
249 typename PointCoordVecType::ComponentType v1 = pts[1] - pts[0];
250 typename PointCoordVecType::ComponentType v2 = pts[3] - pts[0];
251 typename PointCoordVecType::ComponentType v3 = pts[4] - pts[0];
252 OutType volume = Dot(Cross(v1, v2), v3) / OutType(6.0);
253
254 typename PointCoordVecType::ComponentType v4 = pts[1] - pts[2];
255 typename PointCoordVecType::ComponentType v5 = pts[3] - pts[2];
256 typename PointCoordVecType::ComponentType v6 = pts[4] - pts[2];
257 volume += Dot(Cross(v5, v4), v6) / OutType(6.0);
258
259 return volume;
260 }
261 }
262 }
263
264 #endif // vtk_m_exec_CellMeasure_h
265