Viskores  1.0
CellMeasure.h
Go to the documentation of this file.
1 //============================================================================
2 // The contents of this file are covered by the Viskores license. See
3 // LICENSE.txt for details.
4 //
5 // By contributing to this file, all contributors agree to the Developer
6 // Certificate of Origin Version 1.1 (DCO 1.1) as stated in DCO.txt.
7 //============================================================================
8 
9 //============================================================================
10 // Copyright (c) Kitware, Inc.
11 // All rights reserved.
12 // See LICENSE.txt for details.
13 //
14 // This software is distributed WITHOUT ANY WARRANTY; without even
15 // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
16 // PURPOSE. See the above copyright notice for more information.
17 //============================================================================
18 #ifndef viskores_exec_CellMeasure_h
19 #define viskores_exec_CellMeasure_h
20 
24 #include <viskores/CellShape.h>
25 #include <viskores/CellTraits.h>
26 #include <viskores/ErrorCode.h>
27 #include <viskores/VecTraits.h>
30 
31 namespace viskores
32 {
33 namespace exec
34 {
35 
37 template <typename OutType, typename PointCoordVecType, typename CellShapeType>
39  const PointCoordVecType& pts,
40  CellShapeType shape,
42 {
43  (void)numPts;
44  (void)pts;
45  (void)shape;
46  return OutType(0.0);
47 }
48 
49 // ========================= Arc length cells ==================================
51 template <typename OutType, typename PointCoordVecType>
53  const PointCoordVecType& pts,
56 {
57  OutType arcLength(0.0);
58  if (numPts < 2)
59  {
61  }
62  else
63  {
64  arcLength = static_cast<OutType>(Magnitude(pts[1] - pts[0]));
65  for (int ii = 2; ii < numPts; ++ii)
66  {
67  arcLength += static_cast<OutType>(Magnitude(pts[ii] - pts[ii - 1]));
68  }
69  }
70  return arcLength;
71 }
72 
73 // =============================== Area cells ==================================
75 template <typename OutType, typename PointCoordVecType>
77  const PointCoordVecType& pts,
80 {
81  if (numPts != 3)
82  {
84  return OutType(0.0);
85  }
86  typename PointCoordVecType::ComponentType v1 = pts[1] - pts[0];
87  typename PointCoordVecType::ComponentType v2 = pts[2] - pts[0];
88  OutType area = OutType(0.5) * static_cast<OutType>(Magnitude(Cross(v1, v2)));
89  return area;
90 }
91 
93 template <typename OutType, typename PointCoordVecType>
95  const PointCoordVecType& pts,
98 {
99  if (numPts != 4)
100  {
102  return OutType(0.0);
103  }
104 
105  typename PointCoordVecType::ComponentType edges[4] = {
106  pts[1] - pts[0],
107  pts[2] - pts[1],
108  pts[3] - pts[2],
109  pts[0] - pts[3],
110  };
111 
112  if (viskores::MagnitudeSquared(edges[0]) == OutType(0.0) ||
113  viskores::MagnitudeSquared(edges[1]) == OutType(0.0) ||
114  viskores::MagnitudeSquared(edges[2]) == OutType(0.0) ||
115  viskores::MagnitudeSquared(edges[3]) == OutType(0.0))
116  return OutType(0.0);
117 
118  typename PointCoordVecType::ComponentType cornerNormals[4] = {
119  Cross(edges[3], edges[0]),
120  Cross(edges[0], edges[1]),
121  Cross(edges[1], edges[2]),
122  Cross(edges[2], edges[3]),
123  };
124 
125  // principal axes
126  typename PointCoordVecType::ComponentType principalAxes[2] = {
127  edges[0] - edges[2],
128  edges[1] - edges[3],
129  };
130 
131  // Unit normal at the quadrilateral center
132  typename PointCoordVecType::ComponentType unitCenterNormal =
133  Cross(principalAxes[0], principalAxes[1]);
134  Normalize(unitCenterNormal);
135 
136  OutType area =
137  static_cast<OutType>(
138  (Dot(unitCenterNormal, cornerNormals[0]) + Dot(unitCenterNormal, cornerNormals[1]) +
139  Dot(unitCenterNormal, cornerNormals[2]) + Dot(unitCenterNormal, cornerNormals[3]))) *
140  OutType(0.25);
141  return area;
142 }
143 
144 template <typename OutType, typename PointCoordVecType>
146  const PointCoordVecType&,
149 {
151  return OutType(0.0);
152 }
153 
154 
155 // ============================= Volume cells ==================================
157 template <typename OutType, typename PointCoordVecType>
159  const PointCoordVecType& pts,
162 {
163  if (numPts != 4)
164  {
166  return OutType(0.0);
167  }
168 
169  typename PointCoordVecType::ComponentType v1 = pts[1] - pts[0];
170  typename PointCoordVecType::ComponentType v2 = pts[2] - pts[0];
171  typename PointCoordVecType::ComponentType v3 = pts[3] - pts[0];
172  OutType volume = static_cast<OutType>(Dot(Cross(v1, v2), v3)) / OutType(6.0);
173  return volume;
174 }
175 
177 template <typename OutType, typename PointCoordVecType>
179  const PointCoordVecType& pts,
182 {
183  if (numPts != 8)
184  {
186  return OutType(0.0);
187  }
188 
189  auto efg1 = pts[1];
190  efg1 += pts[2];
191  efg1 += pts[5];
192  efg1 += pts[6];
193  efg1 -= pts[0];
194  efg1 -= pts[3];
195  efg1 -= pts[4];
196  efg1 -= pts[7];
197 
198  auto efg2 = pts[2];
199  efg2 += pts[3];
200  efg2 += pts[6];
201  efg2 += pts[7];
202  efg2 -= pts[0];
203  efg2 -= pts[1];
204  efg2 -= pts[4];
205  efg2 -= pts[5];
206 
207  auto efg3 = pts[4];
208  efg3 += pts[5];
209  efg3 += pts[6];
210  efg3 += pts[7];
211  efg3 -= pts[0];
212  efg3 -= pts[1];
213  efg3 -= pts[2];
214  efg3 -= pts[3];
215 
216  OutType volume = static_cast<OutType>(Dot(Cross(efg2, efg3), efg1)) / OutType(64.0);
217  return volume;
218 }
219 
221 template <typename OutType, typename PointCoordVecType>
223  const PointCoordVecType& pts,
226 {
227  if (numPts != 6)
228  {
230  return OutType(0.0);
231  }
232 
233  typename PointCoordVecType::ComponentType v0 = pts[1] - pts[0];
234  typename PointCoordVecType::ComponentType v1 = pts[2] - pts[0];
235  typename PointCoordVecType::ComponentType v2 = pts[3] - pts[0];
236  OutType volume = static_cast<OutType>(Dot(Cross(v0, v1), v2)) / OutType(6.0);
237 
238  typename PointCoordVecType::ComponentType v3 = pts[4] - pts[1];
239  typename PointCoordVecType::ComponentType v4 = pts[5] - pts[1];
240  typename PointCoordVecType::ComponentType v5 = pts[3] - pts[1];
241  volume += static_cast<OutType>(Dot(Cross(v3, v4), v5)) / OutType(6.0);
242 
243  typename PointCoordVecType::ComponentType v6 = pts[5] - pts[1];
244  typename PointCoordVecType::ComponentType v7 = pts[2] - pts[1];
245  typename PointCoordVecType::ComponentType v8 = pts[3] - pts[1];
246  volume += static_cast<OutType>(Dot(Cross(v6, v7), v8)) / OutType(6.0);
247 
248  return volume;
249 }
250 
252 template <typename OutType, typename PointCoordVecType>
254  const PointCoordVecType& pts,
257 {
258  if (numPts != 5)
259  {
261  return OutType(0.0);
262  }
263 
264  typename PointCoordVecType::ComponentType v1 = pts[1] - pts[0];
265  typename PointCoordVecType::ComponentType v2 = pts[3] - pts[0];
266  typename PointCoordVecType::ComponentType v3 = pts[4] - pts[0];
267  OutType volume = static_cast<OutType>(Dot(Cross(v1, v2), v3)) / OutType(6.0);
268 
269  typename PointCoordVecType::ComponentType v4 = pts[1] - pts[2];
270  typename PointCoordVecType::ComponentType v5 = pts[3] - pts[2];
271  typename PointCoordVecType::ComponentType v6 = pts[4] - pts[2];
272  volume += static_cast<OutType>(Dot(Cross(v5, v4), v6)) / OutType(6.0);
273 
274  return volume;
275 }
276 }
277 }
278 
279 #endif // viskores_exec_CellMeasure_h
viskores::CellShapeTagLine
Definition: CellShape.h:158
viskores::CellShapeTagHexahedron
Definition: CellShape.h:167
viskores::Normalize
void Normalize(T &x)
Changes a vector to be normal.
Definition: VectorAnalysis.h:177
viskores::ErrorCode::InvalidCellMetric
@ InvalidCellMetric
A cell metric was requested for a cell that does not support that metric.
viskores::ErrorCode
ErrorCode
Identifies whether an operation was successful or what type of error it had.
Definition: ErrorCode.h:36
viskores::CellShapeTagTriangle
Definition: CellShape.h:160
viskores::ErrorCode::InvalidNumberOfPoints
@ InvalidNumberOfPoints
The wrong number of points was provided for a given cell type.
viskores::IdComponent
viskores::Int32 IdComponent
Base type to use to index small lists.
Definition: Types.h:202
CellShape.h
ErrorCode.h
viskores::Cross
viskores::Vec< typename detail::FloatingPointReturnType< T >::Type, 3 > Cross(const viskores::Vec< T, 3 > &x, const viskores::Vec< T, 3 > &y)
Find the cross product of two vectors.
Definition: VectorAnalysis.h:188
VectorAnalysis.h
viskores::CellShapeTagWedge
Definition: CellShape.h:168
viskores
Groups connected points that have the same field value.
Definition: Atomic.h:27
FunctorBase.h
viskores::CellShapeTagQuad
Definition: CellShape.h:164
viskores::CellShapeTagTetra
Definition: CellShape.h:165
viskores::CellShapeTagPyramid
Definition: CellShape.h:169
viskores::CellShapeTagPolygon
Definition: CellShape.h:162
viskores::Magnitude
detail::FloatingPointReturnType< T >::Type Magnitude(const T &x)
Returns the magnitude of a vector.
Definition: VectorAnalysis.h:108
viskores::MagnitudeSquared
detail::FloatingPointReturnType< T >::Type MagnitudeSquared(const T &x)
Returns the square of the magnitude of a vector.
Definition: VectorAnalysis.h:72
CellTraits.h
viskores::exec::ComputeMeasure
OutType ComputeMeasure(const viskores::IdComponent &, const PointCoordVecType &, viskores::CellShapeTagPolygon, viskores::ErrorCode &ec)
Definition: CellMeasure.h:145
viskores::exec::CellMeasure
OutType CellMeasure(const viskores::IdComponent &numPts, const PointCoordVecType &pts, CellShapeType shape, viskores::ErrorCode &)
By default, cells have zero measure unless this template is specialized below.
Definition: CellMeasure.h:38
VISKORES_EXEC
#define VISKORES_EXEC
Definition: ExportMacros.h:59
VecTraits.h