Compute cell measures (arc length, area, volume).

This commit adds a worklet and a filter to compute
signed integral measures of cells.

It also begins the practic of adding a markdown
changelog entry in `doc/changelog`. See the
changelog entry for details.

docs/changelog/cell-measure.md

@@ -0,0 +1,54 @@
+# Cell measure functions, worklet, and filter
+
+VTK-m now provides free functions, a worklet, and a filter for computing
+the integral measure of a cell (i.e., its arc length, area, or volume).
+
+The free functions are located in `vtkm/exec/CellMeasure.h` and share the
+same signature:
+
+```c++
+  template<typename OutType, typename PointVecType>
+  OutType CellMeasure(
+    const vtkm::IdComponent& numPts,
+    const PointCoordVecType& pts,
+    CellShapeTag,
+    const vtkm::exec::FunctorBase& worklet);
+```
+
+The number of points argument is provided for cell-types such as lines,
+which allow an arbitrary number of points per cell.
+See the worklet for examples of their use.
+
+The worklet is named `vtkm::worklet::CellMeasure` and takes a template
+parameter that is a tag list of measures to include.
+Cells that are not selected by the tag list return a measure of 0.
+Some convenient tag lists are predefined for you:
+
++ `vtkm::ArcLength` will only compute the measure of cells with a 1-dimensional parameter-space.
++ `vtkm::Area` will only compute the measure of cells with a 2-dimensional parameter-space.
++ `vtkm::Volume` will only compute the measure of cells with a 3-dimensional parameter-space.
++ `vtkm::AllMeasures` will compute all of the above.
+
+The filter version, named `vtkm::filter::CellMeasures` – plural since
+it produces a cell-centered array of measures — takes the same template
+parameter and tag lists as the worklet.
+By default, the output array of measure values is named "measure" but
+the filter accepts other names via the `SetCellMeasureName()` method.
+
+The only cell type that is not supported is the polygon;
+you must triangulate polygons before running this filter.
+See the unit tests for examples of how to use the worklet and filter.
+
+The cell measures are all signed: negative measures indicate that the cell is inverted.
+Simplicial cells (points, lines, triangles, tetrahedra) cannot not be inverted
+by definition and thus always return values above or equal to 0.0.
+Negative values indicate either the order in which vertices appear in its connectivity
+array is improper or the relative locations of the vertices in world coordinates
+result in a cell with a negative Jacobian somewhere in its interior.
+Finaly, note that cell measures may return invalid (NaN) or infinite (Inf, -Inf)
+values if the cell is poorly defined, e.g., has coincident vertices
+or a parametric dimension larger than the space spanned by its world-coordinate
+vertices.
+
+The verdict mesh quality library was used as the source of the methods
+for approximating the cell measures.

vtkm/exec/CMakeLists.txt

