Implementation of Taper for Hex and Quad. Untested

vtkm/exec/cellmetrics/CellTaperMetric.h

@@ -41,20 +41,51 @@ namespace vtkm
 {
 namespace exec
 {
+static constexpr FloatType FLOAT_MAX = vtkm::Infinity<FloatType>();
+static constexpr FloatType FLOAT_MIN = vtkm::NegativeInfinity<FloatType>();
+
+template <typename OutType, typename PointCoordVecType, typename CellShapeType>
+VTKM_EXEC OutType CellStretchMetric(const vtkm::IdComponent& numPts,
+                                    const PointCoordVecType& pts,
+                                    CellShapeType shape,
+                                    const vtkm::exec::FunctorBase& worklet)
+{
+    worklet.RaiseError("Shape type template must be specified to compute taper")
+    return OutType(-1.0);
+}
+
 template <typename OutType, typename PointCoordVecType>
 VTKM_EXEC OutType CellTaperMetric(const vtkm::IdComponents &numPts,
                                   const PointCoordVecType& pts,
                                   vtkm::CellShapeHex,
                                   const vtkm::exec::FunctorBase& worklet)
 {
+    FloatDefault X1 = (pts[1] - pts[0]) + (pts[2] - pts[3]) + (pts[5] - pts[4]) + (pts[6] - pts[7]);
+    FloatDefault X2 = (pts[3] - pts[0]) + (pts[2] - pts[1]) + (pts[7] - pts[4]) + (pts[6] - pts[5]);
+    FloatDefault X3 = (pts[4] - pts[0]) + (pts[5] - pts[1]) + (pts[6] - pts[2]) + (pts[7] - pts[3]);
+    if ((X1 < FLOAT_MIN) || (X2 < FLOAT_MIN) || (X3 < FLOAT_MIN))
+	return FLOAT_MAX;
+    FloatDefault X12 = ((pts[2] - pts[3]) - (pts[1] - pts[0])) + ((pts[6] - pts[7]) - (pts[5] - pts[4]));
+    FloatDefault X13 = ((pts[5] - pts[1]) - (pts[4] - pts[0])) + ((pts[6] - pts[2]) - (pts[7] - pts[3]));
+    FloatDefault X23 = ((pts[7] - pts[4]) - (pts[3] - pts[0])) + ((pts[6] - pts[5]) - (pts[2] - pts[1]));
+    FloatDefault T12 = X12 / vtkm::Min(X1,vtkm::Min(X2,X3));
+    FloatDefault T13 = X13 / vtkm::Min(X1,vtkm::Min(X2,X3));
+    FloatDefault T23 = X23 / vtkm::Min(X1,vtkm::Min(X2,X3));
+    return vtkm::Max(T12, vtkm::Max(T13,T23));
 }
 
 template <typename OutType, typename PointCoordVecType>
 VTKM_EXEC OutType CellTaperMetric(const vtkm::IdComponents &numPts,
                                   const PointCoordVecType& pts,
-                                  vtkm::CellShapeHex,
+                                  vtkm::CellShapeQuad,
                                   const vtkm::exec::FunctorBase& worklet)
 {
+    FloatDefault X1 = (pts[1] - pts[0]) + (pts[2] - pts[3]);
+    FloatDefault X2 = (pts[0] - pts[3]) + (pts[2] - pts[1]);
+    if ((X1 < FLOAT_MIN) || (X2 < FLOAT_MIN)) 
+	return FLOAT_MAX;
+    FloatDefault X12 = vtkm::Magnitude((pts[0] - pts[1]) + (pts[2] - pts[3]));
+    return X12 / vtkm::Min(X1,X2);
 }
 
 } // exec