@@ -25,6 +25,7 @@ set(headers
   CellFace.h
   CellInside.h
   CellInterpolate.h
+  CellMeasure.h
   ColorTable.h
   ConnectivityExplicit.h
   ConnectivityPermuted.h

vtkm/exec/CellMeasure.h

@@ -0,0 +1,267 @@
+//============================================================================
+//  Copyright (c) Kitware, Inc.
+//  All rights reserved.
+//  See LICENSE.txt for details.
+//  This software is distributed WITHOUT ANY WARRANTY; without even
+//  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+//  PURPOSE.  See the above copyright notice for more information.
+//
+//  Copyright 2014 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
+//  Copyright 2014 UT-Battelle, LLC.
+//  Copyright 2014 Los Alamos National Security.
+//
+//  Under the terms of Contract DE-NA0003525 with NTESS,
+//  the U.S. Government retains certain rights in this software.
+//
+//  Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
+//  Laboratory (LANL), the U.S. Government retains certain rights in
+//  this software.
+//============================================================================
+#ifndef vtk_m_exec_CellMeasure_h
+#define vtk_m_exec_CellMeasure_h
+/*!\file Functions that provide integral measures over cells. */
+
+#include "vtkm/CellShape.h"
+#include "vtkm/CellTraits.h"
+#include "vtkm/VecTraits.h"
+#include "vtkm/VectorAnalysis.h"
+#include "vtkm/exec/FunctorBase.h"
+
+namespace vtkm
+{
+namespace exec
+{
+
+/// By default, cells have zero measure unless this template is specialized below.
+template <typename OutType, typename PointCoordVecType, typename CellShapeType>
+VTKM_EXEC OutType CellMeasure(const vtkm::IdComponent& numPts,
+                              const PointCoordVecType& pts,
+                              CellShapeType shape,
+                              const vtkm::exec::FunctorBase&)
+{
+  (void)numPts;
+  (void)pts;
+  (void)shape;
+  return OutType(0.0);
+}
+
+// ========================= Arc length cells ==================================
+/// Compute the arc length of a poly-line cell.
+template <typename OutType, typename PointCoordVecType>
+VTKM_EXEC OutType CellMeasure(const vtkm::IdComponent& numPts,
+                              const PointCoordVecType& pts,
+                              vtkm::CellShapeTagLine,
+                              const vtkm::exec::FunctorBase& worklet)
+{
+  OutType arcLength(0.0);
+  if (numPts < 2)
+  {
+    worklet.RaiseError("Degenerate line has no arc length.");
+  }
+  else
+  {
+    arcLength = Magnitude(pts[1] - pts[0]);
+    for (int ii = 2; ii < numPts; ++ii)
+    {
+      arcLength += Magnitude(pts[ii] - pts[ii - 1]);
+    }
+  }
+  return arcLength;
+}
+
+// =============================== Area cells ==================================
+/// Compute the area of a triangular cell.
+template <typename OutType, typename PointCoordVecType>
+VTKM_EXEC OutType CellMeasure(const vtkm::IdComponent& numPts,
+                              const PointCoordVecType& pts,
+                              vtkm::CellShapeTagTriangle,
+                              const vtkm::exec::FunctorBase& worklet)
+{
+  if (numPts != 3)
+  {
+    worklet.RaiseError("Area(triangle) requires 3 points.");
+    return OutType(0.0);
+  }
+  typename PointCoordVecType::ComponentType v1 = pts[1] - pts[0];
+  typename PointCoordVecType::ComponentType v2 = pts[2] - pts[0];
+  OutType area = OutType(0.5) * Magnitude(Cross(v1, v2));
+  return area;
+}
+
+/// Compute the area of a quadrilateral cell.
+template <typename OutType, typename PointCoordVecType>
+VTKM_EXEC OutType CellMeasure(const vtkm::IdComponent& numPts,
+                              const PointCoordVecType& pts,
+                              vtkm::CellShapeTagQuad,
+                              const vtkm::exec::FunctorBase& worklet)
+{
+  if (numPts != 4)
+  {
+    worklet.RaiseError("Area(quad) requires 4 points.");
+    return OutType(0.0);
+  }
+
+  typename PointCoordVecType::ComponentType edges[4] = {
+    pts[1] - pts[0], pts[2] - pts[1], pts[3] - pts[2], pts[0] - pts[3],
+  };
+
+  typename PointCoordVecType::ComponentType cornerNormals[4] = {
+    Cross(edges[3], edges[0]),
+    Cross(edges[0], edges[1]),
+    Cross(edges[1], edges[2]),
+    Cross(edges[2], edges[3]),
+  };
+
+  // principal axes
+  typename PointCoordVecType::ComponentType principalAxes[2] = {
+    edges[0] - edges[2], edges[1] - edges[3],
+  };
+
+  // Unit normal at the quadrilateral center
+  typename PointCoordVecType::ComponentType unitCenterNormal =
+    Cross(principalAxes[0], principalAxes[1]);
+  Normalize(unitCenterNormal);
+
+  OutType area =
+    (dot(unitCenterNormal, cornerNormals[0]) + dot(unitCenterNormal, cornerNormals[1]) +
+     dot(unitCenterNormal, cornerNormals[2]) + dot(unitCenterNormal, cornerNormals[3])) *
+    OutType(0.25);
+  return area;
+}
+
+template <typename OutType, typename PointCoordVecType>
+VTKM_EXEC OutType ComputeMeasure(const vtkm::IdComponent&,
+                                 const PointCoordVecType&,
+                                 vtkm::CellShapeTagPolygon,
+                                 const vtkm::exec::FunctorBase& worklet)
+{
+  worklet.RaiseError("CellMeasure does not support area computation for polygons.");
+  return OutType(0.0);
+}
+
+
+// ============================= Volume cells ==================================
+/// Compute the volume of a tetrahedron.
+template <typename OutType, typename PointCoordVecType>
+VTKM_EXEC OutType CellMeasure(const vtkm::IdComponent& numPts,
+                              const PointCoordVecType& pts,
+                              vtkm::CellShapeTagTetra,
+                              const vtkm::exec::FunctorBase& worklet)
+{
+  if (numPts != 4)
+  {
+    worklet.RaiseError("Volume(tetrahedron) requires 4 points.");
+    return OutType(0.0);
+  }
+
+  typename PointCoordVecType::ComponentType v1 = pts[1] - pts[0];
+  typename PointCoordVecType::ComponentType v2 = pts[2] - pts[0];
+  typename PointCoordVecType::ComponentType v3 = pts[3] - pts[0];
+  OutType volume = dot(Cross(v1, v2), v3) / OutType(6.0);
+  return volume;
+}
+
+/// Compute the volume of a hexahedral cell (approximated via triple product of average edge along each parametric axis).
+template <typename OutType, typename PointCoordVecType>
+VTKM_EXEC OutType CellMeasure(const vtkm::IdComponent& numPts,
+                              const PointCoordVecType& pts,
+                              vtkm::CellShapeTagHexahedron,
+                              const vtkm::exec::FunctorBase& worklet)
+{
+  if (numPts != 8)
+  {
+    worklet.RaiseError("Volume(hexahedron) requires 8 points.");
+    return OutType(0.0);
+  }
+
+  auto efg1 = pts[1];
+  efg1 += pts[2];
+  efg1 += pts[5];
+  efg1 += pts[6];
+  efg1 -= pts[0];
+  efg1 -= pts[3];
+  efg1 -= pts[4];
+  efg1 -= pts[7];
+
+  auto efg2 = pts[2];
+  efg2 += pts[3];
+  efg2 += pts[6];
+  efg2 += pts[7];
+  efg2 -= pts[0];
+  efg2 -= pts[1];
+  efg2 -= pts[4];
+  efg2 -= pts[5];
+
+  auto efg3 = pts[4];
+  efg3 += pts[5];
+  efg3 += pts[6];
+  efg3 += pts[7];
+  efg3 -= pts[0];
+  efg3 -= pts[1];
+  efg3 -= pts[2];
+  efg3 -= pts[3];
+
+  OutType volume = dot(Cross(efg2, efg3), efg1) / OutType(64.0);
+  return volume;
+}
+
+/// Compute the volume of a wedge cell (approximated as 3 tetrahedra).
+template <typename OutType, typename PointCoordVecType>
+VTKM_EXEC OutType CellMeasure(const vtkm::IdComponent& numPts,
+                              const PointCoordVecType& pts,
+                              vtkm::CellShapeTagWedge,
+                              const vtkm::exec::FunctorBase& worklet)
+{
+  if (numPts != 6)
+  {
+    worklet.RaiseError("Volume(wedge) requires 6 points.");
+    return OutType(0.0);
+  }
+
+  typename PointCoordVecType::ComponentType v0 = pts[1] - pts[0];
+  typename PointCoordVecType::ComponentType v1 = pts[2] - pts[0];
+  typename PointCoordVecType::ComponentType v2 = pts[3] - pts[0];
+  OutType volume = dot(Cross(v0, v1), v2) / OutType(6.0);
+
+  typename PointCoordVecType::ComponentType v3 = pts[4] - pts[1];
+  typename PointCoordVecType::ComponentType v4 = pts[5] - pts[1];
+  typename PointCoordVecType::ComponentType v5 = pts[3] - pts[1];
+  volume += dot(Cross(v3, v4), v5) / OutType(6.0);
+
+  typename PointCoordVecType::ComponentType v6 = pts[5] - pts[1];
+  typename PointCoordVecType::ComponentType v7 = pts[2] - pts[1];
+  typename PointCoordVecType::ComponentType v8 = pts[3] - pts[1];
+  volume += dot(Cross(v6, v7), v8) / OutType(6.0);
+
+  return volume;
+}
+
+/// Compute the volume of a pyramid (approximated as 2 tetrahedra)
+template <typename OutType, typename PointCoordVecType>
+VTKM_EXEC OutType CellMeasure(const vtkm::IdComponent& numPts,
+                              const PointCoordVecType& pts,
+                              vtkm::CellShapeTagPyramid,
+                              const vtkm::exec::FunctorBase& worklet)
+{
+  if (numPts != 5)
+  {
+    worklet.RaiseError("Volume(pyramid) requires 5 points.");
+    return OutType(0.0);
+  }
+
+  typename PointCoordVecType::ComponentType v1 = pts[1] - pts[0];
+  typename PointCoordVecType::ComponentType v2 = pts[3] - pts[0];
+  typename PointCoordVecType::ComponentType v3 = pts[4] - pts[0];
+  OutType volume = dot(Cross(v1, v2), v3) / OutType(6.0);
+
+  typename PointCoordVecType::ComponentType v4 = pts[1] - pts[2];
+  typename PointCoordVecType::ComponentType v5 = pts[3] - pts[2];
+  typename PointCoordVecType::ComponentType v6 = pts[4] - pts[2];
+  volume += dot(Cross(v5, v4), v6) / OutType(6.0);
+
+  return volume;
+}
+}
+}
+
+#endif // vtk_m_exec_CellMeasure_h

vtkm/filter/CMakeLists.txt

@@ -20,6 +20,7 @@
 
 set(headers
   CellAverage.h
+  CellMeasures.h
   CleanGrid.h
   ClipWithField.h
   ClipWithImplicitFunction.h
@@ -64,6 +65,7 @@ set(headers
 
 set(header_template_sources
   CellAverage.hxx
+  CellMeasures.hxx
   CleanGrid.hxx
   ClipWithField.hxx
   ClipWithImplicitFunction.hxx

vtkm/filter/CellMeasures.h

@@ -0,0 +1,70 @@
+//============================================================================
+//  Copyright (c) Kitware, Inc.
+//  All rights reserved.
+//  See LICENSE.txt for details.
+//  This software is distributed WITHOUT ANY WARRANTY; without even
+//  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+//  PURPOSE.  See the above copyright notice for more information.
+//
+//  Copyright 2014 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
+//  Copyright 2014 UT-Battelle, LLC.
+//  Copyright 2014 Los Alamos National Security.
+//
+//  Under the terms of Contract DE-NA0003525 with NTESS,
+//  the U.S. Government retains certain rights in this software.
+//
+//  Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
+//  Laboratory (LANL), the U.S. Government retains certain rights in
+//  this software.
+//============================================================================
+
+#ifndef vtk_m_filter_CellMeasures_h
+#define vtk_m_filter_CellMeasures_h
+
+#include "vtkm/filter/FilterCell.h"
+#include "vtkm/worklet/CellMeasure.h"
+
+namespace vtkm
+{
+namespace filter
+{
+
+/// \brief Compute the measure of each (3D) cell in a dataset.
+///
+/// CellMeasures is a filter that generates a new cell data array (i.e., one value
+/// specified per cell) holding the signed measure of the cell
+/// or 0 (if measure is not well defined or the cell type is unsupported).
+///
+/// By default, the new cell-data array is named "measure".
+template <typename IntegrationType>
+class CellMeasures : public vtkm::filter::FilterCell<CellMeasures<IntegrationType>>
+{
+public:
+  VTKM_CONT
+  CellMeasures();
+
+  /// Set/Get the name of the cell measure field. If empty, "measure" is used.
+  void SetCellMeasureName(const std::string& name) { this->SetOutputFieldName(name); }
+  const std::string& GetCellMeasureName() const { return this->GetOutputFieldName(); }
+
+  template <typename T, typename StorageType, typename DerivedPolicy, typename DeviceAdapter>
+  VTKM_CONT vtkm::cont::DataSet DoExecute(
+    const vtkm::cont::DataSet& input,
+    const vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>, StorageType>& points,
+    const vtkm::filter::FieldMetadata& fieldMeta,
+    const vtkm::filter::PolicyBase<DerivedPolicy>& policy,
+    const DeviceAdapter& tag);
+};
+
+template <typename IntegrationType>
+class FilterTraits<CellMeasures<IntegrationType>>
+{
+public:
+  using InputFieldTypeList = vtkm::TypeListTagFieldVec3;
+};
+}
+} // namespace vtkm::filter
+
+#include <vtkm/filter/CellMeasures.hxx>
+
+#endif // vtk_m_filter_CellMeasures_h

vtkm/filter/CellMeasures.hxx

@@ -0,0 +1,69 @@
+//============================================================================
+//  Copyright (c) Kitware, Inc.
+//  All rights reserved.
+//  See LICENSE.txt for details.
+//  This software is distributed WITHOUT ANY WARRANTY; without even
+//  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+//  PURPOSE.  See the above copyright notice for more information.
+//
+//  Copyright 2014 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
+//  Copyright 2014 UT-Battelle, LLC.
+//  Copyright 2014 Los Alamos National Security.
+//
+//  Under the terms of Contract DE-NA0003525 with NTESS,
+//  the U.S. Government retains certain rights in this software.
+//
+//  Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
+//  Laboratory (LANL), the U.S. Government retains certain rights in
+//  this software.
+//============================================================================
+
+#include "vtkm/cont/DynamicCellSet.h"
+#include "vtkm/cont/ErrorFilterExecution.h"
+#include "vtkm/filter/internal/CreateResult.h"
+#include "vtkm/worklet/DispatcherMapTopology.h"
+
+namespace vtkm
+{
+namespace filter
+{
+
+//-----------------------------------------------------------------------------
+template <typename IntegrationType>
+inline VTKM_CONT CellMeasures<IntegrationType>::CellMeasures()
+  : vtkm::filter::FilterCell<CellMeasures<IntegrationType>>()
+{
+  this->SetUseCoordinateSystemAsField(true);
+}
+
+//-----------------------------------------------------------------------------
+template <typename IntegrationType>
+template <typename T, typename StorageType, typename DerivedPolicy, typename DeviceAdapter>
+inline VTKM_CONT vtkm::cont::DataSet CellMeasures<IntegrationType>::DoExecute(
+  const vtkm::cont::DataSet& input,
+  const vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>, StorageType>& points,
+  const vtkm::filter::FieldMetadata& fieldMeta,
+  const vtkm::filter::PolicyBase<DerivedPolicy>& policy,
+  const DeviceAdapter&)
+{
+  VTKM_ASSERT(fieldMeta.IsPointField());
+  const auto& cellset = input.GetCellSet(this->GetActiveCellSetIndex());
+  vtkm::cont::ArrayHandle<T> outArray;
+
+  vtkm::worklet::DispatcherMapTopology<vtkm::worklet::CellMeasure<IntegrationType>, DeviceAdapter>
+    dispatcher;
+  dispatcher.Invoke(vtkm::filter::ApplyPolicy(cellset, policy), points, outArray);
+
+  vtkm::cont::DataSet result;
+  std::string outputName = this->GetCellMeasureName();
+  if (outputName.empty())
+  {
+    // Default name is name of input.
+    outputName = "measure";
+  }
+  result = internal::CreateResult(input, outArray, outputName, vtkm::cont::Field::ASSOC_CELL_SET);
+
+  return result;
+}
+}
+} // namespace vtkm::filter

vtkm/filter/SurfaceNormals.h

@@ -85,8 +85,7 @@ template <>
 class FilterTraits<SurfaceNormals>
 {
 public:
-  using InputFieldTypeList =
-    vtkm::ListTagBase<vtkm::Vec<vtkm::Float32, 3>, vtkm::Vec<vtkm::Float64, 3>>;
+  using InputFieldTypeList = vtkm::TypeListTagFieldVec3;
 };
 }
 } // vtkm::filter

vtkm/filter/testing/CMakeLists.txt

@@ -20,6 +20,7 @@
 
 set(unit_tests
   UnitTestCellAverageFilter.cxx
+  UnitTestCellMeasuresFilter.cxx
   UnitTestCleanGrid.cxx
   UnitTestClipWithFieldFilter.cxx
   UnitTestClipWithImplicitFunctionFilter.cxx

vtkm/filter/testing/UnitTestCellMeasuresFilter.cxx

@@ -0,0 +1,109 @@
+//============================================================================
+//  Copyright (c) Kitware, Inc.
+//  All rights reserved.
+//  See LICENSE.txt for details.
+//  This software is distributed WITHOUT ANY WARRANTY; without even
+//  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+//  PURPOSE.  See the above copyright notice for more information.
+//
+//  Copyright 2017 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
+//  Copyright 2017 UT-Battelle, LLC.
+//  Copyright 2017 Los Alamos National Security.
+//
+//  Under the terms of Contract DE-NA0003525 with NTESS,
+//  the U.S. Government retains certain rights in this software.
+//
+//  Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
+//  Laboratory (LANL), the U.S. Government retains certain rights in
+//  this software.
+//============================================================================
+
+#include "vtkm/filter/CellMeasures.h"
+
+#include "vtkm/cont/DynamicArrayHandle.h"
+#include "vtkm/cont/testing/MakeTestDataSet.h"
+#include "vtkm/cont/testing/Testing.h"
+
+#include <vector>
+
+namespace
+{
+
+template <typename IntegrationType>
+void TestCellMeasuresFilter(vtkm::cont::DataSet& dataset,
+                            const char* msg,
+                            const std::vector<vtkm::Float32>& expected,
+                            const IntegrationType&)
+{
+  std::cout << "Testing CellMeasures Filter on " << msg << "\n";
+
+  vtkm::filter::CellMeasures<IntegrationType> vols;
+  vtkm::cont::DataSet outputData = vols.Execute(dataset);
+
+  VTKM_TEST_ASSERT(vols.GetCellMeasureName().empty(), "Default output field name should be empty.");
+  VTKM_TEST_ASSERT(outputData.GetNumberOfCellSets() == 1,
+                   "Wrong number of cellsets in the output dataset");
+  VTKM_TEST_ASSERT(outputData.GetNumberOfCoordinateSystems() == 1,
+                   "Wrong number of coordinate systems in the output dataset");
+  VTKM_TEST_ASSERT(outputData.GetCellSet().GetNumberOfCells() ==
+                     static_cast<vtkm::Id>(expected.size()),
+                   "Wrong number of cells in the output dataset");
+
+  // Check that the empty measure name above produced a field with the expected name.
+  vols.SetCellMeasureName("measure");
+  vtkm::cont::DynamicArrayHandle temp = outputData.GetField(vols.GetCellMeasureName()).GetData();
+  VTKM_TEST_ASSERT(temp.GetNumberOfValues() == static_cast<vtkm::Id>(expected.size()),
+                   "Output field could not be found or was improper.");
+
+  vtkm::cont::ArrayHandle<vtkm::Float32> resultArrayHandle;
+  temp.CopyTo(resultArrayHandle);
+  VTKM_TEST_ASSERT(resultArrayHandle.GetNumberOfValues() == static_cast<vtkm::Id>(expected.size()),
+                   "Wrong number of entries in the output dataset");
+
+  for (unsigned int i = 0; i < static_cast<unsigned int>(expected.size()); ++i)
+  {
+    VTKM_TEST_ASSERT(
+      test_equal(resultArrayHandle.GetPortalConstControl().Get(vtkm::Id(i)), expected[i]),
+      "Wrong result for CellMeasure filter");
+  }
+}
+
+void TestCellMeasures()
+{
+  using vtkm::Volume;
+  using vtkm::AllMeasures;
+
+  vtkm::cont::testing::MakeTestDataSet factory;
+  vtkm::cont::DataSet data;
+
+  data = factory.Make3DExplicitDataSet2();
+  TestCellMeasuresFilter(data, "explicit dataset 2", { -1.f }, AllMeasures());
+
+  data = factory.Make3DExplicitDataSet3();
+  TestCellMeasuresFilter(data, "explicit dataset 3", { -1.f / 6.f }, AllMeasures());
+
+  data = factory.Make3DExplicitDataSet4();
+  TestCellMeasuresFilter(data, "explicit dataset 4", { -1.f, -1.f }, AllMeasures());
+
+  data = factory.Make3DExplicitDataSet5();
+  TestCellMeasuresFilter(
+    data, "explicit dataset 5", { 1.f, 1.f / 3.f, 1.f / 6.f, -1.f / 2.f }, AllMeasures());
+
+  data = factory.Make3DExplicitDataSet6();
+  TestCellMeasuresFilter(data,
+                         "explicit dataset 6 (only volume)",
+                         { 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.083426f, 0.25028f },
+                         Volume());
+  TestCellMeasuresFilter(
+    data,
+    "explicit dataset 6 (all)",
+    { 0.999924f, 0.999924f, 0.f, 0.f, 3.85516f, 1.00119f, 0.083426f, 0.25028f },
+    AllMeasures());
+}
+
+} // anonymous namespace
+
+int UnitTestCellMeasuresFilter(int, char* [])
+{
+  return vtkm::cont::testing::Testing::Run(TestCellMeasures);
+}

vtkm/worklet/CMakeLists.txt

@@ -22,6 +22,7 @@ set(headers
   AverageByKey.h
   CellAverage.h
   CellDeepCopy.h
+  CellMeasure.h
   Clip.h
   ContourTreeUniform.h
   CosmoTools.h

vtkm/worklet/CellMeasure.h

@@ -0,0 +1,151 @@
+//============================================================================
+//  Copyright (c) Kitware, Inc.
+//  All rights reserved.
+//  See LICENSE.txt for details.
+//  This software is distributed WITHOUT ANY WARRANTY; without even
+//  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+//  PURPOSE.  See the above copyright notice for more information.
+//
+//  Copyright 2014 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
+//  Copyright 2014 UT-Battelle, LLC.
+//  Copyright 2014 Los Alamos National Security.
+//
+//  Under the terms of Contract DE-NA0003525 with NTESS,
+//  the U.S. Government retains certain rights in this software.
+//
+//  Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
+//  Laboratory (LANL), the U.S. Government retains certain rights in
+//  this software.
+//============================================================================
+
+#ifndef vtk_m_worklet_CellMeasure_h
+#define vtk_m_worklet_CellMeasure_h
+
+#include "vtkm/worklet/WorkletMapTopology.h"
+
+#include "vtkm/exec/CellMeasure.h"
+
+namespace vtkm
+{
+
+// Tags used to choose which types of integration are performed on cells:
+struct IntegrateOver
+{
+};
+struct IntegrateOverCurve : IntegrateOver
+{
+};
+struct IntegrateOverSurface : IntegrateOver
+{
+};
+struct IntegrateOverSolid : IntegrateOver
+{
+};
+
+// Lists of acceptable types of integration
+struct ArcLength : vtkm::ListTagBase<IntegrateOverCurve>
+{
+};
+struct Area : vtkm::ListTagBase<IntegrateOverSurface>
+{
+};
+struct Volume : vtkm::ListTagBase<IntegrateOverSolid>
+{
+};
+struct AllMeasures : vtkm::ListTagBase<IntegrateOverSolid, IntegrateOverSurface, IntegrateOverCurve>
+{
+};
+
+namespace worklet
+{
+
+/**\brief Simple functor that returns the spatial integral of each cell as a cell field.
+  *
+  * The integration is done over the spatial extent of the cell and thus units
+  * are either null, arc length, area, or volume depending on whether the parametric
+  * dimension of the cell is 0 (vertices), 1 (curves), 2 (surfaces), or 3 (volumes).
+  * The template parameter of this class configures which types of cells (based on their
+  * parametric dimensions) should be integrated. Other cells will report 0.
+  *
+  * Note that the integrals are signed; inverted cells will report negative values.
+  */
+template <typename IntegrationTypeList>
+class CellMeasure : public vtkm::worklet::WorkletMapPointToCell
+{
+public:
+  typedef void ControlSignature(CellSetIn cellset,
+                                FieldInPoint<Vec3> pointCoords,
+                                FieldOutCell<Scalar> volumesOut);
+  typedef void ExecutionSignature(CellShape, PointCount, _2, _3);
+  using InputDomain = _1;
+
+  template <typename CellShape, typename PointCoordVecType, typename OutType>
+  VTKM_EXEC void operator()(CellShape shape,
+                            const vtkm::IdComponent& numPoints,
+                            const PointCoordVecType& pts,
+                            OutType& volume) const
+  {
+    switch (shape.Id)
+    {
+      vtkmGenericCellShapeMacro(volume =
+                                  this->ComputeMeasure<OutType>(numPoints, pts, CellShapeTag()));
+      default:
+        this->RaiseError("Asked for volume of unknown cell shape.");
+        volume = OutType(0.0);
+    }
+  }
+
+protected:
+  template <typename OutType, typename PointCoordVecType, typename CellShapeType>
+  VTKM_EXEC OutType ComputeMeasure(const vtkm::IdComponent& numPts,
+                                   const PointCoordVecType& pts,
+                                   CellShapeType) const
+  {
+#if defined(VTKM_MSVC)
+#pragma warning(push)
+#pragma warning(disable : 4068) //unknown pragma
+#endif
+#ifdef __NVCC__
+#pragma push
+#pragma diag_suppress = code_is_unreachable
+#endif
+    switch (vtkm::CellTraits<CellShapeType>::TOPOLOGICAL_DIMENSIONS)
+    {
+      case 0:
+        // Fall through to return 0 measure.
+        break;
+      case 1:
+        if (vtkm::ListContains<IntegrationTypeList, IntegrateOverCurve>::value)
+        {
+          return vtkm::exec::CellMeasure<OutType>(numPts, pts, CellShapeType(), *this);
+        }
+        break;
+      case 2:
+        if (vtkm::ListContains<IntegrationTypeList, IntegrateOverSurface>::value)
+        {
+          return vtkm::exec::CellMeasure<OutType>(numPts, pts, CellShapeType(), *this);
+        }
+        break;
+      case 3:
+        if (vtkm::ListContains<IntegrationTypeList, IntegrateOverSolid>::value)
+        {
+          return vtkm::exec::CellMeasure<OutType>(numPts, pts, CellShapeType(), *this);
+        }
+        break;
+      default:
+        // Fall through to return 0 measure.
+        break;
+    }
+    return OutType(0.0);
+#ifdef __NVCC__
+#pragma pop
+#endif
+#if defined(VTKM_MSVC)
+#pragma warning(pop)
+#endif
+  }
+};
+}
+} // namespace vtkm::worklet
+
+#endif // vtk_m_worklet_CellMeasure_h

vtkm/worklet/testing/CMakeLists.txt

@@ -25,6 +25,7 @@ set(unit_tests
   UnitTestCellGradient.cxx
   UnitTestCellSetConnectivity.cxx
   UnitTestCellSetDualGraph.cxx
+  UnitTestCellMeasure.cxx
   UnitTestClipping.cxx
   UnitTestContourTreeUniform.cxx
   UnitTestCosmoTools.cxx

vtkm/worklet/testing/UnitTestCellMeasure.cxx

@@ -0,0 +1,126 @@
+//============================================================================
+//  Copyright (c) Kitware, Inc.
+//  All rights reserved.
+//  See LICENSE.txt for details.
+//  This software is distributed WITHOUT ANY WARRANTY; without even
+//  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+//  PURPOSE.  See the above copyright notice for more information.
+//
+//  Copyright 2014 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
+//  Copyright 2014 UT-Battelle, LLC.
+//  Copyright 2014 Los Alamos National Security.
+//
+//  Under the terms of Contract DE-NA0003525 with NTESS,
+//  the U.S. Government retains certain rights in this software.
+//
+//  Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
+//  Laboratory (LANL), the U.S. Government retains certain rights in
+//  this software.
+//============================================================================
+
+#include "vtkm/worklet/CellMeasure.h"
+#include "vtkm/worklet/DispatcherMapTopology.h"
+
+#include "vtkm/cont/testing/MakeTestDataSet.h"
+#include "vtkm/cont/testing/Testing.h"
+
+namespace
+{
+
+void TestCellMeasureUniform3D()
+{
+  std::cout << "Testing CellMeasure Worklet on 3D structured data" << std::endl;
+
+  vtkm::cont::testing::MakeTestDataSet testDataSet;
+  vtkm::cont::DataSet dataSet = testDataSet.Make3DUniformDataSet0();
+
+  vtkm::cont::ArrayHandle<vtkm::Float32> result;
+
+  vtkm::worklet::DispatcherMapTopology<vtkm::worklet::CellMeasure<vtkm::Volume>> dispatcher;
+  dispatcher.Invoke(dataSet.GetCellSet(), dataSet.GetCoordinateSystem(), result);
+
+  vtkm::Float32 expected[4] = { 1.f, 1.f, 1.f, 1.f };
+  for (int i = 0; i < 4; ++i)
+  {
+    VTKM_TEST_ASSERT(test_equal(result.GetPortalConstControl().Get(vtkm::Id(i)), expected[i]),
+                     "Wrong result for CellMeasure worklet on 3D uniform data");
+  }
+}
+
+template <typename IntegrationType>
+void TestCellMeasureWorklet(vtkm::cont::DataSet& dataset,
+                            const char* msg,
+                            const std::vector<vtkm::Float32>& expected,
+                            const IntegrationType&)
+{
+  std::cout << "Testing CellMeasures Filter on " << msg << "\n";
+
+  vtkm::cont::ArrayHandle<vtkm::Float32> result;
+
+  vtkm::worklet::DispatcherMapTopology<vtkm::worklet::CellMeasure<IntegrationType>> dispatcher;
+  dispatcher.Invoke(dataset.GetCellSet(), dataset.GetCoordinateSystem(), result);
+
+  VTKM_TEST_ASSERT(result.GetNumberOfValues() == static_cast<vtkm::Id>(expected.size()),
+                   "Wrong number of values in the output array");
+
+  for (unsigned int i = 0; i < static_cast<unsigned int>(expected.size()); ++i)
+  {
+    VTKM_TEST_ASSERT(test_equal(result.GetPortalConstControl().Get(vtkm::Id(i)), expected[i]),
+                     "Wrong result for CellMeasure filter");
+  }
+}
+
+void TestCellMeasure()
+{
+  using vtkm::ArcLength;
+  using vtkm::Area;
+  using vtkm::Volume;
+  using vtkm::AllMeasures;
+
+  TestCellMeasureUniform3D();
+
+  vtkm::cont::testing::MakeTestDataSet factory;
+  vtkm::cont::DataSet data;
+
+  data = factory.Make3DExplicitDataSet2();
+  TestCellMeasureWorklet(data, "explicit dataset 2", { -1.f }, Volume());
+
+  data = factory.Make3DExplicitDataSet3();
+  TestCellMeasureWorklet(data, "explicit dataset 3", { -1.f / 6.f }, Volume());
+
+  data = factory.Make3DExplicitDataSet4();
+  TestCellMeasureWorklet(data, "explicit dataset 4", { -1.f, -1.f }, Volume());
+
+  data = factory.Make3DExplicitDataSet5();
+  TestCellMeasureWorklet(
+    data, "explicit dataset 5", { 1.f, 1.f / 3.f, 1.f / 6.f, -1.f / 2.f }, Volume());
+
+  data = factory.Make3DExplicitDataSet6();
+  TestCellMeasureWorklet(
+    data,
+    "explicit dataset 6 (all)",
+    { 0.999924f, 0.999924f, 0.f, 0.f, 3.85516f, 1.00119f, 0.083426f, 0.25028f },
+    AllMeasures());
+  TestCellMeasureWorklet(data,
+                         "explicit dataset 6 (arc length)",
+                         { 0.999924f, 0.999924f, 0.f, 0.f, 0.0f, 0.0f, 0.0f, 0.0f },
+                         ArcLength());
+  TestCellMeasureWorklet(data,
+                         "explicit dataset 6 (area)",
+                         { 0.0f, 0.0f, 0.f, 0.f, 3.85516f, 1.00119f, 0.0f, 0.0f },
+                         Area());
+  TestCellMeasureWorklet(data,
+                         "explicit dataset 6 (volume)",
+                         { 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.083426f, 0.25028f },
+                         Volume());
+  TestCellMeasureWorklet(data,
+                         "explicit dataset 6 (empty)",
+                         { 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.0f, 0.0f },
+                         vtkm::ListTagBase<>());
+}
+}
+
+int UnitTestCellMeasure(int, char* [])
+{
+  return vtkm::cont::testing::Testing::Run(TestCellMeasure);
+}