Merge branch 'upstream-lcl' into update-vtkc-to-lcl

* upstream-lcl:
  lcl 2019-10-07 (49094fec)

vtkm/thirdparty/lcl/vtkmlcl/LICENSE.md

@@ -0,0 +1,43 @@
+# Lightweight Cell Library License Version 1.0 #
+
+Copyright (c) 2019
+Kitware Inc.,
+National Technology & Engineering Solutions of Sandia, LLC (NTESS),
+UT-Battelle, LLC.,
+Los Alamos National Security, LLC.,
+All rights reserved.
+
+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.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+ * Redistributions of source code must retain the above copyright
+   notice, this list of conditions and the following disclaimer.
+
+ * Redistributions in binary form must reproduce the above copyright
+   notice, this list of conditions and the following disclaimer in the
+   documentation and/or other materials provided with the
+   distribution.
+
+ * Neither the name of Kitware nor the names of any contributors may
+   be used to endorse or promote products derived from this software
+   without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHORS OR
+CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

vtkm/thirdparty/lcl/vtkmlcl/README.md

@@ -0,0 +1,66 @@
+# Lightweight Cell Library #
+
+Lightweight Cell Library is a collection of cell types and cell functionality
+that is designed be used scientific visualization libraries of
+VTK-m, and VTK.
+
+You can find out more about the design of Lightweight Cell Library in [DESIGN.md].
+
+## Contributing ##
+
+There are many ways to contribute to [LCL]:
+
+  + Submit new or add to discussions of a feature requests or bugs on the
+    [LCL Issue Tracker].
+
+  + Submit a Pull Request
+      + See [CONTRIBUTING.md] for detailed instructions on how to create a
+        Pull Request.
+      + See the [LCL Coding Conventions] that must be followed for
+        contributed code.
+
+## Compiler Requirements ##
+
+  + C++11 Compiler. Lightweight Cell Library has been confirmed to work with
+    the following
+      + GCC 4.8+
+      + Clang 3.3+
+      + XCode 5.0+
+      + MSVC 2015+
+
+## Example ##
+
+Below is a simple example of using Lightweight Cell Library to get derivatives
+and parametric coordinates for different cell types:
+
+```cpp
+#incude <lcl/lcl.h>
+
+std::array<float, 3> pcoords;
+auto status = lcl::parametricCenter(lcl::Hexahedron{}, pcoords);
+
+std::vector<std::array<float, 3>> points = { ... };
+std::vector<std::array<double, 4>> field = { ... };
+std::array<double, 4> derivs[3];
+status = lcl::derivative(lcl::Hexahedron{},
+                          lcl::makeFieldAccessorNestedSOAConst(points, 3),
+                          lcl::makeFieldAccessorNestedSOAConst(field, 4),
+                          pcoords,
+                          derivs[0],
+                          derivs[1],
+                          derivs[2]);
+```
+
+## License ##
+
+Lightweight Cell Library is distributed under the OSI-approved BSD 3-clause
+License. See [LICENSE.md] for details.
+
+
+[LCL]:                                       https://gitlab.kitware.com/vtk/lcl/
+[LCL Issue Tracker]:                   https://gitlab.kitware.com/vtk/lcl/issues
+
+[CONTRIBUTING.md]:          CONTRIBUTING.md
+[DESIGN.md]:                docs/Design.md
+[LICENSE.md]:               LICENSE.md
+[LCL Coding Conventions]: docs/CodingConventions.md

vtkm/thirdparty/lcl/vtkmlcl/lcl/ErrorCode.h

@@ -0,0 +1,69 @@
+//============================================================================
+//  Copyright (c) Kitware, Inc.
+//  All rights reserved.
+//  See LICENSE.md 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.
+//============================================================================
+#ifndef lcl_ErrorCode_h
+#define lcl_ErrorCode_h
+
+#include <lcl/internal/Config.h>
+
+#include <cstdint>
+
+namespace lcl
+{
+
+enum class ErrorCode : std::int32_t
+{
+  SUCCESS = 0,
+  INVALID_SHAPE_ID,
+  INVALID_NUMBER_OF_POINTS,
+  WRONG_SHAPE_ID_FOR_TAG_TYPE,
+  INVALID_POINT_ID,
+  SOLUTION_DID_NOT_CONVERGE,
+  MATRIX_LUP_FACTORIZATION_FAILED,
+  DEGENERATE_CELL_DETECTED
+};
+
+LCL_EXEC
+inline const char* errorString(ErrorCode code) noexcept
+{
+  switch (code)
+  {
+    case ErrorCode::SUCCESS:
+      return "Success";
+    case ErrorCode::INVALID_SHAPE_ID:
+      return "Invalid shape id";
+    case ErrorCode::INVALID_NUMBER_OF_POINTS:
+      return "Invalid number of points";
+    case ErrorCode::WRONG_SHAPE_ID_FOR_TAG_TYPE:
+      return "Wrong shape id for tag type";
+    case ErrorCode::INVALID_POINT_ID:
+      return "Invalid point id";
+    case ErrorCode::SOLUTION_DID_NOT_CONVERGE:
+      return "Solution did not converge";
+    case ErrorCode::MATRIX_LUP_FACTORIZATION_FAILED:
+      return "LUP factorization failed";
+    case ErrorCode::DEGENERATE_CELL_DETECTED:
+      return "Degenerate cell detected";
+  }
+
+  return "Invalid error";
+}
+
+} // lcl
+
+#define LCL_RETURN_ON_ERROR(call)                                                                 \
+  {                                                                                                \
+    auto status = call;                                                                            \
+    if (status != lcl::ErrorCode::SUCCESS)                                                        \
+    {                                                                                              \
+      return status;                                                                               \
+    }                                                                                              \
+  }
+
+#endif // lcl_ErrorCode_h

vtkm/thirdparty/lcl/vtkmlcl/lcl/FieldAccessor.h

@@ -0,0 +1,367 @@
+//============================================================================
+//  Copyright (c) Kitware, Inc.
+//  All rights reserved.
+//  See LICENSE.md 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.
+//============================================================================
+#ifndef lcl_FieldAccessor_h
+#define lcl_FieldAccessor_h
+
+#include <lcl/internal/Config.h>
+
+#include <utility>
+
+namespace lcl
+{
+namespace internal
+{
+//----------------------------------------------------------------------------
+/// Does T support the indexing operator?
+/// Examples:
+///   std::array
+///   std::vector
+///   T*
+///   custom classes with operator[](integer_type)
+///
+
+struct NotVecTag {};
+
+template <typename T> inline decltype(std::declval<T>()[0]) VecCheck(T);
+inline NotVecTag VecCheck(...);
+
+template <typename T>
+using IsVecType =
+  std::integral_constant<bool, !std::is_same<decltype(VecCheck(std::declval<T>())),
+                                             internal::NotVecTag>::value>;
+
+//----------------------------------------------------------------------------
+/// Does T support operator()(integer_type)?
+
+struct NotFuncTag {};
+
+template <typename T> inline decltype(std::declval<T>()(0)) FuncCheck(T);
+inline NotFuncTag FuncCheck(...);
+
+template <typename T>
+using IsFuncType =
+  std::integral_constant<bool, !std::is_same<decltype(FuncCheck(std::declval<T>())),
+                                             internal::NotFuncTag>::value>;
+
+//----------------------------------------------------------------------------
+template <typename T> typename T::size_type IndexTypeImpl(T);
+IdComponent IndexTypeImpl(...);
+
+template <typename VecT>
+using IndexType = decltype(IndexTypeImpl(std::declval<VecT>()));
+
+//----------------------------------------------------------------------------
+struct VecTypeIndexer {};
+struct VecTypeFunctor {};
+struct VecTypeScalar  {};
+
+template <typename T>
+using VecTypeTag =
+  typename std::conditional<IsVecType<T>::value,
+                            VecTypeIndexer,
+                            typename std::conditional<IsFuncType<T>::value,
+                                                      VecTypeFunctor,
+                                                      VecTypeScalar>::type>::type;
+
+//----------------------------------------------------------------------------
+template <typename T>
+LCL_EXEC auto componentImpl(VecTypeIndexer, T&& vec, int idx) -> decltype(vec[0])
+{
+  return vec[static_cast<IndexType<T>>(idx)];
+}
+
+template <typename T>
+LCL_EXEC auto componentImpl(VecTypeFunctor, T&& vec, int idx) -> decltype(vec(0))
+{
+  return vec(idx);
+}
+
+template <typename T>
+LCL_EXEC T& componentImpl(VecTypeScalar, T&& vec, int)
+{
+  return vec;
+}
+
+} // namespace internal
+
+//----------------------------------------------------------------------------
+template <typename T>
+LCL_EXEC auto component(T&& vec, int idx)
+  -> decltype(internal::componentImpl<T>(internal::VecTypeTag<T>{}, std::forward<T>(vec), 0))
+{
+  return internal::componentImpl<T>(internal::VecTypeTag<T>{}, std::forward<T>(vec), idx);
+}
+
+//----------------------------------------------------------------------------
+template <typename T>
+using ComponentType = typename std::decay<decltype(component(std::declval<T>(), 0))>::type;
+
+///============================================================================
+/// Since there are different ways fields maybe represented in the clients of
+/// this libarary, LCL relies on helper classes that implement the
+/// \c FieldAccessor "concept" to access the elements of a field.
+///
+/// These classes should wrap the field and provide the follolwing interface:
+///
+/// template <typename FieldType>
+/// class FieldAccessor
+/// {
+/// public:
+///   /// An alias for the component type of the field
+///   using ValueType = ...;
+///
+///   /// Return the number of components
+///   int getNumberOfComponents() const;
+///
+///   /// Set the value to `tuple` and `comp` to `value`.
+///   void setValue(int tuple, int comp, const ValueType& value) const;
+///
+///   /// Get the value at `tuple` and `comp`.
+///   ValueType getValue(int tuple, int comp) const;
+///
+///   /// Set the tuple at index `tuple`. It is recomended to make this a
+///   /// template function and use `lcl::component` to access
+///   /// the components of `value`.
+///   template <typename VecType>
+///   void setTuple(int tuple, const VecType& value) const;
+///
+///   /// Get the tuple at index `tuple`. It is recomended to make this a
+///   /// template function and use `lcl::component` to access
+///   /// the components of `value`.
+///   template <typename VecType>
+///   void getTuple(int tuple, VecType& value) const;
+/// }
+///
+/// The set functions are optional and such a class would act as a const variant
+/// of \c FieldAccessor.
+///
+
+///----------------------------------------------------------------------------
+template <typename FieldType>
+class FieldAccessorNestedSOA
+{
+public:
+  using ValueType = typename std::decay<decltype(component(std::declval<FieldType>()[0], 0))>::type;
+
+  LCL_EXEC FieldAccessorNestedSOA(FieldType& field, int numberOfComponents = 1)
+    : Field(&field), NumberOfComponents(numberOfComponents)
+  {
+  }
+
+  LCL_EXEC int getNumberOfComponents() const
+  {
+    return this->NumberOfComponents;
+  }
+
+  LCL_EXEC void setValue(int tuple, int comp, const ValueType& value) const
+  {
+    component((*this->Field)[static_cast<IdxType>(tuple)], comp) = value;
+  }
+
+  LCL_EXEC ValueType getValue(int tuple, int comp) const
+  {
+    return component((*this->Field)[static_cast<IdxType>(tuple)], comp);
+  }
+
+  template <typename VecType>
+  LCL_EXEC void setTuple(int tuple, const VecType& value) const
+  {
+    for (int i = 0; i < this->NumberOfComponents; ++i)
+    {
+      component((*this->Field)[static_cast<IdxType>(tuple)], i) =
+        static_cast<ValueType>(component(value, i));
+    }
+  }
+
+  template <typename VecType>
+  LCL_EXEC void getTuple(int tuple, VecType& value) const
+  {
+    for (int i = 0; i < this->NumberOfComponents; ++i)
+    {
+      component(value, i) = static_cast<ComponentType<VecType>>(
+        component((*this->Field)[static_cast<IdxType>(tuple)], i));
+    }
+  }
+
+private:
+  using IdxType = internal::IndexType<FieldType>;
+
+  FieldType* Field;
+  int NumberOfComponents;
+};
+
+template <typename FieldType>
+class FieldAccessorNestedSOAConst
+{
+public:
+  using ValueType = typename std::decay<decltype(component(std::declval<FieldType>()[0], 0))>::type;
+
+  LCL_EXEC FieldAccessorNestedSOAConst(const FieldType& field, int numberOfComponents = 1)
+    : Field(&field), NumberOfComponents(numberOfComponents)
+  {
+  }
+
+  LCL_EXEC int getNumberOfComponents() const
+  {
+    return this->NumberOfComponents;
+  }
+
+  LCL_EXEC ValueType getValue(int tuple, int comp) const
+  {
+    return component((*this->Field)[static_cast<IdxType>(tuple)], comp);
+  }
+
+  template <typename VecType>
+  LCL_EXEC void getTuple(int tuple, VecType& value) const
+  {
+    for (int i = 0; i < this->NumberOfComponents; ++i)
+    {
+      component(value, i) = static_cast<ComponentType<VecType>>(
+        component((*this->Field)[static_cast<IdxType>(tuple)], i));
+    }
+  }
+
+private:
+  using IdxType = internal::IndexType<FieldType>;
+
+  const FieldType* Field;
+  int NumberOfComponents;
+};
+
+template <typename FieldType>
+LCL_EXEC
+FieldAccessorNestedSOA<FieldType> makeFieldAccessorNestedSOA(FieldType& field,
+                                                             int numberOfComponents = 1)
+{
+  return FieldAccessorNestedSOA<FieldType>(field, numberOfComponents);
+}
+
+template <typename FieldType>
+LCL_EXEC
+FieldAccessorNestedSOAConst<FieldType> makeFieldAccessorNestedSOAConst(const FieldType& field,
+                                                                       int numberOfComponents = 1)
+{
+  return FieldAccessorNestedSOAConst<FieldType>(field, numberOfComponents);
+}
+
+///----------------------------------------------------------------------------
+template <typename FieldType>
+class FieldAccessorFlatSOA
+{
+public:
+  using ValueType = typename std::decay<decltype(std::declval<FieldType>()[0])>::type;
+
+  LCL_EXEC FieldAccessorFlatSOA(FieldType& field, int numberOfComponents = 1)
+    : Field(&field), NumberOfComponents(numberOfComponents)
+  {
+  }
+
+  LCL_EXEC int getNumberOfComponents() const
+  {
+    return this->NumberOfComponents;
+  }
+
+  LCL_EXEC void setValue(int tuple, int comp, const ValueType& value) const
+  {
+    auto FlatIdx = static_cast<IdxType>(tuple * this->NumberOfComponents + comp);
+    (*this->Field)[FlatIdx] = value;
+  }
+
+  LCL_EXEC ValueType getValue(int tuple, int comp) const
+  {
+    auto FlatIdx = static_cast<IdxType>(tuple * this->NumberOfComponents + comp);
+    return (*this->Field)[FlatIdx];
+  }
+
+  template <typename VecType>
+  LCL_EXEC void setTuple(int tuple, const VecType& value) const
+  {
+    auto start = static_cast<IdxType>(tuple * this->NumberOfComponents);
+    for (int i = 0; i < this->NumberOfComponents; ++i)
+    {
+      (*this->Field)[start++] = static_cast<ValueType>(component(value, i));
+    }
+  }
+
+  template <typename VecType>
+  LCL_EXEC void getTuple(int tuple, VecType& value) const
+  {
+    auto start = static_cast<IdxType>(tuple * this->NumberOfComponents);
+    for (int i = 0; i < this->NumberOfComponents; ++i)
+    {
+      component(value, i) = static_cast<ComponentType<VecType>>((*this->Field)[start++]);
+    }
+  }
+
+private:
+  using IdxType = internal::IndexType<FieldType>;
+
+  FieldType* Field;
+  int NumberOfComponents;
+};
+
+template <typename FieldType>
+class FieldAccessorFlatSOAConst
+{
+public:
+  using ValueType = typename std::decay<decltype(std::declval<FieldType>()[0])>::type;
+
+  LCL_EXEC FieldAccessorFlatSOAConst(const FieldType& field, int numberOfComponents = 1)
+    : Field(&field), NumberOfComponents(numberOfComponents)
+  {
+  }
+
+  LCL_EXEC int getNumberOfComponents() const
+  {
+    return this->NumberOfComponents;
+  }
+
+  LCL_EXEC ValueType getValue(int tuple, int comp) const
+  {
+    auto FlatIdx = static_cast<IdxType>(tuple * this->NumberOfComponents + comp);
+    return (*this->Field)[FlatIdx];
+  }
+
+  template <typename VecType>
+  LCL_EXEC void getTuple(int tuple, VecType& value) const
+  {
+    auto start = static_cast<IdxType>(tuple * this->NumberOfComponents);
+    for (int i = 0; i < this->NumberOfComponents; ++i)
+    {
+      component(value, i) = static_cast<ComponentType<VecType>>((*this->Field)[start++]);
+    }
+  }
+
+private:
+  using IdxType = internal::IndexType<FieldType>;
+
+  const FieldType* Field;
+  int NumberOfComponents;
+};
+
+template <typename FieldType>
+LCL_EXEC
+FieldAccessorFlatSOA<FieldType> makeFieldAccessorFlatSOA(FieldType& field,
+                                                         int numberOfComponents)
+{
+  return FieldAccessorFlatSOA<FieldType>(field, numberOfComponents);
+}
+
+template <typename FieldType>
+LCL_EXEC
+FieldAccessorFlatSOAConst<FieldType> makeFieldAccessorFlatSOAConst(const FieldType& field,
+                                                                   int numberOfComponents)
+{
+  return FieldAccessorFlatSOAConst<FieldType>(field, numberOfComponents);
+}
+
+} // namespace lcl
+
+#endif // lcl_FieldAccessor_h

vtkm/thirdparty/lcl/vtkmlcl/lcl/Hexahedron.h

@@ -0,0 +1,256 @@
+//============================================================================
+//  Copyright (c) Kitware, Inc.
+//  All rights reserved.
+//  See LICENSE.md 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.
+//============================================================================
+#ifndef lcl_Hexahedron_h
+#define lcl_Hexahedron_h
+
+#include <lcl/ErrorCode.h>
+#include <lcl/Shapes.h>
+
+#include <lcl/internal/Common.h>
+
+namespace lcl
+{
+
+class Hexahedron : public Cell
+{
+public:
+  constexpr LCL_EXEC Hexahedron() : Cell(ShapeId::HEXAHEDRON, 8) {}
+  constexpr LCL_EXEC explicit Hexahedron(const Cell& cell) : Cell(cell) {}
+};
+
+LCL_EXEC inline lcl::ErrorCode validate(Hexahedron tag) noexcept
+{
+  if (tag.shape() != ShapeId::HEXAHEDRON && tag.shape() != ShapeId::VOXEL)
+  {
+    return ErrorCode::WRONG_SHAPE_ID_FOR_TAG_TYPE;
+  }
+  if (tag.numberOfPoints() != 8)
+  {
+    return ErrorCode::INVALID_NUMBER_OF_POINTS;
+  }
+
+  return ErrorCode::SUCCESS;
+}
+
+template<typename CoordType>
+LCL_EXEC inline lcl::ErrorCode parametricCenter(Hexahedron, CoordType&& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  component(pcoords, 0) = 0.5f;
+  component(pcoords, 1) = 0.5f;
+  component(pcoords, 2) = 0.5f;
+  return ErrorCode::SUCCESS;
+}
+
+template<typename CoordType>
+LCL_EXEC inline lcl::ErrorCode parametricPoint(
+  Hexahedron, IdComponent pointId, CoordType&& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  switch (pointId)
+  {
+    case 0:
+      component(pcoords, 0) = 0.0f;
+      component(pcoords, 1) = 0.0f;
+      component(pcoords, 2) = 0.0f;
+      break;
+    case 1:
+      component(pcoords, 0) = 1.0f;
+      component(pcoords, 1) = 0.0f;
+      component(pcoords, 2) = 0.0f;
+      break;
+    case 2:
+      component(pcoords, 0) = 1.0f;
+      component(pcoords, 1) = 1.0f;
+      component(pcoords, 2) = 0.0f;
+      break;
+    case 3:
+      component(pcoords, 0) = 0.0f;
+      component(pcoords, 1) = 1.0f;
+      component(pcoords, 2) = 0.0f;
+      break;
+    case 4:
+      component(pcoords, 0) = 0.0f;
+      component(pcoords, 1) = 0.0f;
+      component(pcoords, 2) = 1.0f;
+      break;
+    case 5:
+      component(pcoords, 0) = 1.0f;
+      component(pcoords, 1) = 0.0f;
+      component(pcoords, 2) = 1.0f;
+      break;
+    case 6:
+      component(pcoords, 0) = 1.0f;
+      component(pcoords, 1) = 1.0f;
+      component(pcoords, 2) = 1.0f;
+      break;
+    case 7:
+      component(pcoords, 0) = 0.0f;
+      component(pcoords, 1) = 1.0f;
+      component(pcoords, 2) = 1.0f;
+      break;
+    default:
+      return ErrorCode::INVALID_POINT_ID;
+  }
+
+  return ErrorCode::SUCCESS;
+}
+
+template<typename CoordType>
+LCL_EXEC inline ComponentType<CoordType> parametricDistance(Hexahedron, const CoordType& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+  return internal::findParametricDistance(pcoords, 3);
+}
+
+template<typename CoordType>
+LCL_EXEC inline bool cellInside(Hexahedron, const CoordType& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  using T = ComponentType<CoordType>;
+
+  constexpr T eps = 1e-6f;
+  return component(pcoords, 0) >= -eps && component(pcoords, 0) <= (T{1} + eps) &&
+         component(pcoords, 1) >= -eps && component(pcoords, 1) <= (T{1} + eps) &&
+         component(pcoords, 2) >= -eps && component(pcoords, 2) <= (T{1} + eps);
+}
+
+template <typename Values, typename CoordType, typename Result>
+LCL_EXEC inline lcl::ErrorCode interpolate(
+  Hexahedron,
+  const Values& values,
+  const CoordType& pcoords,
+  Result&& result) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  using T = internal::ClosestFloatType<typename Values::ValueType>;
+
+  for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
+  {
+    auto vbf = internal::lerp(static_cast<T>(values.getValue(0, c)),
+                              static_cast<T>(values.getValue(1, c)),
+                              static_cast<T>(component(pcoords, 0)));
+    auto vbb = internal::lerp(static_cast<T>(values.getValue(3, c)),
+                              static_cast<T>(values.getValue(2, c)),
+                              static_cast<T>(component(pcoords, 0)));
+    auto vtf = internal::lerp(static_cast<T>(values.getValue(4, c)),
+                              static_cast<T>(values.getValue(5, c)),
+                              static_cast<T>(component(pcoords, 0)));
+    auto vtb = internal::lerp(static_cast<T>(values.getValue(7, c)),
+                              static_cast<T>(values.getValue(6, c)),
+                              static_cast<T>(component(pcoords, 0)));
+    auto vb = internal::lerp(vbf, vbb, static_cast<T>(component(pcoords, 1)));
+    auto vt = internal::lerp(vtf, vtb, static_cast<T>(component(pcoords, 1)));
+    auto v = internal::lerp(vb, vt, static_cast<T>(component(pcoords, 2)));
+    component(result, c) = static_cast<ComponentType<Result>>(v);
+  }
+
+  return ErrorCode::SUCCESS;
+}
+
+namespace internal
+{
+
+template <typename Values, typename CoordType, typename Result>
+LCL_EXEC inline void parametricDerivative(Hexahedron,
+                                           const Values& values,
+                                           IdComponent comp,
+                                           const CoordType& pcoords,
+                                           Result&& result) noexcept
+{
+  using T = internal::ClosestFloatType<typename Values::ValueType>;
+  T p0 = static_cast<T>(component(pcoords, 0));
+  T p1 = static_cast<T>(component(pcoords, 1));
+  T p2 = static_cast<T>(component(pcoords, 2));
+  T rm = T{1} - p0;
+  T sm = T{1} - p1;
+  T tm = T{1} - p2;
+
+  T dr = (static_cast<T>(values.getValue(0, comp)) * -sm * tm) +
+         (static_cast<T>(values.getValue(1, comp)) *  sm * tm) +
+         (static_cast<T>(values.getValue(2, comp)) *  p1 * tm) +
+         (static_cast<T>(values.getValue(3, comp)) * -p1 * tm) +
+         (static_cast<T>(values.getValue(4, comp)) * -sm * p2) +
+         (static_cast<T>(values.getValue(5, comp)) *  sm * p2) +
+         (static_cast<T>(values.getValue(6, comp)) *  p1 * p2) +
+         (static_cast<T>(values.getValue(7, comp)) * -p1 * p2);
+
+  T ds = (static_cast<T>(values.getValue(0, comp)) * -rm * tm) +
+         (static_cast<T>(values.getValue(1, comp)) * -p0 * tm) +
+         (static_cast<T>(values.getValue(2, comp)) *  p0 * tm) +
+         (static_cast<T>(values.getValue(3, comp)) *  rm * tm) +
+         (static_cast<T>(values.getValue(4, comp)) * -rm * p2) +
+         (static_cast<T>(values.getValue(5, comp)) * -p0 * p2) +
+         (static_cast<T>(values.getValue(6, comp)) *  p0 * p2) +
+         (static_cast<T>(values.getValue(7, comp)) *  rm * p2);
+
+  T dt = (static_cast<T>(values.getValue(0, comp)) * -rm * sm) +
+         (static_cast<T>(values.getValue(1, comp)) * -p0 * sm) +
+         (static_cast<T>(values.getValue(2, comp)) * -p0 * p1) +
+         (static_cast<T>(values.getValue(3, comp)) * -rm * p1) +
+         (static_cast<T>(values.getValue(4, comp)) *  rm * sm) +
+         (static_cast<T>(values.getValue(5, comp)) *  p0 * sm) +
+         (static_cast<T>(values.getValue(6, comp)) *  p0 * p1) +
+         (static_cast<T>(values.getValue(7, comp)) *  rm * p1);
+
+  component(result, 0) = static_cast<ComponentType<Result>>(dr);
+  component(result, 1) = static_cast<ComponentType<Result>>(ds);
+  component(result, 2) = static_cast<ComponentType<Result>>(dt);
+}
+
+} // internal
+
+template <typename Points, typename Values, typename CoordType, typename Result>
+LCL_EXEC inline lcl::ErrorCode derivative(
+  Hexahedron,
+  const Points& points,
+  const Values& values,
+  const CoordType& pcoords,
+  Result&& dx,
+  Result&& dy,
+  Result&& dz) noexcept
+{
+  return internal::derivative3D(Hexahedron{},
+                                points,
+                                values,
+                                pcoords,
+                                std::forward<Result>(dx),
+                                std::forward<Result>(dy),
+                                std::forward<Result>(dz));
+}
+
+template <typename Points, typename PCoordType, typename WCoordType>
+LCL_EXEC inline lcl::ErrorCode parametricToWorld(
+  Hexahedron,
+  const Points& points,
+  const PCoordType& pcoords,
+  WCoordType&& wcoords) noexcept
+{
+  return interpolate(Hexahedron{}, points, pcoords, std::forward<WCoordType>(wcoords));
+}
+
+template <typename Points, typename WCoordType, typename PCoordType>
+LCL_EXEC inline lcl::ErrorCode worldToParametric(
+  Hexahedron,
+  const Points& points,
+  const WCoordType& wcoords,
+  PCoordType&& pcoords) noexcept
+{
+  return internal::worldToParametric3D(
+    Hexahedron{}, points, wcoords, std::forward<PCoordType>(pcoords));
+}
+
+} // lcl
+
+#endif // lcl_Hexahedron_h

vtkm/thirdparty/lcl/vtkmlcl/lcl/Line.h

@@ -0,0 +1,176 @@
+//============================================================================
+//  Copyright (c) Kitware, Inc.
+//  All rights reserved.
+//  See LICENSE.md 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.
+//============================================================================
+#ifndef lcl_Line_h
+#define lcl_Line_h
+
+#include <lcl/ErrorCode.h>
+#include <lcl/Shapes.h>
+
+#include <lcl/internal/Common.h>
+
+namespace lcl
+{
+
+class Line : public Cell
+{
+public:
+  constexpr LCL_EXEC Line() : Cell(ShapeId::LINE, 2) {}
+  constexpr LCL_EXEC explicit Line(const Cell& cell) : Cell(cell) {}
+};
+
+LCL_EXEC inline lcl::ErrorCode validate(Line tag) noexcept
+{
+  if (tag.shape() != ShapeId::LINE)
+  {
+    return ErrorCode::WRONG_SHAPE_ID_FOR_TAG_TYPE;
+  }
+  if (tag.numberOfPoints() != 2)
+  {
+    return ErrorCode::INVALID_NUMBER_OF_POINTS;
+  }
+
+  return ErrorCode::SUCCESS;
+}
+
+template<typename CoordType>
+LCL_EXEC inline lcl::ErrorCode parametricCenter(Line, CoordType&& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  component(pcoords, 0) = 0.5f;
+  return ErrorCode::SUCCESS;
+}
+
+template<typename CoordType>
+LCL_EXEC inline lcl::ErrorCode parametricPoint(
+  Line, IdComponent pointId, CoordType&& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  switch (pointId)
+  {
+    case 0:
+      component(pcoords, 0) = 0.0f;
+      return ErrorCode::SUCCESS;
+    case 1:
+      component(pcoords, 0) = 1.0f;
+      return ErrorCode::SUCCESS;
+    default:
+      return ErrorCode::INVALID_POINT_ID;
+  }
+}
+
+template<typename CoordType>
+LCL_EXEC inline ComponentType<CoordType> parametricDistance(Line, const CoordType& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+  return internal::findParametricDistance(pcoords, 1);
+}
+
+template<typename CoordType>
+LCL_EXEC inline bool cellInside(Line, const CoordType& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  using T = ComponentType<CoordType>;
+  return component(pcoords, 0) >= T{0} && component(pcoords, 0) <= T{1};
+}
+
+template <typename Values, typename CoordType, typename Result>
+LCL_EXEC inline lcl::ErrorCode interpolate(
+  Line,
+  const Values& values,
+  const CoordType& pcoords,
+  Result&& result) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  using ProcessingType = internal::ClosestFloatType<typename Values::ValueType>;
+  using ResultCompType = ComponentType<Result>;
+
+  for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
+  {
+    auto ival = internal::lerp(static_cast<ProcessingType>(values.getValue(0, c)),
+                               static_cast<ProcessingType>(values.getValue(1, c)),
+                               static_cast<ProcessingType>(component(pcoords, 0)));
+    component(result, c) = static_cast<ResultCompType>(ival);
+  }
+
+  return ErrorCode::SUCCESS;
+}
+
+template <typename Points, typename Values, typename CoordType, typename Result>
+LCL_EXEC inline lcl::ErrorCode derivative(
+  Line,
+  const Points& points,
+  const Values& values,
+  const CoordType&,
+  Result&& dx,
+  Result&& dy,
+  Result&& dz) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  using ProcessingType = internal::ClosestFloatType<typename Values::ValueType>;
+  using ResultCompType = ComponentType<Result>;
+
+  ProcessingType dPt[3] = {
+    static_cast<ProcessingType>(points.getValue(1, 0) - points.getValue(0, 0)),
+    static_cast<ProcessingType>(points.getValue(1, 1) - points.getValue(0, 1)),
+    static_cast<ProcessingType>(points.getValue(1, 2) - points.getValue(0, 2)) };
+
+  for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
+  {
+    auto dv = static_cast<ProcessingType>(values.getValue(1, c) - values.getValue(0, c));
+    component(dx, c) = (dPt[0] != 0.0f) ? static_cast<ResultCompType>(dv/dPt[0]) : ResultCompType{0};
+    component(dy, c) = (dPt[1] != 0.0f) ? static_cast<ResultCompType>(dv/dPt[1]) : ResultCompType{0};
+    component(dz, c) = (dPt[2] != 0.0f) ? static_cast<ResultCompType>(dv/dPt[2]) : ResultCompType{0};
+  }
+  return ErrorCode::SUCCESS;
+}
+
+template <typename Points, typename PCoordType, typename WCoordType>
+LCL_EXEC inline lcl::ErrorCode parametricToWorld(
+  Line,
+  const Points& points,
+  const PCoordType& pcoords,
+  WCoordType&& wcoords) noexcept
+{
+  return interpolate(Line{}, points, pcoords, std::forward<WCoordType>(wcoords));
+}
+
+template <typename Points, typename WCoordType, typename PCoordType>
+LCL_EXEC inline lcl::ErrorCode worldToParametric(
+  Line,
+  const Points& points,
+  const WCoordType& wcoords,
+  PCoordType&& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(PCoordType);
+
+  using T = ComponentType<PCoordType>;
+  internal::Vector<T, 3> p0(static_cast<T>(points.getValue(0, 0)),
+                            static_cast<T>(points.getValue(0, 1)),
+                            static_cast<T>(points.getValue(0, 2)));
+  internal::Vector<T, 3> p1(static_cast<T>(points.getValue(1, 0)),
+                            static_cast<T>(points.getValue(1, 1)),
+                            static_cast<T>(points.getValue(1, 2)));
+  internal::Vector<T, 3> wc(static_cast<T>(component(wcoords, 0)),
+                            static_cast<T>(component(wcoords, 1)),
+                            static_cast<T>(component(wcoords, 2)));
+  auto v1 = p1 - p0;
+  auto v2 = wc - p0;
+  component(pcoords, 0) = internal::dot(v1, v2) / internal::dot(v1, v1);
+  return ErrorCode::SUCCESS;
+}
+
+} //namespace lcl
+
+#endif //lcl_Line_h

vtkm/thirdparty/lcl/vtkmlcl/lcl/Pixel.h

@@ -0,0 +1,178 @@
+//============================================================================
+//  Copyright (c) Kitware, Inc.
+//  All rights reserved.
+//  See LICENSE.md 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.
+//============================================================================
+#ifndef lcl_Pixel_h
+#define lcl_Pixel_h
+
+#include <lcl/ErrorCode.h>
+#include <lcl/Quad.h>
+
+#include <lcl/internal/Common.h>
+
+namespace lcl
+{
+
+class Pixel : public Quad
+{
+public:
+  constexpr LCL_EXEC Pixel() : Quad(Cell(ShapeId::PIXEL, 4)) {}
+  constexpr LCL_EXEC explicit Pixel(const Cell& cell) : Quad(cell) {}
+};
+
+namespace internal
+{
+
+template <typename Points, typename T>
+LCL_EXEC inline int getPixelSpacing(const Points& points, T spacing[3])
+{
+  int zeros = 0;
+  for (int i = 0; i < 3; ++i)
+  {
+    spacing[i] = static_cast<T>(points.getValue(2, i) - points.getValue(0, i));
+    if (spacing[i] == T{0})
+    {
+      zeros |= 1 << i;
+    }
+  }
+  return zeros;
+}
+
+} // internal
+
+template <typename Points, typename Values, typename CoordType, typename Result>
+LCL_EXEC inline lcl::ErrorCode derivative(
+  Pixel,
+  const Points& points,
+  const Values& values,
+  const CoordType& pcoords,
+  Result&& dx,
+  Result&& dy,
+  Result&& dz) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  using ProcessingType = internal::ClosestFloatType<typename Values::ValueType>;
+  using ResultCompType = ComponentType<Result>;
+
+  ProcessingType spacing[3];
+  int zeros = internal::getPixelSpacing(points, spacing);
+
+  for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
+  {
+    ProcessingType dvdp[2];
+    internal::parametricDerivative(Quad{}, values, c, pcoords, dvdp);
+    switch (zeros)
+    {
+      case 1: // yz plane
+        component(dx, c) = ResultCompType{0};
+        component(dy, c) = static_cast<ResultCompType>(dvdp[0] / spacing[1]);
+        component(dz, c) = static_cast<ResultCompType>(dvdp[1] / spacing[2]);
+        break;
+      case 2: // xz plane
+        component(dx, c) = static_cast<ResultCompType>(dvdp[0] / spacing[0]);
+        component(dy, c) = ResultCompType{0};
+        component(dz, c) = static_cast<ResultCompType>(dvdp[1] / spacing[2]);
+        break;
+      case 4: // xy plane
+        component(dx, c) = static_cast<ResultCompType>(dvdp[0] / spacing[0]);
+        component(dy, c) = static_cast<ResultCompType>(dvdp[1] / spacing[1]);
+        component(dz, c) = ResultCompType{0};
+        break;
+      default:
+        return ErrorCode::DEGENERATE_CELL_DETECTED;
+    }
+  }
+
+  return ErrorCode::SUCCESS;
+}
+
+template <typename Points, typename PCoordType, typename WCoordType>
+LCL_EXEC lcl::ErrorCode parametricToWorld(
+  Pixel,
+  const Points& points,
+  const PCoordType& pcoords,
+  WCoordType&& wcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(PCoordType);
+
+  using T = typename Points::ValueType;
+
+  T spacing[3];
+  auto zeros = internal::getPixelSpacing(points, spacing);
+
+  switch (zeros)
+  {
+    case 1: // yz plane
+      component(wcoords, 0) = points.getValue(0, 0);
+      component(wcoords, 1) = points.getValue(0, 1) +
+                              (spacing[1] * static_cast<T>(component(pcoords, 0)));
+      component(wcoords, 2) = points.getValue(0, 2) +
+                              (spacing[2] * static_cast<T>(component(pcoords, 1)));
+      return ErrorCode::SUCCESS;
+    case 2: // xz plane
+      component(wcoords, 0) = points.getValue(0, 0) +
+                              (spacing[0] * static_cast<T>(component(pcoords, 0)));
+      component(wcoords, 1) = points.getValue(0, 1);
+      component(wcoords, 2) = points.getValue(0, 2) +
+                              (spacing[2] * static_cast<T>(component(pcoords, 1)));
+      return ErrorCode::SUCCESS;
+    case 4: // xy plane
+      component(wcoords, 0) = points.getValue(0, 0) +
+                              (spacing[0] * static_cast<T>(component(pcoords, 0)));
+      component(wcoords, 1) = points.getValue(0, 1) +
+                              (spacing[1] * static_cast<T>(component(pcoords, 1)));
+      component(wcoords, 2) = points.getValue(0, 2);
+      return ErrorCode::SUCCESS;
+    default:
+      return ErrorCode::DEGENERATE_CELL_DETECTED;
+  }
+}
+
+template <typename Points, typename WCoordType, typename PCoordType>
+LCL_EXEC lcl::ErrorCode worldToParametric(
+  Pixel,
+  const Points& points,
+  const WCoordType& wcoords,
+  PCoordType&& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(PCoordType);
+
+  using T = ComponentType<PCoordType>;
+
+  T spacing[3];
+  int zeros = internal::getPixelSpacing(points, spacing);
+
+  switch (zeros)
+  {
+  case 1: // yz plane
+    component(pcoords, 0) = static_cast<T>(component(wcoords, 1) - points.getValue(0, 1)) /
+                            spacing[1];
+    component(pcoords, 1) = static_cast<T>(component(wcoords, 2) - points.getValue(0, 2)) /
+                            spacing[2];
+    return ErrorCode::SUCCESS;
+  case 2: // xz plane
+    component(pcoords, 0) = static_cast<T>(component(wcoords, 0) - points.getValue(0, 0)) /
+                            spacing[0];
+    component(pcoords, 1) = static_cast<T>(component(wcoords, 2) - points.getValue(0, 2)) /
+                            spacing[2];
+    return ErrorCode::SUCCESS;
+  case 4: // xy plane
+    component(pcoords, 0) = static_cast<T>(component(wcoords, 0) - points.getValue(0, 0)) /
+                            spacing[0];
+    component(pcoords, 1) = static_cast<T>(component(wcoords, 1) - points.getValue(0, 1)) /
+                            spacing[1];
+    return ErrorCode::SUCCESS;
+  default:
+    return ErrorCode::DEGENERATE_CELL_DETECTED;
+  }
+}
+
+} //namespace lcl
+
+#endif //lcl_Pixel_h

vtkm/thirdparty/lcl/vtkmlcl/lcl/Polygon.h

@@ -0,0 +1,576 @@
+//============================================================================
+//  Copyright (c) Kitware, Inc.
+//  All rights reserved.
+//  See LICENSE.md 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.
+//============================================================================
+#ifndef lcl_Polygon_h
+#define lcl_Polygon_h
+
+#include <lcl/ErrorCode.h>
+#include <lcl/Quad.h>
+#include <lcl/Shapes.h>
+#include <lcl/Triangle.h>
+
+#include <lcl/internal/Common.h>
+
+namespace lcl
+{
+
+/// \c Polygon with 3 and 4 points behave exactly as \c Triangle and \c Quad
+/// respectively. For 5 or more points, the points are arranged such that
+/// they are on the circle circumscribed in the
+/// unit square from 0 to 1. That is, the point are on the circle centered at
+/// coordinate 0.5,0.5 with radius 0.5. The polygon is divided into regions
+/// defined by the triangle fan formed by the points around the center. This
+/// is C0 continuous but not necessarily C1 continuous. It is also possible to
+/// have a non 1 to 1 mapping between parametric coordinates world coordinates
+/// if the polygon is not planar or convex.
+class Polygon : public Cell
+{
+public:
+  constexpr LCL_EXEC Polygon() : Cell(ShapeId::POLYGON, 3) {}
+  constexpr LCL_EXEC explicit Polygon(lcl::IdComponent numPoints)
+    : Cell(ShapeId::POLYGON, numPoints)
+  {
+  }
+  constexpr LCL_EXEC explicit Polygon(const Cell& cell) : Cell(cell) {}
+};
+
+LCL_EXEC inline lcl::ErrorCode validate(Polygon tag) noexcept
+{
+  if (tag.shape() != ShapeId::POLYGON)
+  {
+    return ErrorCode::WRONG_SHAPE_ID_FOR_TAG_TYPE;
+  }
+  if (tag.numberOfPoints() < 3)
+  {
+    return ErrorCode::INVALID_NUMBER_OF_POINTS;
+  }
+
+  return ErrorCode::SUCCESS;
+}
+
+template<typename CoordType>
+LCL_EXEC inline lcl::ErrorCode parametricCenter(Polygon tag, CoordType&& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  switch (tag.numberOfPoints())
+  {
+    case 3:
+      return parametricCenter(Triangle{}, pcoords);
+    case 4:
+      return parametricCenter(Quad{}, pcoords);
+    default:
+      component(pcoords, 0) = 0.5f;
+      component(pcoords, 1) = 0.5f;
+      return ErrorCode::SUCCESS;
+  }
+}
+
+template<typename CoordType>
+LCL_EXEC inline lcl::ErrorCode parametricPoint(
+  Polygon tag, IdComponent pointId, CoordType&& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  if (pointId < 0 || pointId >= tag.numberOfPoints())
+  {
+    return ErrorCode::INVALID_POINT_ID;
+  }
+
+  switch (tag.numberOfPoints())
+  {
+    case 3:
+      return parametricPoint(Triangle{}, pointId, pcoords);
+    case 4:
+      return parametricPoint(Quad{}, pointId, pcoords);
+    default:
+    {
+      using T = ComponentType<CoordType>;
+      constexpr double two_pi = 2.0 * 3.14159265359;
+      auto angle = (static_cast<T>(pointId) * static_cast<T>(two_pi)) / static_cast<T>(tag.numberOfPoints());
+      component(pcoords, 0) = 0.5f * (LCL_MATH_CALL(cos, (angle)) + 1.0f);
+      component(pcoords, 1) = 0.5f * (LCL_MATH_CALL(sin, (angle)) + 1.0f);
+      return ErrorCode::SUCCESS;
+    }
+  }
+}
+
+template<typename CoordType>
+LCL_EXEC inline ComponentType<CoordType> parametricDistance(Polygon tag, const CoordType& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  switch (tag.numberOfPoints())
+  {
+    case 3:
+      return parametricDistance(Triangle{}, pcoords);
+    default:
+      return internal::findParametricDistance(pcoords, 2);
+  }
+}
+
+template<typename CoordType>
+LCL_EXEC inline bool cellInside(Polygon tag, const CoordType& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  using T = ComponentType<CoordType>;
+  switch (tag.numberOfPoints())
+  {
+    case 3:
+      return cellInside(Triangle{}, pcoords);
+    case 4:
+      return cellInside(Quad{}, pcoords);
+    default:
+      break;
+  }
+
+  constexpr T epsilon = std::is_same<T, float>::value ? T(1e-5f) : T(1e-9f);
+
+  auto x = component(pcoords, 0) - T(0.5f);
+  auto y = component(pcoords, 1) - T(0.5f);
+  auto dist2 = (x * x) + (y * y);
+  if (dist2 > 0.25f) // definitely outside
+  {
+    return false;
+  }
+  else if (LCL_MATH_CALL(abs, (x)) < (T(4) * epsilon) && LCL_MATH_CALL(abs, (y)) < (T(4) * epsilon))
+  {
+    return true; // at the center
+  }
+  else
+  {
+    constexpr double two_pi = 2.0 * 3.14159265359;
+    T deltaAngle = static_cast<T>(two_pi) / static_cast<T>(tag.numberOfPoints());
+    T apothem = 0.5f * LCL_MATH_CALL(cos, (deltaAngle/2.0f));
+    if (dist2 <= (apothem * apothem)) // inside in-circle
+    {
+      return true;
+    }
+
+    // compute distance at which the line, from the center, through the given point, intersects
+    // the polygon edge
+    T angle = LCL_MATH_CALL(atan2, (y), (x));
+    if (angle < T(0))
+    {
+      angle += static_cast<T>(two_pi);
+    }
+
+    T a2 = angle - (LCL_MATH_CALL(floor, (angle / deltaAngle)) * deltaAngle);
+    T maxDist = apothem / LCL_MATH_CALL(cos, (LCL_MATH_CALL(abs, (deltaAngle/2.0f - a2))));
+    return dist2 <= (maxDist * maxDist);
+  }
+}
+
+namespace internal
+{
+
+template <typename CoordType>
+LCL_EXEC inline lcl::ErrorCode polygonToSubTrianglePCoords(
+  Polygon tag,
+  const CoordType& polygonPC,
+  IdComponent& p0,
+  IdComponent& p1,
+  ComponentType<CoordType> trianglePC[2]) noexcept
+{
+  using T = ComponentType<CoordType>;
+
+  constexpr T epsilon = std::is_same<T, float>::value ? T(1e-5f) : T(1e-9f);
+
+  // Find the sub-triangle containing pcoords
+  auto x = component(polygonPC, 0) - T(0.5f);
+  auto y = component(polygonPC, 1) - T(0.5f);
+  if (LCL_MATH_CALL(abs, (x)) < (T(4) * epsilon) && LCL_MATH_CALL(abs, (y)) < (T(4) * epsilon))
+  {
+    // we are at the center
+    p0 = 0;
+    p1 = 1;
+    trianglePC[0] = trianglePC[1] = T(0);
+    return ErrorCode::SUCCESS;
+  }
+
+  constexpr double two_pi = 2.0 * 3.14159265359;
+  T angle = LCL_MATH_CALL(atan2, (y), (x));
+  if (angle < T(0))
+  {
+    angle += static_cast<T>(two_pi);
+  }
+  T deltaAngle = static_cast<T>(two_pi) / static_cast<T>(tag.numberOfPoints());
+
+  p0 = static_cast<IdComponent>(LCL_MATH_CALL(floor, (angle / deltaAngle)));
+  p1 = (p0 + 1) % tag.numberOfPoints();
+
+  // Build triangle with polygon pcoords as its wcoords
+  T triPts[9] = { T(0), T(0), T(0), T(0), T(0), T(0), T(0), T(0), T(0) };
+  LCL_RETURN_ON_ERROR(parametricCenter(tag, triPts))
+  LCL_RETURN_ON_ERROR(parametricPoint(tag, p0, triPts + 3))
+  LCL_RETURN_ON_ERROR(parametricPoint(tag, p1, triPts + 6))
+
+  // Find the parametric coord on the triangle
+  T triWC[3] = { component(polygonPC, 0), component(polygonPC, 1), T(0) };
+  LCL_RETURN_ON_ERROR(worldToParametric(Triangle{}, makeFieldAccessorFlatSOAConst(triPts, 3), triWC, trianglePC))
+
+  return ErrorCode::SUCCESS;
+}
+
+template <typename Values>
+LCL_EXEC inline typename Values::ValueType polygonInterpolateComponentAtCenter(
+  Polygon tag, const Values& values, IdComponent comp) noexcept
+{
+  using T = internal::ClosestFloatType<typename Values::ValueType>;
+
+  auto weight = T{1} / static_cast<T>(tag.numberOfPoints());
+  auto result = static_cast<T>(values.getValue(0, comp));
+  for (IdComponent i = 1; i < tag.numberOfPoints(); ++i)
+  {
+    result += static_cast<T>(values.getValue(i, comp));
+  }
+  result *= weight;
+
+  return static_cast<typename Values::ValueType>(result);
+}
+
+} // namespace internal
+
+template <typename Values, typename CoordType, typename Result>
+LCL_EXEC lcl::ErrorCode interpolate(
+  Polygon tag,
+  const Values& values,
+  const CoordType& pcoords,
+  Result&& result) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  switch (tag.numberOfPoints())
+  {
+    case 3:
+      return interpolate(Triangle{}, values, pcoords, std::forward<Result>(result));
+    case 4:
+      return interpolate(Quad{}, values, pcoords, std::forward<Result>(result));
+    default:
+      break;
+  }
+
+  using ResultCompType = ComponentType<Result>;
+  using ProcessingType = internal::ClosestFloatType<typename Values::ValueType>;
+
+  IdComponent p0, p1;
+  ComponentType<CoordType> triPc[2];
+  LCL_RETURN_ON_ERROR(internal::polygonToSubTrianglePCoords(tag, pcoords, p0, p1, triPc))
+
+  // compute polygon interpolation from triangle weights
+  for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
+  {
+    ProcessingType triVals[3];
+    triVals[0] = static_cast<ProcessingType>(internal::polygonInterpolateComponentAtCenter(tag, values, c));
+    triVals[1] = static_cast<ProcessingType>(values.getValue(p0, c));
+    triVals[2] = static_cast<ProcessingType>(values.getValue(p1, c));
+    ResultCompType val = 0;
+    LCL_RETURN_ON_ERROR(interpolate(Triangle{}, makeFieldAccessorNestedSOA(triVals), triPc, &val))
+    component(result, c) = val;
+  }
+
+  return ErrorCode::SUCCESS;
+}
+
+namespace internal
+{
+// To find the gradient in a polygon (of 5 or more points), we will extract a small triangle near
+// the desired parameteric coordinates (pcoords). We return the field values (outField) and world
+// coordinates (outWCoords) for this triangle, which is all that is needed to find the gradient
+// in a triangle.
+//
+// The trangle will be "pointing" away from the center of the polygon, and pcoords will be placed
+// at the apex of the triangle. This is because if pcoords is at or near the edge of the polygon,
+// we do not want to push any of the points over the edge, and it is not trivial to determine
+// exactly where the edge of the polygon is.
+template <typename CoordType>
+LCL_EXEC inline void polygonGetTriangleAroundPCoords(
+  const CoordType& pcoords, ComponentType<CoordType> pc1[2], ComponentType<CoordType> pc2[2]) noexcept
+{
+  using T = ComponentType<CoordType>;
+
+  // Find the unit vector pointing from the center of the polygon to pcoords
+  Vector<T, 2> radialVector(component(pcoords, 0) - 0.5f, component(pcoords, 1) - 0.5f);
+  auto magSqr = dot(radialVector, radialVector);
+  if (magSqr > 8.0f * 1e-4f)
+  {
+    radialVector /= LCL_MATH_CALL(sqrt, (magSqr));
+  }
+  else
+  {
+    // pcoords is in the center of the polygon. Just point in an arbitrary direction
+    radialVector[0] = T(1);
+    radialVector[1] = T(0);
+  }
+
+  // We want the two points away from pcoords to be back toward the center but moved at 45 degrees
+  // off the radius. Simple geometry shows us that the (not quite unit) vectors of those two
+  // directions are (-radialVector[1] - radialVector[0], radialVector[0] - radialVector[1]) and
+  // (radialVector[1] - radialVector[0], -radialVector[0] - radialVector[1]).
+  //
+  //  *\ (-radialVector[1], radialVector[0])                                           //
+  //  |  \                                                                             //
+  //  |    \ (-radialVector[1] - radialVector[0], radialVector[0] - radialVector[1])   //
+  //  |      \                                                                         //
+  //  +-------* radialVector                                                           //
+  //  |      /                                                                         //
+  //  |    / (radialVector[1] - radialVector[0], -radialVector[0] - radialVector[1])   //
+  //  |  /                                                                             //
+  //  */ (radialVector[1], -radialVector[0])                                           //
+
+  // This scaling value is somewhat arbitrary. It is small enough to be "close" to the selected
+  // point and small enough to be guaranteed to be inside the polygon, but large enough to
+  // get an accurate gradient.
+  static constexpr T scale = 0.05f;
+
+  pc1[0] = pcoords[0] + scale * (-radialVector[1] - radialVector[0]);
+  pc1[1] = pcoords[1] + scale * (radialVector[0] - radialVector[1]);
+
+  pc2[0] = pcoords[0] + scale * (radialVector[1] - radialVector[0]);
+  pc2[1] = pcoords[1] + scale * (-radialVector[0] - radialVector[1]);
+}
+
+} // namespace internal
+
+template <typename Points, typename Values, typename CoordType, typename Result>
+LCL_EXEC inline lcl::ErrorCode derivative(
+  Polygon tag,
+  const Points& points,
+  const Values& values,
+  const CoordType& pcoords,
+  Result&& dx,
+  Result&& dy,
+  Result&& dz) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  switch (tag.numberOfPoints())
+  {
+    case 3:
+      return derivative(Triangle{},
+                        points,
+                        values,
+                        pcoords,
+                        std::forward<Result>(dx),
+                        std::forward<Result>(dy),
+                        std::forward<Result>(dz));
+    case 4:
+      return derivative(Quad{},
+                        points,
+                        values,
+                        pcoords,
+                        std::forward<Result>(dx),
+                        std::forward<Result>(dy),
+                        std::forward<Result>(dz));
+    default:
+      break;
+  }
+
+  using ResultCompType = ComponentType<Result>;
+  using ProcessingType = internal::ClosestFloatType<typename Values::ValueType>;
+
+  // Get the parametric coordinates of a small triangle, with pcoords as one of the vertices
+  ComponentType<CoordType> ptPc1[2], ptPc2[2];
+  internal::polygonGetTriangleAroundPCoords(pcoords, ptPc1, ptPc2);
+
+  // Compute world coordinates of the points of the triangle
+  internal::Vector<ProcessingType, 3> triPts[3];
+  LCL_RETURN_ON_ERROR(interpolate(tag, points, pcoords, triPts[0]))
+  LCL_RETURN_ON_ERROR(interpolate(tag, points, ptPc1, triPts[1]))
+  LCL_RETURN_ON_ERROR(interpolate(tag, points, ptPc2, triPts[2]))
+
+  // Compute the derivative on the triangle
+  //----------------------------------------
+  // 2-D coordinate system on the triangle's plane
+  internal::Space2D<ProcessingType> triSpace(triPts[0], triPts[1], triPts[2]);
+  internal::Vector<ProcessingType, 2> pts2d[3];
+  for (int i = 0; i < 3; ++i)
+  {
+    pts2d[i] = triSpace.to2DPoint(triPts[i]);
+  }
+
+  // pre-compute once
+  internal::Matrix<ProcessingType, 2, 2> jacobian;
+  internal::jacobian2D(Triangle{}, makeFieldAccessorNestedSOA(pts2d, 2), nullptr, jacobian);
+  internal::Matrix<ProcessingType, 2, 2> invJacobian;
+  LCL_RETURN_ON_ERROR(internal::matrixInverse(jacobian, invJacobian))
+
+  // Compute sub-triangle information of the three vertices of the derivation triangle to
+  // reduce the amount of redundant computations in the loop.
+  IdComponent subP1P2[3][2];
+  ComponentType<CoordType> pcs[3][2];
+  internal::polygonToSubTrianglePCoords(tag, pcoords, subP1P2[0][0], subP1P2[0][1], pcs[0]);
+  internal::polygonToSubTrianglePCoords(tag, ptPc1, subP1P2[1][0], subP1P2[1][1], pcs[1]);
+  internal::polygonToSubTrianglePCoords(tag, ptPc2, subP1P2[2][0], subP1P2[2][1], pcs[2]);
+
+  for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
+  {
+    // Interpolate component values at the vertices of the derivation triangle.
+    auto vCenter = static_cast<ProcessingType>(internal::polygonInterpolateComponentAtCenter(tag, values, c));
+    ProcessingType triVals[3];
+    for (int i = 0; i < 3; ++i)
+    {
+      ProcessingType field[3] = {vCenter,
+                                 static_cast<ProcessingType>(values.getValue(subP1P2[i][0], c)),
+                                 static_cast<ProcessingType>(values.getValue(subP1P2[i][1], c))};
+      LCL_RETURN_ON_ERROR(interpolate(Triangle{}, makeFieldAccessorNestedSOA(field), pcs[i], triVals + i))
+    }
+
+    // Compute derivative in the triangle
+    internal::Vector<ProcessingType, 2> dvdp;
+    parametricDerivative(Triangle{}, makeFieldAccessorNestedSOA(triVals), 0, nullptr, dvdp);
+    auto d2D = matrixMultiply(dvdp, invJacobian);
+    auto d3D = triSpace.to3DVec(d2D);
+
+    component(dx, c) = static_cast<ResultCompType>(d3D[0]);
+    component(dy, c) = static_cast<ResultCompType>(d3D[1]);
+    component(dz, c) = static_cast<ResultCompType>(d3D[2]);
+  }
+
+  return ErrorCode::SUCCESS;
+}
+
+template <typename Points, typename PCoordType, typename WCoordType>
+LCL_EXEC inline lcl::ErrorCode parametricToWorld(
+  Polygon tag,
+  const Points& points,
+  const PCoordType& pcoords,
+  WCoordType&& wcoords) noexcept
+{
+  return interpolate(tag, points, pcoords, std::forward<WCoordType>(wcoords));
+}
+
+template <typename Points, typename WCoordType, typename PCoordType>
+LCL_EXEC inline lcl::ErrorCode worldToParametric(
+  Polygon tag,
+  const Points& points,
+  const WCoordType& wcoords,
+  PCoordType&& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(PCoordType);
+
+  switch (tag.numberOfPoints())
+  {
+    case 3:
+      return worldToParametric(Triangle{}, points, wcoords, pcoords);
+    case 4:
+      return worldToParametric(Quad{}, points, wcoords, pcoords);
+    default:
+      break;
+  }
+
+  using T = ComponentType<WCoordType>;
+  auto numPoints = tag.numberOfPoints();
+
+  // Find the position of the center point.
+  internal::Vector<T, 3> wcoordCenter{T(0), T(0), T(0)};
+  for (IdComponent pointIndex = 0; pointIndex < numPoints; ++pointIndex)
+  {
+    wcoordCenter[0] += points.getValue(pointIndex, 0);
+    wcoordCenter[1] += points.getValue(pointIndex, 1);
+    wcoordCenter[2] += points.getValue(pointIndex, 2);
+  }
+  wcoordCenter /= static_cast<T>(numPoints);
+
+  // Find the normal vector to the polygon. If the polygon is planar, convex,
+  // and in general position, any three points will give a normal in the same
+  // direction. Although not perfectly robust, we can reduce the effect of
+  // non-planar, non-convex, or degenerate polygons by picking three points
+  // topologically far from each other. Note that we do not care about the
+  // length of the normal in this case.
+  internal::Vector<T, 3> polygonNormal;
+  {
+    internal::Vector<T, 3> v1p1, v1p2;
+    points.getTuple(0, v1p1);
+    points.getTuple(numPoints / 3, v1p2);
+
+    internal::Vector<T, 3> v2p1, v2p2;
+    points.getTuple(1, v2p1);
+    points.getTuple(2 * numPoints / 3, v2p2);
+
+    polygonNormal = internal::cross(v1p2 - v1p1, v2p2 - v2p1);
+  }
+
+  // Find which triangle wcoords is located in. We do this by defining the
+  // equations for the planes through the radial edges and perpendicular to the
+  // polygon. The point is in the triangle if it is on the correct side of both
+  // planes.
+  internal::Vector<T, 3> wc{ component(wcoords, 0), component(wcoords, 1), component(wcoords, 2) };
+  IdComponent firstPointIndex;
+  IdComponent secondPointIndex = 0;
+  internal::Vector<T, 3> firstPoint, secondPoint;
+  bool foundTriangle = false;
+  for (firstPointIndex = 0; firstPointIndex < numPoints - 1; ++firstPointIndex)
+  {
+    points.getTuple(firstPointIndex, firstPoint);
+    auto vecInPlane = firstPoint - wcoordCenter;
+
+    auto planeNormal = internal::cross(polygonNormal, vecInPlane);
+    auto planeOffset = internal::dot(planeNormal, wcoordCenter);
+    if (internal::dot(planeNormal, wc) < planeOffset)
+    {
+      // wcoords on wrong side of plane, thus outside of triangle
+      continue;
+    }
+
+    secondPointIndex = firstPointIndex + 1;
+    points.getTuple(secondPointIndex, secondPoint);
+    vecInPlane = secondPoint - wcoordCenter;
+
+    planeNormal = internal::cross(polygonNormal, vecInPlane);
+    planeOffset = internal::dot(planeNormal, wcoordCenter);
+    if (internal::dot(planeNormal, wc) > planeOffset)
+    {
+      // wcoords on wrong side of plane, thus outside of triangle
+      continue;
+    }
+
+    foundTriangle = true;
+    break;
+  }
+  if (!foundTriangle)
+  {
+    // wcoord was outside of all triangles we checked. It must be inside the
+    // one triangle we did not check (the one between the first and last
+    // polygon points).
+    firstPointIndex = numPoints - 1;
+    points.getTuple(firstPointIndex, firstPoint);
+    secondPointIndex = 0;
+    points.getTuple(secondPointIndex, secondPoint);
+  }
+
+  // Build a structure containing the points of the triangle wcoords is in and
+  // use the triangle version of this function to find the parametric
+  // coordinates.
+  internal::Vector<T, 3> triangleWCoords[3] = { wcoordCenter, firstPoint, secondPoint };
+  internal::Vector<T, 3> trianglePCoords;
+  LCL_RETURN_ON_ERROR(worldToParametric(
+    Triangle{}, makeFieldAccessorNestedSOA(triangleWCoords, 3), wc, trianglePCoords))
+
+  // trianglePCoords is in the triangle's parameter space rather than the
+  // polygon's parameter space. We can find the polygon's parameter space by
+  // repurposing parametricToWorld by using the
+  // polygon parametric coordinates as a proxy for world coordinates.
+  LCL_RETURN_ON_ERROR(parametricCenter(tag, triangleWCoords[0]))
+  LCL_RETURN_ON_ERROR(parametricPoint(tag, firstPointIndex, triangleWCoords[1]))
+  LCL_RETURN_ON_ERROR(parametricPoint(tag, secondPointIndex, triangleWCoords[2]))
+  triangleWCoords[0][2] = triangleWCoords[1][2] = triangleWCoords[2][2] = T(0);
+  LCL_RETURN_ON_ERROR(
+    parametricToWorld(Triangle{}, makeFieldAccessorNestedSOA(triangleWCoords, 3), trianglePCoords, wc))
+
+  component(pcoords, 0) = static_cast<ComponentType<PCoordType>>(wc[0]);
+  component(pcoords, 1) = static_cast<ComponentType<PCoordType>>(wc[1]);
+
+  return ErrorCode::SUCCESS;
+}
+
+} //namespace lcl
+
+#endif //lcl_Polygon_h

vtkm/thirdparty/lcl/vtkmlcl/lcl/Pyramid.h

@@ -0,0 +1,295 @@
+//============================================================================
+//  Copyright (c) Kitware, Inc.
+//  All rights reserved.
+//  See LICENSE.md 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.
+//============================================================================
+#ifndef lcl_Pyramid_h
+#define lcl_Pyramid_h
+
+#include <lcl/ErrorCode.h>
+#include <lcl/Shapes.h>
+
+#include <lcl/internal/Common.h>
+
+namespace lcl
+{
+
+class Pyramid : public Cell
+{
+public:
+  constexpr LCL_EXEC Pyramid() : Cell(ShapeId::PYRAMID, 5) {}
+  constexpr LCL_EXEC explicit Pyramid(const Cell& cell) : Cell(cell) {}
+};
+
+LCL_EXEC inline lcl::ErrorCode validate(Pyramid tag) noexcept
+{
+  if (tag.shape() != ShapeId::PYRAMID)
+  {
+    return ErrorCode::WRONG_SHAPE_ID_FOR_TAG_TYPE;
+  }
+  if (tag.numberOfPoints() != 5)
+  {
+    return ErrorCode::INVALID_NUMBER_OF_POINTS;
+  }
+
+  return ErrorCode::SUCCESS;
+}
+
+template<typename CoordType>
+LCL_EXEC inline lcl::ErrorCode parametricCenter(Pyramid, CoordType&& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  component(pcoords, 0) = 0.5f;
+  component(pcoords, 1) = 0.5f;
+  component(pcoords, 2) = 0.2f;
+  return ErrorCode::SUCCESS;
+}
+
+template<typename CoordType>
+LCL_EXEC inline lcl::ErrorCode parametricPoint(
+  Pyramid, IdComponent pointId, CoordType&& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  switch (pointId)
+  {
+    case 0:
+      component(pcoords, 0) = 0.0f;
+      component(pcoords, 1) = 0.0f;
+      component(pcoords, 2) = 0.0f;
+      break;
+    case 1:
+      component(pcoords, 0) = 1.0f;
+      component(pcoords, 1) = 0.0f;
+      component(pcoords, 2) = 0.0f;
+      break;
+    case 2:
+      component(pcoords, 0) = 1.0f;
+      component(pcoords, 1) = 1.0f;
+      component(pcoords, 2) = 0.0f;
+      break;
+    case 3:
+      component(pcoords, 0) = 0.0f;
+      component(pcoords, 1) = 1.0f;
+      component(pcoords, 2) = 0.0f;
+      break;
+    case 4:
+      component(pcoords, 0) = 0.5f;
+      component(pcoords, 1) = 0.5f;
+      component(pcoords, 2) = 1.0f;
+      break;
+    default:
+      return ErrorCode::INVALID_POINT_ID;
+  }
+
+  return ErrorCode::SUCCESS;
+}
+
+template<typename CoordType>
+LCL_EXEC inline ComponentType<CoordType> parametricDistance(Pyramid, const CoordType& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+  return internal::findParametricDistance(pcoords, 3);
+}
+
+template<typename CoordType>
+LCL_EXEC inline bool cellInside(Pyramid, const CoordType& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  using T = ComponentType<CoordType>;
+
+  constexpr T eps = 0.001f;
+  return component(pcoords, 0) >= -eps && component(pcoords, 0) <= (T{1} + eps) &&
+         component(pcoords, 1) >= -eps && component(pcoords, 1) <= (T{1} + eps) &&
+         component(pcoords, 2) >= -eps && component(pcoords, 2) <= (T{1} + eps);
+}
+
+template <typename Values, typename CoordType, typename Result>
+LCL_EXEC inline lcl::ErrorCode interpolate(
+  Pyramid,
+  const Values& values,
+  const CoordType& pcoords,
+  Result&& result) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  using T = internal::ClosestFloatType<typename Values::ValueType>;
+
+  for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
+  {
+    auto baseV0 = internal::lerp(static_cast<T>(values.getValue(0, c)),
+                                 static_cast<T>(values.getValue(1, c)),
+                                 static_cast<T>(component(pcoords, 0)));
+    auto baseV1 = internal::lerp(static_cast<T>(values.getValue(3, c)),
+                                 static_cast<T>(values.getValue(2, c)),
+                                 static_cast<T>(component(pcoords, 0)));
+
+    auto baseV = internal::lerp(baseV0, baseV1, static_cast<T>(component(pcoords, 1)));
+
+    auto v = internal::lerp(baseV,
+                            static_cast<T>(values.getValue(4, c)),
+                            static_cast<T>(component(pcoords, 2)));
+
+    component(result, c) = static_cast<ComponentType<Result>>(v);
+  }
+
+  return ErrorCode::SUCCESS;
+}
+
+namespace internal
+{
+
+template <typename Values, typename CoordType, typename Result>
+LCL_EXEC inline void parametricDerivative(Pyramid,
+                                           const Values& values,
+                                           IdComponent comp,
+                                           const CoordType& pcoords,
+                                           Result&& result) noexcept
+{
+  using T = internal::ClosestFloatType<typename Values::ValueType>;
+  T p0 = static_cast<T>(component(pcoords, 0));
+  T p1 = static_cast<T>(component(pcoords, 1));
+  T p2 = static_cast<T>(component(pcoords, 2));
+  T rm = T{1} - p0;
+  T sm = T{1} - p1;
+  T tm = T{1} - p2;
+
+  T dr = (static_cast<T>(values.getValue(0, comp)) * -sm * tm) +
+         (static_cast<T>(values.getValue(1, comp)) *  sm * tm) +
+         (static_cast<T>(values.getValue(2, comp)) *  p1 * tm) +
+         (static_cast<T>(values.getValue(3, comp)) * -p1 * tm);
+
+  T ds = (static_cast<T>(values.getValue(0, comp)) * -rm * tm) +
+         (static_cast<T>(values.getValue(1, comp)) * -p0 * tm) +
+         (static_cast<T>(values.getValue(2, comp)) *  p0 * tm) +
+         (static_cast<T>(values.getValue(3, comp)) *  rm * tm);
+
+  T dt = (static_cast<T>(values.getValue(0, comp)) * -rm * sm) +
+         (static_cast<T>(values.getValue(1, comp)) * -p0 * sm) +
+         (static_cast<T>(values.getValue(2, comp)) * -p0 * p1) +
+         (static_cast<T>(values.getValue(3, comp)) * -rm * p1) +
+         (static_cast<T>(values.getValue(4, comp)));
+
+  component(result, 0) = static_cast<ComponentType<Result>>(dr);
+  component(result, 1) = static_cast<ComponentType<Result>>(ds);
+  component(result, 2) = static_cast<ComponentType<Result>>(dt);
+}
+
+} // internal
+
+template <typename Points, typename Values, typename CoordType, typename Result>
+LCL_EXEC inline lcl::ErrorCode derivative(
+  Pyramid,
+  const Points& points,
+  const Values& values,
+  const CoordType& pcoords,
+  Result&& dx,
+  Result&& dy,
+  Result&& dz) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  using ProcessingType = internal::ClosestFloatType<typename Values::ValueType>;
+  using ResultCompType = ComponentType<Result>;
+
+  if (component(pcoords, 2) > ComponentType<CoordType>(.999f))
+  {
+    // If we are at the apex of the pyramid we need to do something special.
+    // As we approach the apex, the derivatives of the parametric shape
+    // functions in x and y go to 0 while the inverse of the Jacobian
+    // also goes to 0.  This results in 0/0 but using l'Hopital's rule
+    // we could actually compute the value of the limit, if we had a
+    // functional expression to compute the gradient.  We're on a computer
+    // so we don't but we can cheat and do a linear extrapolation of the
+    // derivatives which really ends up as the same thing.
+    internal::Matrix<ProcessingType, 3, 3> j, ij1, ij2;
+
+    const ComponentType<CoordType> pc1[3] = {0.5f, 0.5f, (2.0f * 0.998f) - component(pcoords, 2)};
+    internal::jacobian3D(Pyramid{}, points, pc1, j);
+    LCL_RETURN_ON_ERROR(internal::matrixInverse(j, ij1))
+
+    const ComponentType<CoordType> pc2[3] = {0.5f, 0.5f, 0.998f};
+    internal::jacobian3D(Pyramid{}, points, pc2, j);
+    LCL_RETURN_ON_ERROR(internal::matrixInverse(j, ij2))
+
+    for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
+    {
+      internal::Vector<ProcessingType, 3> dvdp;
+
+      internal::parametricDerivative(Pyramid{}, values, c, pc1, dvdp);
+      auto d1 = matrixMultiply(dvdp, ij1);
+
+      internal::parametricDerivative(Pyramid{}, values, c, pc2, dvdp);
+      auto d2 = matrixMultiply(dvdp, ij2);
+
+      component(dx, c) = static_cast<ResultCompType>((d2[0] * 2.0f) - d1[0]);
+      component(dy, c) = static_cast<ResultCompType>((d2[1] * 2.0f) - d1[1]);
+      component(dz, c) = static_cast<ResultCompType>((d2[2] * 2.0f) - d1[2]);
+    }
+  }
+  else
+  {
+    return internal::derivative3D(Pyramid{},
+                                  points,
+                                  values,
+                                  pcoords,
+                                  std::forward<Result>(dx),
+                                  std::forward<Result>(dy),
+                                  std::forward<Result>(dz));
+  }
+
+  return ErrorCode::SUCCESS;
+}
+
+template <typename Points, typename PCoordType, typename WCoordType>
+LCL_EXEC inline lcl::ErrorCode parametricToWorld(
+  Pyramid,
+  const Points& points,
+  const PCoordType& pcoords,
+  WCoordType&& wcoords) noexcept
+{
+  return interpolate(Pyramid{}, points, pcoords, std::forward<WCoordType>(wcoords));
+}
+
+template <typename Points, typename WCoordType, typename PCoordType>
+LCL_EXEC inline lcl::ErrorCode worldToParametric(
+  Pyramid,
+  const Points& points,
+  const WCoordType& wcoords,
+  PCoordType&& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(PCoordType);
+
+  using TIn = typename Points::ValueType;
+  using TOut = ComponentType<PCoordType>;
+
+  internal::Vector<TIn, 3> wcVec{component(wcoords, 0), component(wcoords, 1), component(wcoords, 2)};
+
+  // Newton's method fails if the wcoord is too close to the apex. Just return the pcoords at the
+  // apex for those cases.
+  internal::Vector<TOut, 3> pcBaseCenter(0.5f, 0.5f, 0.0f);
+  internal::Vector<TIn, 3> apex, wcBaseCenter;
+  points.getTuple(4, apex);
+  LCL_RETURN_ON_ERROR(parametricToWorld(Pyramid{}, points, pcBaseCenter, wcBaseCenter))
+  auto apexToBase = wcBaseCenter - apex;
+  auto apexToWc = wcVec - apex;
+  auto dist2ApexToBase = internal::dot(apexToBase, apexToBase);
+  auto dist2ApexToWC = internal::dot(apexToWc, apexToWc);
+  if (dist2ApexToWC <= (1e-6f * dist2ApexToBase))
+  {
+    return parametricPoint(Pyramid{}, 4, pcoords);
+  }
+
+  return internal::worldToParametric3D(
+    Pyramid{}, points, wcoords, std::forward<PCoordType>(pcoords));
+}
+
+} // lcl
+
+#endif // lcl_Pyramid_h

vtkm/thirdparty/lcl/vtkmlcl/lcl/Quad.h

@@ -0,0 +1,197 @@
+//============================================================================
+//  Copyright (c) Kitware, Inc.
+//  All rights reserved.
+//  See LICENSE.md 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.
+//============================================================================
+#ifndef lcl_Quad_h
+#define lcl_Quad_h
+
+#include <lcl/ErrorCode.h>
+#include <lcl/Shapes.h>
+
+#include <lcl/internal/Common.h>
+
+namespace lcl
+{
+
+class Quad : public Cell
+{
+public:
+  constexpr LCL_EXEC Quad() : Cell(ShapeId::QUAD, 4) {}
+  constexpr LCL_EXEC explicit Quad(const Cell& cell) : Cell(cell) {}
+};
+
+LCL_EXEC inline lcl::ErrorCode validate(Quad tag) noexcept
+{
+  if (tag.shape() != ShapeId::QUAD && tag.shape() != ShapeId::PIXEL)
+  {
+    return ErrorCode::WRONG_SHAPE_ID_FOR_TAG_TYPE;
+  }
+  if (tag.numberOfPoints() != 4)
+  {
+    return ErrorCode::INVALID_NUMBER_OF_POINTS;
+  }
+
+  return ErrorCode::SUCCESS;
+}
+
+template<typename CoordType>
+LCL_EXEC inline lcl::ErrorCode parametricCenter(Quad, CoordType&& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  component(pcoords, 0) = 0.5f;
+  component(pcoords, 1) = 0.5f;
+  return ErrorCode::SUCCESS;
+}
+
+template<typename CoordType>
+LCL_EXEC inline lcl::ErrorCode parametricPoint(
+  Quad, IdComponent pointId, CoordType&& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  switch (pointId)
+  {
+    case 0:
+      component(pcoords, 0) = 0.0f;
+      component(pcoords, 1) = 0.0f;
+      break;
+    case 1:
+      component(pcoords, 0) = 1.0f;
+      component(pcoords, 1) = 0.0f;
+      break;
+    case 2:
+      component(pcoords, 0) = 1.0f;
+      component(pcoords, 1) = 1.0f;
+      break;
+    case 3:
+      component(pcoords, 0) = 0.0f;
+      component(pcoords, 1) = 1.0f;
+      break;
+    default:
+      return ErrorCode::INVALID_POINT_ID;
+  }
+
+  return ErrorCode::SUCCESS;
+}
+
+template<typename CoordType>
+LCL_EXEC inline ComponentType<CoordType> parametricDistance(Quad, const CoordType& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+  return internal::findParametricDistance(pcoords, 2);
+}
+
+template<typename CoordType>
+LCL_EXEC inline bool cellInside(Quad, const CoordType& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  using T = ComponentType<CoordType>;
+
+  constexpr T eps = 0.001f;
+  return component(pcoords, 0) >= -eps && component(pcoords, 0) <= (T{1} + eps) &&
+         component(pcoords, 1) >= -eps && component(pcoords, 1) <= (T{1} + eps);
+}
+
+template <typename Values, typename CoordType, typename Result>
+LCL_EXEC inline lcl::ErrorCode interpolate(
+  Quad,
+  const Values& values,
+  const CoordType& pcoords,
+  Result&& result) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  using T = internal::ClosestFloatType<typename Values::ValueType>;
+
+  for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
+  {
+    auto v0 = internal::lerp(static_cast<T>(values.getValue(0, c)),
+                             static_cast<T>(values.getValue(1, c)),
+                             static_cast<T>(component(pcoords, 0)));
+    auto v1 = internal::lerp(static_cast<T>(values.getValue(3, c)),
+                             static_cast<T>(values.getValue(2, c)),
+                             static_cast<T>(component(pcoords, 0)));
+    auto v = internal::lerp(v0, v1, static_cast<T>(component(pcoords, 1)));
+    component(result, c) = static_cast<ComponentType<Result>>(v);
+  }
+
+  return ErrorCode::SUCCESS;
+}
+
+namespace internal
+{
+
+template <typename Values, typename CoordType, typename Result>
+LCL_EXEC inline void parametricDerivative(
+  Quad, const Values& values, IdComponent comp, const CoordType& pcoords, Result&& result) noexcept
+{
+  using T = internal::ClosestFloatType<typename Values::ValueType>;
+  T p0 = static_cast<T>(component(pcoords, 0));
+  T p1 = static_cast<T>(component(pcoords, 1));
+  T rm = T{1} - p0;
+  T sm = T{1} - p1;
+
+  T dr = (static_cast<T>(values.getValue(0, comp)) * -sm) +
+         (static_cast<T>(values.getValue(1, comp)) *  sm) +
+         (static_cast<T>(values.getValue(2, comp)) *  p1) +
+         (static_cast<T>(values.getValue(3, comp)) * -p1);
+  T ds = (static_cast<T>(values.getValue(0, comp)) * -rm) +
+         (static_cast<T>(values.getValue(1, comp)) * -p0) +
+         (static_cast<T>(values.getValue(2, comp)) *  p0) +
+         (static_cast<T>(values.getValue(3, comp)) *  rm);
+
+  component(result, 0) = static_cast<ComponentType<Result>>(dr);
+  component(result, 1) = static_cast<ComponentType<Result>>(ds);
+}
+
+} // internal
+
+template <typename Points, typename Values, typename CoordType, typename Result>
+LCL_EXEC inline lcl::ErrorCode derivative(
+  Quad,
+  const Points& points,
+  const Values& values,
+  const CoordType& pcoords,
+  Result&& dx,
+  Result&& dy,
+  Result&& dz) noexcept
+{
+  return internal::derivative2D(Quad{},
+                                points,
+                                values,
+                                pcoords,
+                                std::forward<Result>(dx),
+                                std::forward<Result>(dy),
+                                std::forward<Result>(dz));
+}
+
+template <typename Points, typename PCoordType, typename WCoordType>
+LCL_EXEC inline lcl::ErrorCode parametricToWorld(
+  Quad,
+  const Points& points,
+  const PCoordType& pcoords,
+  WCoordType&& wcoords) noexcept
+{
+  return interpolate(Quad{}, points, pcoords, std::forward<WCoordType>(wcoords));
+}
+
+template <typename Points, typename WCoordType, typename PCoordType>
+LCL_EXEC inline lcl::ErrorCode worldToParametric(
+  Quad,
+  const Points& points,
+  const WCoordType& wcoords,
+  PCoordType&& pcoords) noexcept
+{
+  return internal::worldToParametric2D(Quad{}, points, wcoords, std::forward<PCoordType>(pcoords));
+}
+
+} //namespace lcl
+
+#endif //lcl_Quad_h

vtkm/thirdparty/lcl/vtkmlcl/lcl/Shapes.h

@@ -0,0 +1,142 @@
+//============================================================================
+//  Copyright (c) Kitware, Inc.
+//  All rights reserved.
+//  See LICENSE.md 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.
+//============================================================================
+#ifndef lcl_Shapes_h
+#define lcl_Shapes_h
+
+#include <lcl/internal/Config.h>
+
+#include <cstdint>
+
+namespace lcl
+{
+
+enum ShapeId : IdShape
+{
+  // Linear cells
+  EMPTY            = 0,
+  VERTEX           = 1,
+  //POLY_VERTEX      = 2,
+  LINE             = 3,
+  //POLY_LINE        = 4,
+  TRIANGLE         = 5,
+  //TRIANGLE_STRIP   = 6,
+  POLYGON          = 7,
+  PIXEL            = 8,
+  QUAD             = 9,
+  TETRA            = 10,
+  VOXEL            = 11,
+  HEXAHEDRON       = 12,
+  WEDGE            = 13,
+  PYRAMID          = 14,
+
+  NUMBER_OF_CELL_SHAPES
+};
+
+class Cell
+{
+public:
+  constexpr LCL_EXEC Cell() : Shape(ShapeId::EMPTY), NumberOfPoints(0) {}
+  constexpr LCL_EXEC Cell(IdShape shape, IdComponent numPoints)
+    : Shape(shape), NumberOfPoints(numPoints)
+  {
+  }
+
+  constexpr LCL_EXEC IdShape shape() const noexcept { return this->Shape; }
+  constexpr LCL_EXEC IdComponent numberOfPoints() const noexcept { return this->NumberOfPoints; }
+
+protected:
+  IdShape Shape;
+  IdComponent NumberOfPoints;
+};
+
+/// \brief Check if a shape id is valid
+/// \param[in]  shapeId  The id to check.
+/// \return              true if id is a valid shape id.
+///
+constexpr inline LCL_EXEC bool isValidShape(IdShape shapeId)
+{
+  return (shapeId >= ShapeId::EMPTY) && (shapeId < ShapeId::NUMBER_OF_CELL_SHAPES);
+}
+
+/// \brief Returns the dimensionality of a cell
+/// \param[in]  shapeId  The shape id of the cell.
+/// \return              The dimensionality of the cell, -1 for invalid shapes.
+///
+inline LCL_EXEC int dimension(IdShape shapeId)
+{
+  switch (shapeId)
+  {
+    case VERTEX:
+      return 0;
+    case LINE:
+      return 1;
+    case TRIANGLE:
+    case POLYGON:
+    case PIXEL:
+    case QUAD:
+      return 2;
+    case TETRA:
+    case VOXEL:
+    case HEXAHEDRON:
+    case WEDGE:
+    case PYRAMID:
+      return 3;
+    case EMPTY:
+    default:
+      return -1;
+  }
+}
+
+/// \brief Returns the dimensionality of a cell
+/// \param[in]  cell  The cell.
+/// \return           The dimensionality of the cell, -1 for invalid shapes.
+///
+inline LCL_EXEC int dimension(Cell cell)
+{
+  return dimension(cell.shape());
+}
+
+// forward declare cell tags
+class Vertex;
+class Line;
+class Triangle;
+class Polygon;
+class Pixel;
+class Quad;
+class Tetra;
+class Voxel;
+class Hexahedron;
+class Wedge;
+class Pyramid;
+
+} //namespace lcl
+
+#define lclGenericCellShapeMacroCase(cellId, cell, call)                                          \
+  case cellId:                                                                                     \
+  {                                                                                                \
+    using CellTag = cell;                                                                          \
+    call;                                                                                          \
+  }                                                                                                \
+  break
+
+#define lclGenericCellShapeMacro(call)                                                            \
+  lclGenericCellShapeMacroCase(lcl::ShapeId::VERTEX,     lcl::Vertex,     call);                \
+  lclGenericCellShapeMacroCase(lcl::ShapeId::LINE,       lcl::Line,       call);                \
+  lclGenericCellShapeMacroCase(lcl::ShapeId::TRIANGLE,   lcl::Triangle,   call);                \
+  lclGenericCellShapeMacroCase(lcl::ShapeId::POLYGON,    lcl::Polygon,    call);                \
+  lclGenericCellShapeMacroCase(lcl::ShapeId::PIXEL,      lcl::Pixel,      call);                \
+  lclGenericCellShapeMacroCase(lcl::ShapeId::QUAD,       lcl::Quad,       call);                \
+  lclGenericCellShapeMacroCase(lcl::ShapeId::TETRA,      lcl::Tetra,      call);                \
+  lclGenericCellShapeMacroCase(lcl::ShapeId::VOXEL,      lcl::Voxel,      call);                \
+  lclGenericCellShapeMacroCase(lcl::ShapeId::HEXAHEDRON, lcl::Hexahedron, call);                \
+  lclGenericCellShapeMacroCase(lcl::ShapeId::WEDGE,      lcl::Wedge,      call);                \
+  lclGenericCellShapeMacroCase(lcl::ShapeId::PYRAMID,    lcl::Pyramid,    call)
+
+#endif //lcl_Shapes_h

vtkm/thirdparty/lcl/vtkmlcl/lcl/Tetra.h

@@ -0,0 +1,219 @@
+//============================================================================
+//  Copyright (c) Kitware, Inc.
+//  All rights reserved.
+//  See LICENSE.md 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.
+//============================================================================
+#ifndef lcl_Tetra_h
+#define lcl_Tetra_h
+
+#include <lcl/ErrorCode.h>
+#include <lcl/Shapes.h>
+
+#include <lcl/internal/Common.h>
+
+namespace lcl
+{
+
+class Tetra : public Cell
+{
+public:
+  constexpr LCL_EXEC Tetra() : Cell(ShapeId::TETRA, 4) {}
+  constexpr LCL_EXEC explicit Tetra(const Cell& cell) : Cell(cell) {}
+};
+
+LCL_EXEC inline lcl::ErrorCode validate(Tetra tag) noexcept
+{
+  if (tag.shape() != ShapeId::TETRA)
+  {
+    return ErrorCode::WRONG_SHAPE_ID_FOR_TAG_TYPE;
+  }
+  if (tag.numberOfPoints() != 4)
+  {
+    return ErrorCode::INVALID_NUMBER_OF_POINTS;
+  }
+
+  return ErrorCode::SUCCESS;
+}
+
+template<typename CoordType>
+LCL_EXEC inline lcl::ErrorCode parametricCenter(Tetra, CoordType&& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  component(pcoords, 0) = 0.25f;
+  component(pcoords, 1) = 0.25f;
+  component(pcoords, 2) = 0.25f;
+  return ErrorCode::SUCCESS;
+}
+
+template<typename CoordType>
+LCL_EXEC inline lcl::ErrorCode parametricPoint(
+  Tetra, IdComponent pointId, CoordType&& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  switch (pointId)
+  {
+    case 0:
+      component(pcoords, 0) = 0.0f;
+      component(pcoords, 1) = 0.0f;
+      component(pcoords, 2) = 0.0f;
+      break;
+    case 1:
+      component(pcoords, 0) = 1.0f;
+      component(pcoords, 1) = 0.0f;
+      component(pcoords, 2) = 0.0f;
+      break;
+    case 2:
+      component(pcoords, 0) = 0.0f;
+      component(pcoords, 1) = 1.0f;
+      component(pcoords, 2) = 0.0f;
+      break;
+    case 3:
+      component(pcoords, 0) = 0.0f;
+      component(pcoords, 1) = 0.0f;
+      component(pcoords, 2) = 1.0f;
+      break;
+    default:
+      return ErrorCode::INVALID_POINT_ID;
+  }
+
+  return ErrorCode::SUCCESS;
+}
+
+template<typename CoordType>
+LCL_EXEC inline ComponentType<CoordType> parametricDistance(Tetra, const CoordType& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  ComponentType<CoordType> weights[4];
+  weights[0] = ComponentType<CoordType>{1} - component(pcoords, 0) - component(pcoords, 1) - component(pcoords, 2);
+  weights[1] = component(pcoords, 0);
+  weights[2] = component(pcoords, 1);
+  weights[3] = component(pcoords, 2);
+  return internal::findParametricDistance(weights, 4);
+}
+
+template<typename CoordType>
+LCL_EXEC inline bool cellInside(Tetra, const CoordType& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  using T = ComponentType<CoordType>;
+
+  constexpr T eps = 0.001f;
+  return component(pcoords, 0) >= -eps &&
+         component(pcoords, 1) >= -eps &&
+         component(pcoords, 2) >= -eps &&
+         (component(pcoords, 0) + component(pcoords, 1) + component(pcoords, 2)) <= (T{1} + eps);
+}
+
+template <typename Values, typename CoordType, typename Result>
+LCL_EXEC inline lcl::ErrorCode interpolate(
+  Tetra,
+  const Values& values,
+  const CoordType& pcoords,
+  Result&& result) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  using T = internal::ClosestFloatType<typename Values::ValueType>;
+
+  auto w0 = T(1) - static_cast<T>(component(pcoords, 0) + component(pcoords, 1) + component(pcoords, 2));
+  auto w1 = static_cast<T>(component(pcoords, 0));
+  auto w2 = static_cast<T>(component(pcoords, 1));
+  auto w3 = static_cast<T>(component(pcoords, 2));
+
+  for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
+  {
+    auto v = static_cast<T>(values.getValue(0, c)) * w0 +
+             static_cast<T>(values.getValue(1, c)) * w1 +
+             static_cast<T>(values.getValue(2, c)) * w2 +
+             static_cast<T>(values.getValue(3, c)) * w3;
+    component(result, c) = static_cast<ComponentType<Result>>(v);
+  }
+  return ErrorCode::SUCCESS;
+}
+
+namespace internal
+{
+
+template <typename Values, typename CoordType, typename Result>
+LCL_EXEC inline void parametricDerivative(
+  Tetra, const Values& values, IdComponent comp, const CoordType&, Result&& result) noexcept
+{
+  component(result, 0) = static_cast<ComponentType<Result>>(values.getValue(1, comp) -
+                                                            values.getValue(0, comp));
+  component(result, 1) = static_cast<ComponentType<Result>>(values.getValue(2, comp) -
+                                                            values.getValue(0, comp));
+  component(result, 2) = static_cast<ComponentType<Result>>(values.getValue(3, comp) -
+                                                            values.getValue(0, comp));
+}
+
+} // internal
+
+template <typename Points, typename Values, typename CoordType, typename Result>
+LCL_EXEC inline lcl::ErrorCode derivative(
+  Tetra,
+  const Points& points,
+  const Values& values,
+  const CoordType& pcoords,
+  Result&& dx,
+  Result&& dy,
+  Result&& dz) noexcept
+{
+  return internal::derivative3D(Tetra{},
+                                points,
+                                values,
+                                pcoords,
+                                std::forward<Result>(dx),
+                                std::forward<Result>(dy),
+                                std::forward<Result>(dz));
+}
+
+template <typename Points, typename PCoordType, typename WCoordType>
+LCL_EXEC inline lcl::ErrorCode parametricToWorld(
+  Tetra,
+  const Points& points,
+  const PCoordType& pcoords,
+  WCoordType&& wcoords) noexcept
+{
+  return interpolate(Tetra{}, points, pcoords, std::forward<WCoordType>(wcoords));
+}
+
+template <typename Points, typename WCoordType, typename PCoordType>
+LCL_EXEC inline lcl::ErrorCode worldToParametric(
+  Tetra,
+  const Points& points,
+  const WCoordType& wcoords,
+  PCoordType&& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(PCoordType);
+
+  using T = ComponentType<PCoordType>;
+
+  internal::Matrix<T, 3, 3> A;
+  internal::Vector<T, 3> b, x;
+  for (int i = 0; i < 3; ++i)
+  {
+    for (int j = 0; j < 3; ++j)
+    {
+      A(j, i) = static_cast<T>(points.getValue(i + 1, j) - points.getValue(0, j));
+    }
+    b[i] = static_cast<T>(component(wcoords, i) - points.getValue(0, i));
+  }
+  LCL_RETURN_ON_ERROR(internal::solveLinearSystem(A, b, x))
+
+  component(pcoords, 0) = x[0];
+  component(pcoords, 1) = x[1];
+  component(pcoords, 2) = x[2];
+  return ErrorCode::SUCCESS;
+}
+
+} // lcl
+
+#endif // lcl_Tetra_h

vtkm/thirdparty/lcl/vtkmlcl/lcl/Triangle.h

@@ -0,0 +1,271 @@
+//============================================================================
+//  Copyright (c) Kitware, Inc.
+//  All rights reserved.
+//  See LICENSE.md 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.
+//============================================================================
+#ifndef lcl_Triangle_h
+#define lcl_Triangle_h
+
+#include <lcl/ErrorCode.h>
+#include <lcl/Shapes.h>
+
+#include <lcl/internal/Common.h>
+
+namespace lcl
+{
+
+class Triangle : public Cell
+{
+public:
+  constexpr LCL_EXEC Triangle() : Cell(ShapeId::TRIANGLE, 3) {}
+  constexpr LCL_EXEC explicit Triangle(const Cell& cell) : Cell(cell) {}
+};
+
+LCL_EXEC inline lcl::ErrorCode validate(Triangle tag) noexcept
+{
+  if (tag.shape() != ShapeId::TRIANGLE)
+  {
+    return ErrorCode::WRONG_SHAPE_ID_FOR_TAG_TYPE;
+  }
+  if (tag.numberOfPoints() != 3)
+  {
+    return ErrorCode::INVALID_NUMBER_OF_POINTS;
+  }
+
+  return ErrorCode::SUCCESS;
+}
+
+template<typename CoordType>
+LCL_EXEC inline lcl::ErrorCode parametricCenter(Triangle, CoordType&& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  using T = ComponentType<CoordType>;
+  component(pcoords, 0) = T(1)/T(3);
+  component(pcoords, 1) = T(1)/T(3);
+  component(pcoords, 2) = T(0);
+  return ErrorCode::SUCCESS;
+}
+
+template<typename CoordType>
+LCL_EXEC inline lcl::ErrorCode parametricPoint(
+  Triangle, IdComponent pointId, CoordType&& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  component(pcoords, 2) = 0.0f;
+  switch (pointId)
+  {
+    case 0:
+      component(pcoords, 0) = 0.0f;
+      component(pcoords, 1) = 0.0f;
+      break;
+    case 1:
+      component(pcoords, 0) = 1.0f;
+      component(pcoords, 1) = 0.0f;
+      break;
+    case 2:
+      component(pcoords, 0) = 0.0f;
+      component(pcoords, 1) = 1.0f;
+      break;
+    default:
+      return ErrorCode::INVALID_POINT_ID;
+  }
+
+  return ErrorCode::SUCCESS;
+}
+
+template<typename CoordType>
+LCL_EXEC inline ComponentType<CoordType> parametricDistance(Triangle, const CoordType& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  ComponentType<CoordType> weights[3];
+  weights[0] = ComponentType<CoordType>{1} - component(pcoords, 0) - component(pcoords, 1);
+  weights[1] = component(pcoords, 0);
+  weights[2] = component(pcoords, 1);
+  return internal::findParametricDistance(weights, 3);
+}
+
+template<typename CoordType>
+LCL_EXEC inline bool cellInside(Triangle, const CoordType& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  using T = ComponentType<CoordType>;
+  return component(pcoords, 0) >= T{0} && component(pcoords, 1) >= T{0} &&
+         (component(pcoords, 0) + component(pcoords, 1)) <= T{1};
+}
+
+template <typename Values, typename CoordType, typename Result>
+LCL_EXEC inline lcl::ErrorCode interpolate(
+  Triangle,
+  const Values& values,
+  const CoordType& pcoords,
+  Result&& result) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  using T = internal::ClosestFloatType<typename Values::ValueType>;
+
+  auto w0 = T(1) - static_cast<T>(component(pcoords, 0) + component(pcoords, 1));
+  auto w1 = static_cast<T>(component(pcoords, 0));
+  auto w2 = static_cast<T>(component(pcoords, 1));
+
+  for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
+  {
+    auto v = static_cast<T>(values.getValue(0, c)) * w0 +
+             static_cast<T>(values.getValue(1, c)) * w1 +
+             static_cast<T>(values.getValue(2, c)) * w2;
+    component(result, c) = static_cast<ComponentType<Result>>(v);
+  }
+  return ErrorCode::SUCCESS;
+}
+
+namespace internal
+{
+
+template <typename Values, typename CoordType, typename Result>
+LCL_EXEC inline void parametricDerivative(
+  Triangle, const Values& values, IdComponent comp, const CoordType&, Result&& result) noexcept
+{
+  component(result, 0) = static_cast<ComponentType<Result>>(values.getValue(1, comp) -
+                                                            values.getValue(0, comp));
+  component(result, 1) = static_cast<ComponentType<Result>>(values.getValue(2, comp) -
+                                                            values.getValue(0, comp));
+}
+
+} // internal
+
+template <typename Points, typename Values, typename CoordType, typename Result>
+LCL_EXEC inline lcl::ErrorCode derivative(
+  Triangle,
+  const Points& points,
+  const Values& values,
+  const CoordType& pcoords,
+  Result&& dx,
+  Result&& dy,
+  Result&& dz) noexcept
+{
+  return internal::derivative2D(Triangle{},
+                                points,
+                                values,
+                                pcoords,
+                                std::forward<Result>(dx),
+                                std::forward<Result>(dy),
+                                std::forward<Result>(dz));
+}
+
+template <typename Points, typename PCoordType, typename WCoordType>
+LCL_EXEC inline lcl::ErrorCode parametricToWorld(
+  Triangle,
+  const Points& points,
+  const PCoordType& pcoords,
+  WCoordType&& wcoords) noexcept
+{
+  return interpolate(Triangle{}, points, pcoords, std::forward<WCoordType>(wcoords));
+}
+
+//-----------------------------------------------------------------------------
+// The following implementation is lifted from VTK-m
+//-----------------------------------------------------------------------------
+// We will solve the world to parametric coordinates problem geometrically.
+// Consider the parallelogram formed by wcoords and p0 of the triangle and
+// the two adjacent edges. This parallelogram is equivalent to the
+// axis-aligned rectangle anchored at the origin of parametric space.
+//
+//   p2 |\                 (1,0) |\                                        //
+//      | \                      |  \                                      //
+//      |  \                     |    \                                    //
+//     |    \                    |      \                                  //
+//     |     \                   |        \                                //
+//     |      \                  |    (u,v) \                              //
+//    | ---    \                 |-------*    \                            //
+//    |    ---*wcoords           |       |      \                          //
+//    |       |  \               |       |        \                        //
+// p0 *---    |   \        (0,0) *------------------\ (1,0)                //
+//        ---|     \                                                       //
+//           x--    \                                                      //
+//              ---  \                                                     //
+//                 ---\ p1                                                 //
+//
+// In this diagram, the distance between p0 and the point marked x divided by
+// the length of the edge it is on is equal, by proportionality, to the u
+// parametric coordinate. (The v coordinate follows the other edge
+// accordingly.) Thus, if we can find the intersection at x (or more
+// specifically the distance between p0 and x), then we can find that
+// parametric coordinate.
+//
+// Because the triangle is in 3-space, we are actually going to intersect the
+// edge with a plane that is parallel to the opposite edge of p0 and
+// perpendicular to the triangle. This is partially because it is easy to
+// find the intersection between a plane and a line and partially because the
+// computation will work for points not on the plane. (The result is
+// equivalent to a point projected on the plane.)
+//
+// First, we define an implicit plane as:
+//
+// Dot((p - wcoords), planeNormal) = 0
+//
+// where planeNormal is the normal to the plane (easily computed from the
+// triangle), and p is any point in the plane. Next, we define the parametric
+// form of the line:
+//
+// p(d) = (p1 - p0)d + p0
+//
+// Where d is the fraction of distance from p0 toward p1. Note that d is
+// actually equal to the parametric coordinate we are trying to find. Once we
+// compute it, we are done. We can skip the part about finding the actual
+// coordinates of the intersection.
+//
+// Solving for the intersection is as simple as substituting the line's
+// definition of p(d) into p for the plane equation. With some basic algebra
+// you get:
+//
+// d = Dot((wcoords - p0), planeNormal)/Dot((p1-p0), planeNormal)
+//
+// From here, the u coordinate is simply d. The v coordinate follows
+// similarly.
+//
+template <typename Points, typename WCoordType, typename PCoordType>
+LCL_EXEC inline lcl::ErrorCode worldToParametric(
+  Triangle,
+  const Points& points,
+  const WCoordType& wcoords,
+  PCoordType&& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(PCoordType);
+
+  using TIn = typename Points::ValueType;
+  using TOut = ComponentType<PCoordType>;
+
+  internal::Vector<TIn, 3> pts[3];
+  for (int i = 0; i < 3; ++i)
+  {
+    points.getTuple(i, pts[i]);
+  }
+
+  internal::Vector<TIn, 3> wc(component(wcoords, 0), component(wcoords, 1), component(wcoords, 2));
+
+  auto triangleNormal = internal::cross(pts[1] - pts[0], pts[2] - pts[0]);
+  for (IdComponent i = 0; i < 2; ++i)
+  {
+    auto& p0 = pts[0];
+    auto& p1 = pts[i + 1];
+    auto& p2 = pts[2 - i];
+    auto planeNormal = internal::cross(triangleNormal, p2 - p0);
+
+    component(pcoords, i) = static_cast<TOut>(internal::dot(wc - p0, planeNormal)) /
+                            static_cast<TOut>(internal::dot(p1 - p0, planeNormal));
+  }
+
+  return ErrorCode::SUCCESS;
+}
+
+} //namespace lcl
+
+#endif //lcl_Triangle_h

vtkm/thirdparty/lcl/vtkmlcl/lcl/Vertex.h

@@ -0,0 +1,131 @@
+//============================================================================
+//  Copyright (c) Kitware, Inc.
+//  All rights reserved.
+//  See LICENSE.md 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.
+//============================================================================
+#ifndef lcl_Vertex_h
+#define lcl_Vertex_h
+
+#include <lcl/ErrorCode.h>
+#include <lcl/Shapes.h>
+
+namespace lcl
+{
+
+class Vertex : public Cell
+{
+public:
+  constexpr LCL_EXEC Vertex() : Cell(ShapeId::VERTEX, 1) {}
+  constexpr LCL_EXEC explicit Vertex(const Cell& cell) : Cell(cell) {}
+};
+
+LCL_EXEC inline lcl::ErrorCode validate(Vertex tag) noexcept
+{
+  if (tag.shape() != ShapeId::VERTEX)
+  {
+    return ErrorCode::WRONG_SHAPE_ID_FOR_TAG_TYPE;
+  }
+  if (tag.numberOfPoints() != 1)
+  {
+    return ErrorCode::INVALID_NUMBER_OF_POINTS;
+  }
+
+  return ErrorCode::SUCCESS;
+}
+
+template<typename CoordType>
+LCL_EXEC inline lcl::ErrorCode parametricCenter(Vertex, CoordType&&) noexcept
+{
+  return ErrorCode::SUCCESS;
+}
+
+template<typename CoordType>
+LCL_EXEC inline lcl::ErrorCode parametricPoint(Vertex, IdComponent pointId, CoordType&&) noexcept
+{
+  if (pointId == 0)
+  {
+    return ErrorCode::SUCCESS;
+  }
+  else
+  {
+    return ErrorCode::INVALID_POINT_ID;
+  }
+}
+
+template<typename CoordType>
+LCL_EXEC inline ComponentType<CoordType> parametricDistance(Vertex, const CoordType&) noexcept
+{
+  return ComponentType<CoordType>{1};
+}
+
+template<typename CoordType>
+LCL_EXEC inline bool cellInside(Vertex, const CoordType&) noexcept
+{
+  return false;
+}
+
+template <typename Values, typename CoordType, typename Result>
+LCL_EXEC inline lcl::ErrorCode interpolate(
+  Vertex,
+  const Values& values,
+  const CoordType&,
+  Result&& result) noexcept
+{
+  using T = ComponentType<Result>;
+  for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
+  {
+    component(result, c) = static_cast<T>(values.getValue(0, c));
+  }
+
+  return ErrorCode::SUCCESS;
+}
+
+template <typename Points, typename Values, typename CoordType, typename Result>
+LCL_EXEC inline lcl::ErrorCode derivative(
+  Vertex,
+  const Points&,
+  const Values& values,
+  const CoordType&,
+  Result&& dx,
+  Result&& dy,
+  Result&& dz) noexcept
+{
+  using T = ComponentType<Result>;
+  for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
+  {
+    component(dx, c) = component(dy, c) = component(dz, c) = T{0};
+  }
+  return ErrorCode::SUCCESS;
+}
+
+template <typename Points, typename PCoordType, typename WCoordType>
+LCL_EXEC inline  lcl::ErrorCode parametricToWorld(
+  Vertex,
+  const Points& points,
+  const PCoordType&,
+  WCoordType&& wcoords) noexcept
+{
+  using T = ComponentType<WCoordType>;
+  component(wcoords, 0) = static_cast<T>(points.getValue(0, 0));
+  component(wcoords, 1) = static_cast<T>(points.getValue(0, 1));
+  component(wcoords, 2) = static_cast<T>(points.getValue(0, 2));
+  return ErrorCode::SUCCESS;
+}
+
+template <typename Points, typename WCoordType, typename PCoordType>
+LCL_EXEC inline lcl::ErrorCode worldToParametric(
+  Vertex,
+  const Points&,
+  const WCoordType&,
+  PCoordType&&) noexcept
+{
+  return ErrorCode::SUCCESS;
+}
+
+} //namespace lcl
+
+#endif //lcl_Vertex_h

vtkm/thirdparty/lcl/vtkmlcl/lcl/Voxel.h

@@ -0,0 +1,136 @@
+//============================================================================
+//  Copyright (c) Kitware, Inc.
+//  All rights reserved.
+//  See LICENSE.md 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.
+//============================================================================
+#ifndef lcl_Voxel_h
+#define lcl_Voxel_h
+
+#include <lcl/ErrorCode.h>
+#include <lcl/Hexahedron.h>
+
+#include <lcl/internal/Common.h>
+
+namespace lcl
+{
+
+class Voxel : public Hexahedron
+{
+public:
+  constexpr LCL_EXEC Voxel() : Hexahedron(Cell(ShapeId::VOXEL, 8)) {}
+  constexpr LCL_EXEC explicit Voxel(const Cell& cell) : Hexahedron(cell) {}
+};
+
+namespace internal
+{
+
+template <typename Points, typename T>
+LCL_EXEC inline int getVoxelSpacing(const Points& points, T spacing[3])
+{
+  int zeros = 0;
+  for (int i = 0; i < 3; ++i)
+  {
+    spacing[i] = static_cast<T>(points.getValue(6, i) - points.getValue(0, i));
+    if (spacing[i] == T(0))
+    {
+      zeros |= 1 << i;
+    }
+  }
+  return zeros;
+}
+
+} // internal
+
+template <typename Points, typename Values, typename CoordType, typename Result>
+LCL_EXEC inline lcl::ErrorCode derivative(
+  Voxel,
+  const Points& points,
+  const Values& values,
+  const CoordType& pcoords,
+  Result&& dx,
+  Result&& dy,
+  Result&& dz) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  using ProcessingType = internal::ClosestFloatType<typename Values::ValueType>;
+  using ResultCompType = ComponentType<Result>;
+
+  ProcessingType spacing[3];
+  if (internal::getVoxelSpacing(points, spacing) != 0)
+  {
+    return ErrorCode::DEGENERATE_CELL_DETECTED;
+  }
+
+  for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
+  {
+    ProcessingType dvdp[3];
+    internal::parametricDerivative(Hexahedron{}, values, c, pcoords, dvdp);
+
+    component(dx, c) = static_cast<ResultCompType>(dvdp[0] / spacing[0]);
+    component(dy, c) = static_cast<ResultCompType>(dvdp[1] / spacing[1]);
+    component(dz, c) = static_cast<ResultCompType>(dvdp[2] / spacing[2]);
+  }
+
+  return ErrorCode::SUCCESS;
+}
+
+template <typename Points, typename PCoordType, typename WCoordType>
+LCL_EXEC lcl::ErrorCode parametricToWorld(
+  Voxel,
+  const Points& points,
+  const PCoordType& pcoords,
+  WCoordType&& wcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(PCoordType);
+
+  using T = ComponentType<WCoordType>;
+
+  T spacing[3];
+  if (internal::getVoxelSpacing(points, spacing) != 0)
+  {
+    return ErrorCode::DEGENERATE_CELL_DETECTED;
+  }
+
+  component(wcoords, 0) = points.getValue(0, 0) +
+                          (spacing[0] * static_cast<T>(component(pcoords, 0)));
+  component(wcoords, 1) = points.getValue(0, 1) +
+                          (spacing[1] * static_cast<T>(component(pcoords, 1)));
+  component(wcoords, 2) = points.getValue(0, 2) +
+                          (spacing[2] * static_cast<T>(component(pcoords, 2)));
+  return ErrorCode::SUCCESS;
+}
+
+template <typename Points, typename WCoordType, typename PCoordType>
+LCL_EXEC lcl::ErrorCode worldToParametric(
+  Voxel,
+  const Points& points,
+  const WCoordType& wcoords,
+  PCoordType&& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(PCoordType);
+
+  using T = ComponentType<PCoordType>;
+
+  T spacing[3];
+  if (internal::getVoxelSpacing(points, spacing) != 0)
+  {
+    return ErrorCode::DEGENERATE_CELL_DETECTED;
+  }
+
+  component(pcoords, 0) = static_cast<T>(component(wcoords, 0) - points.getValue(0, 0)) /
+                          spacing[0];
+  component(pcoords, 1) = static_cast<T>(component(wcoords, 1) - points.getValue(0, 1)) /
+                          spacing[1];
+  component(pcoords, 2) = static_cast<T>(component(wcoords, 2) - points.getValue(0, 2)) /
+                          spacing[2];
+  return ErrorCode::SUCCESS;
+}
+
+} //namespace lcl
+
+#endif //lcl_Voxel_h

vtkm/thirdparty/lcl/vtkmlcl/lcl/Wedge.h

@@ -0,0 +1,230 @@
+//============================================================================
+//  Copyright (c) Kitware, Inc.
+//  All rights reserved.
+//  See LICENSE.md 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.
+//============================================================================
+#ifndef lcl_Wedge_h
+#define lcl_Wedge_h
+
+#include <lcl/ErrorCode.h>
+#include <lcl/Shapes.h>
+
+#include <lcl/internal/Common.h>
+
+namespace lcl
+{
+
+class Wedge : public Cell
+{
+public:
+  constexpr LCL_EXEC Wedge() : Cell(ShapeId::WEDGE, 6) {}
+  constexpr LCL_EXEC explicit Wedge(const Cell& cell) : Cell(cell) {}
+};
+
+LCL_EXEC inline lcl::ErrorCode validate(Wedge tag) noexcept
+{
+  if (tag.shape() != ShapeId::WEDGE)
+  {
+    return ErrorCode::WRONG_SHAPE_ID_FOR_TAG_TYPE;
+  }
+  if (tag.numberOfPoints() != 6)
+  {
+    return ErrorCode::INVALID_NUMBER_OF_POINTS;
+  }
+
+  return ErrorCode::SUCCESS;
+}
+
+template<typename CoordType>
+LCL_EXEC inline lcl::ErrorCode parametricCenter(Wedge, CoordType&& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  using T = ComponentType<CoordType>;
+  component(pcoords, 0) = T(1)/T(3);
+  component(pcoords, 1) = T(1)/T(3);
+  component(pcoords, 2) = 0.5f;
+  return ErrorCode::SUCCESS;
+}
+
+template<typename CoordType>
+LCL_EXEC inline lcl::ErrorCode parametricPoint(
+  Wedge, IdComponent pointId, CoordType&& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  switch (pointId)
+  {
+    case 0:
+      component(pcoords, 0) = 0.0f;
+      component(pcoords, 1) = 0.0f;
+      component(pcoords, 2) = 0.0f;
+      break;
+    case 1:
+      component(pcoords, 0) = 1.0f;
+      component(pcoords, 1) = 0.0f;
+      component(pcoords, 2) = 0.0f;
+      break;
+    case 2:
+      component(pcoords, 0) = 0.0f;
+      component(pcoords, 1) = 1.0f;
+      component(pcoords, 2) = 0.0f;
+      break;
+    case 3:
+      component(pcoords, 0) = 0.0f;
+      component(pcoords, 1) = 0.0f;
+      component(pcoords, 2) = 1.0f;
+      break;
+    case 4:
+      component(pcoords, 0) = 1.0f;
+      component(pcoords, 1) = 0.0f;
+      component(pcoords, 2) = 1.0f;
+      break;
+    case 5:
+      component(pcoords, 0) = 0.0f;
+      component(pcoords, 1) = 1.0f;
+      component(pcoords, 2) = 1.0f;
+      break;
+    default:
+      return ErrorCode::INVALID_POINT_ID;
+  }
+
+  return ErrorCode::SUCCESS;
+}
+
+template<typename CoordType>
+LCL_EXEC inline ComponentType<CoordType> parametricDistance(Wedge, const CoordType& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+  return internal::findParametricDistance(pcoords, 3);
+}
+
+template<typename CoordType>
+LCL_EXEC inline bool cellInside(Wedge, const CoordType& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  using T = ComponentType<CoordType>;
+
+  constexpr T eps = 0.001f;
+  return component(pcoords, 0) >= -eps &&
+         component(pcoords, 1) >= -eps &&
+         component(pcoords, 2) >= -eps &&
+         (component(pcoords, 0) + component(pcoords, 1)) <= (T{1} + eps) &&
+         component(pcoords, 2) <= (T{1} + eps);
+}
+
+template <typename Values, typename CoordType, typename Result>
+LCL_EXEC lcl::ErrorCode interpolate(
+  Wedge,
+  const Values& values,
+  const CoordType& pcoords,
+  Result&& result) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  using T = internal::ClosestFloatType<typename Values::ValueType>;
+
+  auto p0 = static_cast<T>(component(pcoords, 0));
+  auto p1 = static_cast<T>(component(pcoords, 1));
+  auto p2 = static_cast<T>(component(pcoords, 2));
+  auto sm = T{1} - p0 - p1;
+
+  for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
+  {
+    T t1 = static_cast<T>(values.getValue(0, c)) * sm +
+           static_cast<T>(values.getValue(1, c)) * p0 +
+           static_cast<T>(values.getValue(2, c)) * p1;
+    T t2 = static_cast<T>(values.getValue(3, c)) * sm +
+           static_cast<T>(values.getValue(4, c)) * p0 +
+           static_cast<T>(values.getValue(5, c)) * p1;
+    component(result, c) = static_cast<ComponentType<Result>>(internal::lerp(t1, t2, p2));
+  }
+
+  return ErrorCode::SUCCESS;
+}
+
+namespace internal
+{
+
+template <typename Values, typename CoordType, typename Result>
+LCL_EXEC inline void parametricDerivative(
+  Wedge, const Values& values, IdComponent comp, const CoordType& pcoords, Result&& result) noexcept
+{
+  using T = internal::ClosestFloatType<typename Values::ValueType>;
+  auto p0 = static_cast<T>(component(pcoords, 0));
+  auto p1 = static_cast<T>(component(pcoords, 1));
+  auto p2 = static_cast<T>(component(pcoords, 2));
+  auto rm = T{1} - p2;
+  auto sm = T{1} - p0 - p1;
+
+  T dr = (static_cast<T>(values.getValue(0, comp)) * -rm) +
+         (static_cast<T>(values.getValue(1, comp)) *  rm) +
+         (static_cast<T>(values.getValue(3, comp)) * -p2) +
+         (static_cast<T>(values.getValue(4, comp)) *  p2);
+
+  T ds = (static_cast<T>(values.getValue(0, comp)) * -rm) +
+         (static_cast<T>(values.getValue(2, comp)) *  rm) +
+         (static_cast<T>(values.getValue(3, comp)) * -p2) +
+         (static_cast<T>(values.getValue(5, comp)) *  p2);
+
+  T dt = (static_cast<T>(values.getValue(0, comp)) * -sm) +
+         (static_cast<T>(values.getValue(1, comp)) * -p0) +
+         (static_cast<T>(values.getValue(2, comp)) * -p1) +
+         (static_cast<T>(values.getValue(3, comp)) *  sm) +
+         (static_cast<T>(values.getValue(4, comp)) *  p0) +
+         (static_cast<T>(values.getValue(5, comp)) *  p1);
+
+  component(result, 0) = static_cast<ComponentType<Result>>(dr);
+  component(result, 1) = static_cast<ComponentType<Result>>(ds);
+  component(result, 2) = static_cast<ComponentType<Result>>(dt);
+}
+
+} // internal
+
+template <typename Points, typename Values, typename CoordType, typename Result>
+LCL_EXEC inline lcl::ErrorCode derivative(
+  Wedge,
+  const Points& points,
+  const Values& values,
+  const CoordType& pcoords,
+  Result&& dx,
+  Result&& dy,
+  Result&& dz) noexcept
+{
+  return internal::derivative3D(Wedge{},
+                                points,
+                                values,
+                                pcoords,
+                                std::forward<Result>(dx),
+                                std::forward<Result>(dy),
+                                std::forward<Result>(dz));
+}
+
+template <typename Points, typename PCoordType, typename WCoordType>
+LCL_EXEC inline lcl::ErrorCode parametricToWorld(
+  Wedge,
+  const Points& points,
+  const PCoordType& pcoords,
+  WCoordType&& wcoords) noexcept
+{
+  return interpolate(Wedge{}, points, pcoords, std::forward<WCoordType>(wcoords));
+}
+
+template <typename Points, typename WCoordType, typename PCoordType>
+LCL_EXEC inline lcl::ErrorCode worldToParametric(
+  Wedge,
+  const Points& points,
+  const WCoordType& wcoords,
+  PCoordType&& pcoords) noexcept
+{
+  return internal::worldToParametric3D(Wedge{}, points, wcoords, std::forward<PCoordType>(pcoords));
+}
+
+} // lcl
+
+#endif // lcl_Wedge_h

vtkm/thirdparty/lcl/vtkmlcl/lcl/internal/Common.h

@@ -0,0 +1,316 @@
+//============================================================================
+//  Copyright (c) Kitware, Inc.
+//  All rights reserved.
+//  See LICENSE.md 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.
+//============================================================================
+#ifndef lcl_internal_Common_h
+#define lcl_internal_Common_h
+
+#include <lcl/FieldAccessor.h>
+#include <lcl/internal/Math.h>
+
+#include <cstdint>
+#include <type_traits>
+
+#define LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(PCoordType)                   \
+  static_assert(std::is_floating_point<ComponentType<PCoordType>>::value,      \
+                "parametric coordinates should be of floating point type");
+
+namespace lcl
+{
+namespace internal
+{
+
+///=========================================================================
+template<typename CoordType>
+LCL_EXEC inline ComponentType<CoordType> findParametricDistance(
+  const CoordType& pvals, IdComponent numVals) noexcept
+{
+  using T = ComponentType<CoordType>;
+
+  T pDistMax = 0.0f;
+  for (IdComponent i = 0; i < numVals; ++i)
+  {
+    ComponentType<CoordType> pDist = 0.0f;
+    if (component(pvals, i) < T(0))
+    {
+      pDist = -1.0f * component(pvals, i);
+    }
+    else if (component(pvals, i) > T(1))
+    {
+      pDist = component(pvals, i) - 1.0f;
+    }
+    if (pDist > pDistMax)
+    {
+      pDistMax = pDist;
+    }
+  }
+
+  return pDistMax;
+}
+
+///=========================================================================
+/// Forward declaration
+#define FORWARD_DECLAR_PARAMETRIC_DERIVATIVE(tag)                                                  \
+template <typename Values, typename CoordType, typename Result>                                    \
+LCL_EXEC inline void parametricDerivative(                                                        \
+  lcl::tag, const Values& values, IdComponent comp, const CoordType&, Result&& result) noexcept
+
+FORWARD_DECLAR_PARAMETRIC_DERIVATIVE(Triangle);
+FORWARD_DECLAR_PARAMETRIC_DERIVATIVE(Quad);
+FORWARD_DECLAR_PARAMETRIC_DERIVATIVE(Tetra);
+FORWARD_DECLAR_PARAMETRIC_DERIVATIVE(Hexahedron);
+FORWARD_DECLAR_PARAMETRIC_DERIVATIVE(Wedge);
+FORWARD_DECLAR_PARAMETRIC_DERIVATIVE(Pyramid);
+
+#undef FORWARD_DECLAR_PARAMETRIC_DERIVATIVE
+
+///=========================================================================
+template <typename T>
+class Space2D
+{
+public:
+  explicit LCL_EXEC Space2D(const Vector<T, 3>& origin, const Vector<T, 3>& p1, const Vector<T, 3>& p2) noexcept
+  {
+    this->Origin = origin;
+    this->XAxis = p1 - origin;
+    auto normal = internal::cross(this->XAxis, p2 - origin);
+    this->YAxis = internal::cross(normal, this->XAxis);
+
+    internal::normalize(this->XAxis);
+    internal::normalize(this->YAxis);
+  }
+
+  LCL_EXEC Vector<T, 2> to2DPoint(Vector<T, 3> pt) const noexcept
+  {
+    pt -= this->Origin;
+    return Vector<T, 2>{ internal::dot(pt, this->XAxis), internal::dot(pt, this->YAxis) };
+  }
+
+  LCL_EXEC Vector<T, 3> to3DVec(const Vector<T, 2>& vec) const noexcept
+  {
+    return (this->XAxis * vec[0]) + (this->YAxis * vec[1]);
+  }
+
+private:
+  Vector<T, 3> Origin;
+  Vector<T, 3> XAxis, YAxis;
+};
+
+template <typename CellTag, typename Points, typename PCoords, typename T>
+LCL_EXEC inline void jacobian2D(
+  CellTag tag, const Points& points, const PCoords& pcoords, Matrix<T, 2, 2>& jacobian) noexcept
+{
+  T pd[2];
+  parametricDerivative(tag, points, 0, pcoords, pd);
+  jacobian(0, 0) = pd[0];
+  jacobian(0, 1) = pd[1];
+  parametricDerivative(tag, points, 1, pcoords, pd);
+  jacobian(1, 0) = pd[0];
+  jacobian(1, 1) = pd[1];
+}
+
+template <typename CellTag, typename Points, typename Values, typename CoordType, typename Result>
+LCL_EXEC inline lcl::ErrorCode derivative2D(
+  CellTag tag,
+  const Points& points,
+  const Values& values,
+  const CoordType& pcoords,
+  Result&& dx,
+  Result&& dy,
+  Result&& dz) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  using ProcessingType = internal::ClosestFloatType<typename Values::ValueType>;
+  using ResultCompType = ComponentType<Result>;
+
+  constexpr IdComponent numPoints = CellTag{}.numberOfPoints();
+
+  Vector<ProcessingType, 3> pts[numPoints];
+  for (int i = 0; i < numPoints; ++i)
+  {
+    points.getTuple(i, pts[i]);
+  }
+
+  // 2-D coordinate system on the cell's plane
+  Space2D<ProcessingType> planeSpace(pts[0], pts[1], pts[numPoints - 1]);
+  Vector<ProcessingType, 2> pts2d[numPoints];
+  for (int i = 0; i < numPoints; ++i)
+  {
+    pts2d[i] = planeSpace.to2DPoint(pts[i]);
+  }
+
+  Matrix<ProcessingType, 2, 2> jacobian;
+  jacobian2D(tag, makeFieldAccessorNestedSOA(pts2d, 2), pcoords, jacobian);
+  Matrix<ProcessingType, 2, 2> invJacobian;
+  LCL_RETURN_ON_ERROR(matrixInverse(jacobian, invJacobian))
+
+  for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
+  {
+    Vector<ProcessingType, 2> dvdp;
+    parametricDerivative(tag, values, c, pcoords, dvdp);
+    auto d2D = matrixMultiply(dvdp, invJacobian);
+    auto d3D = planeSpace.to3DVec(d2D);
+
+    component(dx, c) = static_cast<ResultCompType>(d3D[0]);
+    component(dy, c) = static_cast<ResultCompType>(d3D[1]);
+    component(dz, c) = static_cast<ResultCompType>(d3D[2]);
+  }
+
+  return ErrorCode::SUCCESS;
+}
+
+template <typename CellTag, typename Points, typename WCoordType, typename PCoordType>
+LCL_EXEC inline lcl::ErrorCode worldToParametric2D(
+  CellTag tag,
+  const Points& points,
+  const WCoordType& wcoords,
+  PCoordType&& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(PCoordType);
+
+  using TIn = typename Points::ValueType;
+  using TOut = ComponentType<PCoordType>;
+
+  constexpr IdComponent numPoints = CellTag{}.numberOfPoints();
+
+  Vector<TIn, 3> pts[numPoints];
+  for (int i = 0; i < numPoints; ++i)
+  {
+    points.getTuple(i, pts[i]);
+  }
+
+  // 2-D coordinate system on the cell's plane
+  Space2D<TIn> planeSpace(pts[0], pts[1], pts[numPoints - 1]);
+  Vector<TIn, 2> pts2d[numPoints];
+  for (int i = 0; i < numPoints; ++i)
+  {
+    pts2d[i] = planeSpace.to2DPoint(pts[i]);
+  }
+
+  auto jacobianEvaluator = [&pts2d](const Vector<TOut, 2>& pc, Matrix<TIn, 2, 2>& jacobian) {
+    jacobian2D(CellTag{}, makeFieldAccessorNestedSOA(pts2d, 2), pc, jacobian);
+    return ErrorCode::SUCCESS;
+  };
+
+  auto functionEvaluator = [&points, &planeSpace](const Vector<TOut, 2>& pc, Vector<TIn, 2>& wc) {
+    Vector<TIn, 3> wc3(0);
+    LCL_RETURN_ON_ERROR(parametricToWorld(CellTag{}, points, pc, wc3))
+    wc = planeSpace.to2DPoint(wc3);
+
+    return ErrorCode::SUCCESS;
+  };
+
+  Vector<TIn, 3> wcVec{component(wcoords, 0), component(wcoords, 1), component(wcoords, 2)};
+  Vector<TOut, 2> pcVec;
+  LCL_RETURN_ON_ERROR(parametricCenter(tag, pcVec))
+  auto status = newtonsMethod(
+    jacobianEvaluator, functionEvaluator, planeSpace.to2DPoint(wcVec), pcVec);
+
+  if (status == ErrorCode::SUCCESS || status == ErrorCode::SOLUTION_DID_NOT_CONVERGE)
+  {
+    component(pcoords, 0) = pcVec[0];
+    component(pcoords, 1) = pcVec[1];
+  }
+
+  return status;
+}
+
+///=========================================================================
+template <typename CellTag, typename Points, typename PCoords, typename T>
+LCL_EXEC inline void jacobian3D(
+  CellTag tag, const Points& points, const PCoords& pcoords, Matrix<T, 3, 3>& jacobian) noexcept
+{
+  T pd[3];
+  parametricDerivative(tag, points, 0, pcoords, pd);
+  jacobian(0, 0) = pd[0];
+  jacobian(0, 1) = pd[1];
+  jacobian(0, 2) = pd[2];
+  parametricDerivative(tag, points, 1, pcoords, pd);
+  jacobian(1, 0) = pd[0];
+  jacobian(1, 1) = pd[1];
+  jacobian(1, 2) = pd[2];
+  parametricDerivative(tag, points, 2, pcoords, pd);
+  jacobian(2, 0) = pd[0];
+  jacobian(2, 1) = pd[1];
+  jacobian(2, 2) = pd[2];
+}
+
+template <typename CellTag, typename Points, typename Values, typename CoordType, typename Result>
+LCL_EXEC inline lcl::ErrorCode derivative3D(
+  CellTag tag,
+  const Points& points,
+  const Values& values,
+  const CoordType& pcoords,
+  Result&& dx,
+  Result&& dy,
+  Result&& dz) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  using ProcessingType = internal::ClosestFloatType<typename Values::ValueType>;
+  using ResultCompType = ComponentType<Result>;
+
+  Matrix<ProcessingType, 3, 3> jacobian;
+  jacobian3D(tag, points, pcoords, jacobian);
+  Matrix<ProcessingType, 3, 3> invJacobian;
+  LCL_RETURN_ON_ERROR(matrixInverse(jacobian, invJacobian))
+
+  for (IdComponent c = 0; c < values.getNumberOfComponents(); ++c)
+  {
+    Vector<ProcessingType, 3> dvdp;
+    parametricDerivative(tag, values, c, pcoords, dvdp);
+    auto deriv = matrixMultiply(dvdp, invJacobian);
+    component(dx, c) = static_cast<ResultCompType>(deriv[0]);
+    component(dy, c) = static_cast<ResultCompType>(deriv[1]);
+    component(dz, c) = static_cast<ResultCompType>(deriv[2]);
+  }
+
+  return ErrorCode::SUCCESS;
+}
+
+template <typename CellTag, typename Points, typename WCoordType, typename PCoordType>
+LCL_EXEC inline lcl::ErrorCode worldToParametric3D(
+  CellTag tag,
+  const Points& points,
+  const WCoordType& wcoords,
+  PCoordType&& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(PCoordType);
+
+  using TIn = typename Points::ValueType;
+  using TOut = ComponentType<PCoordType>;
+
+  auto jacobianEvaluator = [tag, &points](const Vector<TOut, 3>& pc, Matrix<TIn, 3, 3>& jacobian) {
+    jacobian3D(tag, points, pc, jacobian);
+    return ErrorCode::SUCCESS;
+  };
+
+  auto functionEvaluator = [tag, &points](const Vector<TOut, 3>& pc, Vector<TIn, 3>& wc) {
+    return parametricToWorld(tag, points, pc, wc);
+  };
+
+  internal::Vector<TIn, 3> wcVec{component(wcoords, 0), component(wcoords, 1), component(wcoords, 2)};
+  internal::Vector<TOut, 3> pcVec;
+  LCL_RETURN_ON_ERROR(parametricCenter(tag, pcVec))
+  auto status = newtonsMethod(jacobianEvaluator, functionEvaluator, wcVec, pcVec);
+
+  if (status == ErrorCode::SUCCESS || status == ErrorCode::SOLUTION_DID_NOT_CONVERGE)
+  {
+    component(pcoords, 0) = pcVec[0];
+    component(pcoords, 1) = pcVec[1];
+    component(pcoords, 2) = pcVec[2];
+  }
+
+  return status;
+}
+
+}
+} // lcl::internal
+
+#endif //lcl_internal_Common_h

vtkm/thirdparty/lcl/vtkmlcl/lcl/internal/Config.h

@@ -0,0 +1,38 @@
+//============================================================================
+//  Copyright (c) Kitware, Inc.
+//  All rights reserved.
+//  See LICENSE.md 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.
+//============================================================================
+#ifndef lcl_internal_Config_h
+#define lcl_internal_Config_h
+
+#include <cstdint>
+#include <type_traits>
+
+#ifdef __CUDACC__
+# define LCL_EXEC __device__ __host__
+#else
+# define LCL_EXEC
+#endif
+
+namespace lcl
+{
+
+namespace internal
+{
+template <typename T>
+using ClosestFloatType =
+  typename std::enable_if<std::is_arithmetic<T>::value,
+                          typename std::conditional<sizeof(T) <= 4, float, double>::type>::type;
+}
+
+using IdShape = std::int8_t;
+using IdComponent = std::int32_t;
+
+} // namespace lcl
+
+#endif // lcl_internal_Config_h

vtkm/thirdparty/lcl/vtkmlcl/lcl/internal/Math.h

@@ -0,0 +1,712 @@
+//============================================================================
+//  Copyright (c) Kitware, Inc.
+//  All rights reserved.
+//  See LICENSE.md 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.
+//============================================================================
+#ifndef lcl_internal_Math_h
+#define lcl_internal_Math_h
+
+#include <cmath>
+
+#ifdef __CUDA_ARCH__
+# define LCL_MATH_CALL(func, ...) func(__VA_ARGS__)
+#else
+# define LCL_MATH_CALL(func, ...) std::func(__VA_ARGS__)
+#endif
+
+namespace lcl
+{
+
+namespace internal
+{
+///=========================================================================
+/// Basic vector class and functions for internal use
+
+template <typename T, IdComponent Dim>
+class Vector
+{
+public:
+#if defined(_MSC_VER) && (_MSC_VER == 1900)
+  // workaround vs2015 bug producing the following error:
+  // error C2476: ‘constexpr’ constructor does not initialize all members
+  LCL_EXEC Vector() noexcept {};
+#else
+  LCL_EXEC
+  constexpr Vector() noexcept {};
+#endif
+
+  LCL_EXEC
+  explicit Vector(const T& val) noexcept
+  {
+    for (auto& c : Data)
+    {
+      c = val;
+    }
+  }
+
+  template <typename... Ts>
+  LCL_EXEC
+  constexpr explicit Vector(const T& c1, const Ts&... cs) noexcept
+    : Data{c1, cs...}
+  {
+    static_assert(sizeof...(Ts) == (Dim - 1), "Invalid number of components passed");
+  }
+
+  static constexpr IdComponent getNumberOfComponents() noexcept
+  {
+    return Dim;
+  }
+
+  LCL_EXEC
+  T& operator[](IdComponent c) noexcept
+  {
+    return this->Data[c];
+  }
+
+  LCL_EXEC
+  constexpr const T& operator[](IdComponent c) const noexcept
+  {
+    return this->Data[c];
+  }
+
+  LCL_EXEC
+  T* data() noexcept
+  {
+    return this->Data;
+  }
+
+  LCL_EXEC
+  constexpr const T* data() const noexcept
+  {
+    return this->Data;
+  }
+
+private:
+  T Data[Dim];
+};
+
+//---------------------------------------------------------------------------
+template <typename T, IdComponent Dim>
+LCL_EXEC
+Vector<T, Dim>& operator+=(Vector<T, Dim>& v1, const Vector<T, Dim>& v2) noexcept
+{
+  for (IdComponent i = 0; i < Dim; ++i)
+  {
+    v1[i] += v2[i];
+  }
+  return v1;
+}
+
+template <typename T, IdComponent Dim>
+LCL_EXEC
+Vector<T, Dim>& operator-=(Vector<T, Dim>& v1, const Vector<T, Dim>& v2) noexcept
+{
+  for (IdComponent i = 0; i < Dim; ++i)
+  {
+    v1[i] -= v2[i];
+  }
+  return v1;
+}
+
+template <typename T, IdComponent Dim>
+LCL_EXEC
+Vector<T, Dim>& operator*=(Vector<T, Dim>& v1, const Vector<T, Dim>& v2) noexcept
+{
+  for (IdComponent i = 0; i < Dim; ++i)
+  {
+    v1[i] *= v2[i];
+  }
+  return v1;
+}
+
+template <typename T, IdComponent Dim>
+LCL_EXEC
+Vector<T, Dim>& operator*=(Vector<T, Dim>& v, const T& s) noexcept
+{
+  for (IdComponent i = 0; i < Dim; ++i)
+  {
+    v[i] *= s;
+  }
+  return v;
+}
+
+template <typename T, IdComponent Dim>
+LCL_EXEC
+Vector<T, Dim>& operator/=(Vector<T, Dim>& v1, const Vector<T, Dim>& v2) noexcept
+{
+  for (IdComponent i = 0; i < Dim; ++i)
+  {
+    v1[i] /= v2[i];
+  }
+  return v1;
+}
+
+template <typename T, IdComponent Dim>
+LCL_EXEC
+Vector<T, Dim>& operator/=(Vector<T, Dim>& v, const T& s) noexcept
+{
+  for (IdComponent i = 0; i < Dim; ++i)
+  {
+    v[i] /= s;
+  }
+  return v;
+}
+
+template <typename T, IdComponent Dim>
+LCL_EXEC
+Vector<T, Dim> operator+(Vector<T, Dim> v1, const Vector<T, Dim>& v2) noexcept
+{
+  return v1 += v2;
+}
+
+template <typename T, IdComponent Dim>
+LCL_EXEC
+Vector<T, Dim> operator-(Vector<T, Dim> v1, const Vector<T, Dim>& v2) noexcept
+{
+  return v1 -= v2;
+}
+
+template <typename T, IdComponent Dim>
+LCL_EXEC
+Vector<T, Dim> operator*(Vector<T, Dim> v1, const Vector<T, Dim>& v2) noexcept
+{
+  return v1 *= v2;
+}
+
+template <typename T, IdComponent Dim>
+LCL_EXEC
+Vector<T, Dim> operator*(Vector<T, Dim> v, const T& s) noexcept
+{
+  return v *= s;
+}
+
+template <typename T, IdComponent Dim>
+LCL_EXEC
+Vector<T, Dim> operator*(const T& s, Vector<T, Dim> v) noexcept
+{
+  return v *= s;
+}
+
+template <typename T, IdComponent Dim>
+LCL_EXEC
+Vector<T, Dim> operator/(Vector<T, Dim> v1, const Vector<T, Dim>& v2) noexcept
+{
+  return v1 /= v2;
+}
+
+template <typename T, IdComponent Dim>
+LCL_EXEC
+Vector<T, Dim> operator/(Vector<T, Dim> v, const T& s) noexcept
+{
+  return v /= s;
+}
+
+//---------------------------------------------------------------------------
+template <typename T, IdComponent Dim>
+LCL_EXEC
+T dot(const Vector<T, Dim>& v1, const Vector<T, Dim>& v2) noexcept
+{
+  T result{};
+  for (IdComponent i = 0; i < Dim; ++i)
+  {
+    result += v1[i] * v2[i];
+  }
+  return result;
+}
+
+template <typename T>
+LCL_EXEC
+Vector<T, 3> cross(const Vector<T, 3>& v1, const Vector<T, 3>& v2) noexcept
+{
+  return Vector<T, 3>(v1[1] * v2[2] - v1[2] * v2[1],
+                      v1[2] * v2[0] - v1[0] * v2[2],
+                      v1[0] * v2[1] - v1[1] * v2[0]);
+}
+
+template <typename T, IdComponent Dim>
+LCL_EXEC
+T magnitude(const Vector<T, Dim>& v)
+{
+  return LCL_MATH_CALL(sqrt, (dot(v, v)));
+}
+
+template <typename T, IdComponent Dim>
+LCL_EXEC
+void normalize(Vector<T, Dim>& v)
+{
+  v /= magnitude(v);
+}
+
+template <typename T, IdComponent Dim>
+LCL_EXEC
+Vector<T, Dim> normal(Vector<T, Dim> v)
+{
+  normalize(v);
+  return v;
+}
+
+///=========================================================================
+/// Basic matrix class and functions for internal use
+
+template <typename T, IdComponent NumberOfRows, IdComponent NumberOfColumns>
+class Matrix
+{
+public:
+#if defined(_MSC_VER) && (_MSC_VER == 1900)
+  // workaround vs2015 bug producing the following error:
+  // error C2476: ‘constexpr’ constructor does not initialize all members
+  LCL_EXEC Matrix() noexcept {};
+#else
+  LCL_EXEC
+  constexpr Matrix() noexcept {};
+#endif
+
+  LCL_EXEC
+  T& operator()(IdComponent r, IdComponent c) noexcept
+  {
+    return this->Columns[c][r];
+  }
+
+  LCL_EXEC
+  constexpr const T& operator()(IdComponent r, IdComponent c) const noexcept
+  {
+    return this->Columns[c][r];
+  }
+
+  LCL_EXEC
+  static constexpr IdComponent getNumberOfRows() noexcept
+  {
+    return NumberOfRows;
+  }
+
+  LCL_EXEC
+  static constexpr IdComponent getNumberOfColumns() noexcept
+  {
+    return NumberOfColumns;
+  }
+
+  LCL_EXEC
+  constexpr const Vector<T, NumberOfRows>& getColumn(IdComponent c) const
+  {
+    return this->Columns[c];
+  }
+
+  LCL_EXEC
+  Vector<T, NumberOfColumns> getRow(IdComponent r) const
+  {
+    Vector<T, NumberOfColumns> row;
+    for (IdComponent i = 0; i < NumberOfColumns; ++i)
+    {
+      row[i] = this->Columns[i][r];
+    }
+    return row;
+  }
+
+  LCL_EXEC
+  Matrix& operator+=(const Matrix& m) noexcept
+  {
+    for (IdComponent i = 0; i < NumberOfColumns; ++i)
+    {
+      this->Columns[i] += m.Columns[i];
+    }
+    return *this;
+  }
+
+  LCL_EXEC
+  Matrix& operator-=(const Matrix& m) noexcept
+  {
+    for (IdComponent i = 0; i < NumberOfColumns; ++i)
+    {
+      this->Columns[i] -= m.Columns[i];
+    }
+    return *this;
+  }
+
+  LCL_EXEC
+  Matrix& operator*=(const T& s) noexcept
+  {
+    for (auto& c : this->Columns)
+    {
+      c *= s;
+    }
+    return *this;
+  }
+
+  LCL_EXEC
+  Matrix& operator/=(const T& s) noexcept
+  {
+    for (auto& c : this->Columns)
+    {
+      c /= s;
+    }
+    return *this;
+  }
+
+private:
+  Vector<T, NumberOfRows> Columns[NumberOfColumns];
+};
+
+//---------------------------------------------------------------------------
+template <typename T, IdComponent NRows, IdComponent NCols>
+LCL_EXEC
+Matrix<T, NRows, NCols> operator+(Matrix<T, NRows, NCols> m1,
+                                  const Matrix<T, NRows, NCols>& m2) noexcept
+{
+  return m1 += m2;
+}
+
+template <typename T, IdComponent NRows, IdComponent NCols>
+LCL_EXEC
+Matrix<T, NRows, NCols> operator-(Matrix<T, NRows, NCols> m1,
+                                  const Matrix<T, NRows, NCols>& m2) noexcept
+{
+  return m1 -= m2;
+}
+
+template <typename T, IdComponent NRows, IdComponent NCols>
+LCL_EXEC
+Matrix<T, NRows, NCols> operator*(Matrix<T, NRows, NCols> m, const T& s) noexcept
+{
+  return m *= s;
+}
+
+template <typename T, IdComponent NRows, IdComponent NCols>
+LCL_EXEC
+Matrix<T, NRows, NCols> operator/(Matrix<T, NRows, NCols> m, const T& s) noexcept
+{
+  return m /= s;
+}
+
+//---------------------------------------------------------------------------
+template <typename T,
+          IdComponent NumRow,
+          IdComponent NumCol,
+          IdComponent NumInternal>
+LCL_EXEC
+Matrix<T, NumRow, NumCol> matrixMultiply(
+  const Matrix<T, NumRow, NumInternal>& leftFactor,
+  const Matrix<T, NumInternal, NumCol>& rightFactor) noexcept
+{
+  Matrix<T, NumRow, NumCol> result;
+  for (IdComponent rowIndex = 0; rowIndex < NumRow; ++rowIndex)
+  {
+    for (IdComponent colIndex = 0; colIndex < NumCol; ++colIndex)
+    {
+      T sum = T(leftFactor(rowIndex, 0) * rightFactor(0, colIndex));
+      for (IdComponent internalIndex = 1; internalIndex < NumInternal; internalIndex++)
+      {
+        sum = T(sum + (leftFactor(rowIndex, internalIndex) * rightFactor(internalIndex, colIndex)));
+      }
+      result(rowIndex, colIndex) = sum;
+    }
+  }
+  return result;
+}
+
+template <typename T, IdComponent NumRow, IdComponent NumCol>
+LCL_EXEC
+Vector<T, NumRow> matrixMultiply(
+  const Matrix<T, NumRow, NumCol>& leftFactor,
+  const Vector<T, NumCol>& rightFactor) noexcept
+{
+  Vector<T, NumRow> result;
+  for (IdComponent rowIndex = 0; rowIndex < NumRow; ++rowIndex)
+  {
+    result[rowIndex] = dot(leftFactor.getRow(rowIndex), rightFactor);
+  }
+  return result;
+}
+
+template <typename T, IdComponent NumRow, IdComponent NumCol>
+LCL_EXEC
+Vector<T, NumCol> matrixMultiply(
+  const Vector<T, NumRow>& leftFactor,
+  const Matrix<T, NumRow, NumCol>& rightFactor) noexcept
+{
+  Vector<T, NumCol> result;
+  for (IdComponent colIndex = 0; colIndex < NumCol; ++colIndex)
+  {
+    result[colIndex] = dot(leftFactor, rightFactor.getColumn(colIndex));
+  }
+  return result;
+}
+
+//---------------------------------------------------------------------------
+// The following code is copied from VTK-m
+
+namespace detail
+{
+
+template <typename T, IdComponent Size>
+LCL_EXEC
+lcl::ErrorCode matrixLUPFactorFindPivot(Matrix<T, Size, Size>& m,
+                                         Vector<IdComponent, Size>& permutation,
+                                         IdComponent topCornerIndex,
+                                         T& inversionParity) noexcept
+{
+  constexpr T epsilon = std::is_floating_point<T>::value ?
+                        (std::is_same<T, float>::value ? T(1e-5f) : T(1e-9)) :
+                        T{};
+
+  IdComponent maxRowIndex = topCornerIndex;
+  T maxValue = LCL_MATH_CALL(abs, (m(maxRowIndex, topCornerIndex)));
+  for (IdComponent rowIndex = topCornerIndex + 1; rowIndex < Size; rowIndex++)
+  {
+    T compareValue = LCL_MATH_CALL(abs, (m(rowIndex, topCornerIndex)));
+    if (maxValue < compareValue)
+    {
+      maxValue = compareValue;
+      maxRowIndex = rowIndex;
+    }
+  }
+
+  if (maxValue < epsilon)
+  {
+    return lcl::ErrorCode::MATRIX_LUP_FACTORIZATION_FAILED;
+  }
+
+  if (maxRowIndex != topCornerIndex)
+  {
+    // Swap rows in matrix.
+    for (IdComponent i = 0; i < Size; ++i)
+    {
+      auto temp = m(topCornerIndex, i);
+      m(topCornerIndex, i) = m(maxRowIndex, i);
+      m(maxRowIndex, i) = temp;
+    }
+
+    // Record change in permutation matrix.
+    IdComponent maxOriginalRowIndex = permutation[maxRowIndex];
+    permutation[maxRowIndex] = permutation[topCornerIndex];
+    permutation[topCornerIndex] = maxOriginalRowIndex;
+
+    // Keep track of inversion parity.
+    inversionParity = -inversionParity;
+  }
+
+  return lcl::ErrorCode::SUCCESS;
+}
+
+// Used with MatrixLUPFactor
+template <typename T, IdComponent Size>
+LCL_EXEC
+void matrixLUPFactorFindUpperTriangleElements(
+  Matrix<T, Size, Size>& m, IdComponent topCornerIndex) noexcept
+{
+  // Compute values for upper triangle on row topCornerIndex
+  for (IdComponent colIndex = topCornerIndex + 1; colIndex < Size; colIndex++)
+  {
+    m(topCornerIndex, colIndex) /= m(topCornerIndex, topCornerIndex);
+  }
+
+  // Update the rest of the matrix for calculations on subsequent rows
+  for (IdComponent rowIndex = topCornerIndex + 1; rowIndex < Size; rowIndex++)
+  {
+    for (IdComponent colIndex = topCornerIndex + 1; colIndex < Size; colIndex++)
+    {
+      m(rowIndex, colIndex) -= m(rowIndex, topCornerIndex) * m(topCornerIndex, colIndex);
+    }
+  }
+}
+
+/// Performs an LUP-factorization on the given matrix using Crout's method. The
+/// LU-factorization takes a matrix A and decomposes it into a lower triangular
+/// matrix L and upper triangular matrix U such that A = LU. The
+/// LUP-factorization also allows permutation of A, which makes the
+/// decomposition always possible so long as A is not singular. In addition to
+/// matrices L and U, LUP also finds permutation matrix P containing all zeros
+/// except one 1 per row and column such that PA = LU.
+///
+/// The result is done in place such that the lower triangular matrix, L, is
+/// stored in the lower-left triangle of A including the diagonal. The upper
+/// triangular matrix, U, is stored in the upper-right triangle of L not
+/// including the diagonal. The diagonal of U in Crout's method is all 1's (and
+/// therefore not explicitly stored).
+///
+/// The permutation matrix P is represented by the permutation vector. If
+/// permutation[i] = j then row j in the original matrix A has been moved to
+/// row i in the resulting matrices. The permutation matrix P can be
+/// represented by a matrix with p_i,j = 1 if permutation[i] = j and 0
+/// otherwise. If using LUP-factorization to compute a determinant, you also
+/// need to know the parity (whether there is an odd or even amount) of
+/// inversions. An inversion is an instance of a smaller number appearing after
+/// a larger number in permutation. Although you can compute the inversion
+/// parity after the fact, this function keeps track of it with much less
+/// compute resources. The parameter inversionParity is set to 1.0 for even
+/// parity and -1.0 for odd parity.
+///
+/// Not all matrices (specifically singular matrices) have an
+/// LUP-factorization. If the LUP-factorization succeeds, valid is set to true.
+/// Otherwise, valid is set to false and the result is indeterminant.
+///
+template <typename T, IdComponent Size>
+LCL_EXEC
+lcl::ErrorCode matrixLUPFactor(Matrix<T, Size, Size>& m,
+                                Vector<IdComponent, Size>& permutation,
+                                T& inversionParity) noexcept
+{
+  // Initialize permutation.
+  for (IdComponent index = 0; index < Size; index++)
+  {
+    permutation[index] = index;
+  }
+  inversionParity = T(1);
+
+  for (IdComponent rowIndex = 0; rowIndex < Size; rowIndex++)
+  {
+    LCL_RETURN_ON_ERROR(matrixLUPFactorFindPivot(m, permutation, rowIndex, inversionParity))
+    matrixLUPFactorFindUpperTriangleElements(m, rowIndex);
+  }
+
+  return lcl::ErrorCode::SUCCESS;
+}
+
+/// Use a previous factorization done with MatrixLUPFactor to solve the
+/// system Ax = b.  Instead of A, this method takes in the LU and P
+/// matrices calculated by MatrixLUPFactor from A. The x matrix is returned.
+///
+template <typename T, IdComponent Size>
+LCL_EXEC
+Vector<T, Size> matrixLUPSolve(
+  const Matrix<T, Size, Size>& LU,
+  const Vector<IdComponent, Size>& permutation,
+  const Vector<T, Size>& b) noexcept
+{
+  // The LUP-factorization gives us PA = LU or equivalently A = inv(P)LU.
+  // Substituting into Ax = b gives us inv(P)LUx = b or LUx = Pb.
+  // Now consider the intermediate vector y = Ux.
+  // Substituting in the previous two equations yields Ly = Pb.
+  // Solving Ly = Pb is easy because L is triangular and P is just a
+  // permutation.
+  Vector<T, Size> y;
+  for (IdComponent rowIndex = 0; rowIndex < Size; rowIndex++)
+  {
+    y[rowIndex] = b[permutation[rowIndex]];
+    // Recall that L is stored in the lower triangle of LU including diagonal.
+    for (IdComponent colIndex = 0; colIndex < rowIndex; colIndex++)
+    {
+      y[rowIndex] -= LU(rowIndex, colIndex) * y[colIndex];
+    }
+    y[rowIndex] /= LU(rowIndex, rowIndex);
+  }
+
+  // Now that we have y, we can easily solve Ux = y for x.
+  Vector<T, Size> x;
+  for (IdComponent rowIndex = Size - 1; rowIndex >= 0; rowIndex--)
+  {
+    // Recall that U is stored in the upper triangle of LU with the diagonal
+    // implicitly all 1's.
+    x[rowIndex] = y[rowIndex];
+    for (IdComponent colIndex = rowIndex + 1; colIndex < Size; colIndex++)
+    {
+      x[rowIndex] -= LU(rowIndex, colIndex) * x[colIndex];
+    }
+  }
+
+  return x;
+}
+
+} // namespace detail
+
+template <typename T, IdComponent Size>
+LCL_EXEC
+lcl::ErrorCode solveLinearSystem(const Matrix<T, Size, Size>& A,
+                                  const Vector<T, Size>& b,
+                                  Vector<T, Size>& x) noexcept
+{
+  // First, we will make an LUP-factorization to help us.
+  Matrix<T, Size, Size> LU = A;
+  Vector<IdComponent, Size> permutation;
+  T inversionParity; // Unused.
+  LCL_RETURN_ON_ERROR(detail::matrixLUPFactor(LU, permutation, inversionParity))
+
+  // Next, use the decomposition to solve the system.
+  x = detail::matrixLUPSolve(LU, permutation, b);
+  return lcl::ErrorCode::SUCCESS;
+}
+
+/// Find and return the inverse of the given matrix. If the matrix is singular,
+/// the inverse will not be correct and valid will be set to false.
+///
+template <typename T, IdComponent Size>
+LCL_EXEC
+lcl::ErrorCode matrixInverse(const Matrix<T, Size, Size>& A, Matrix<T, Size, Size>& invA) noexcept
+{
+  // First, we will make an LUP-factorization to help us.
+  Matrix<T, Size, Size> LU = A;
+  Vector<IdComponent, Size> permutation;
+  T inversionParity; // Unused
+  LCL_RETURN_ON_ERROR(detail::matrixLUPFactor(LU, permutation, inversionParity))
+
+  // We will use the decomposition to solve AX = I for X where X is
+  // clearly the inverse of A.  Our solve method only works for vectors,
+  // so we solve for one column of invA at a time.
+  Vector<T, Size> ICol(T(0));
+  for (IdComponent colIndex = 0; colIndex < Size; colIndex++)
+  {
+    ICol[colIndex] = 1;
+    Vector<T, Size> invACol = detail::matrixLUPSolve(LU, permutation, ICol);
+    ICol[colIndex] = 0;
+    for (IdComponent i = 0; i < Size; ++i)
+    {
+      invA(i, colIndex) = invACol[i];
+    }
+  }
+  return lcl::ErrorCode::SUCCESS;
+}
+
+///=========================================================================
+template <typename JacobianFunctor, typename FunctionFunctor, typename T, IdComponent Size>
+LCL_EXEC
+inline lcl::ErrorCode newtonsMethod(
+  const JacobianFunctor& jacobianEvaluator,
+  const FunctionFunctor& functionEvaluator,
+  const Vector<T, Size>& rhs,
+  Vector<T, Size>& result,
+  T convergeDifference = 1e-3f,
+  int maxIterations = 10) noexcept
+{
+  auto x = result;  // pass initial guess in result
+  bool converged = false;
+
+  for (int i = 0; !converged && i < maxIterations; ++i)
+  {
+    Matrix<T, Size, Size> jacobian;
+    Vector<T, Size> fx(0);
+
+    LCL_RETURN_ON_ERROR(jacobianEvaluator(x, jacobian))
+    LCL_RETURN_ON_ERROR(functionEvaluator(x, fx))
+
+    Vector<T, Size> deltax;
+    LCL_RETURN_ON_ERROR(solveLinearSystem(jacobian, fx - rhs, deltax))
+    x -= deltax;
+
+    converged = true;
+    for (int c = 0; c < Size; ++c)
+    {
+      converged &= (LCL_MATH_CALL(abs, (deltax[c])) < convergeDifference);
+    }
+  }
+
+  result = x;
+  return converged ? ErrorCode::SUCCESS : ErrorCode::SOLUTION_DID_NOT_CONVERGE;
+}
+
+///=========================================================================
+template <typename T>
+LCL_EXEC inline T lerp(T v0, T v1, T t)
+{
+  static_assert(std::is_floating_point<T>::value,
+                "lerp requires floating point parameters");
+
+  return LCL_MATH_CALL(fma, (t), (v1), (LCL_MATH_CALL(fma, (-t), (v0), (v0))));
+}
+
+}
+} // lcl::internal
+
+#endif // lcl_internal_Math_h

vtkm/thirdparty/lcl/vtkmlcl/lcl/lcl.h

@@ -0,0 +1,251 @@
+//============================================================================
+//  Copyright (c) Kitware, Inc.
+//  All rights reserved.
+//  See LICENSE.md 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.
+//============================================================================
+#ifndef lcl_lcl_h
+#define lcl_lcl_h
+
+#include <lcl/Hexahedron.h>
+#include <lcl/Line.h>
+#include <lcl/Pixel.h>
+#include <lcl/Polygon.h>
+#include <lcl/Pyramid.h>
+#include <lcl/Quad.h>
+#include <lcl/Tetra.h>
+#include <lcl/Triangle.h>
+#include <lcl/Vertex.h>
+#include <lcl/Voxel.h>
+#include <lcl/Wedge.h>
+
+#include <utility>
+
+namespace lcl
+{
+
+/// \brief Perform basic checks to validate cell's state.
+/// \param[in]  tag  The cell tag to validate.
+/// \return          lcl::ErrorCode::SUCCESS if valid.
+///
+LCL_EXEC inline lcl::ErrorCode validate(Cell tag) noexcept
+{
+  ErrorCode status = ErrorCode::SUCCESS;
+  switch (tag.shape())
+  {
+    lclGenericCellShapeMacro(status = validate(CellTag{tag}));
+    default:
+      status = ErrorCode::INVALID_SHAPE_ID;
+      break;
+  }
+  return status;
+}
+
+/// \brief Return center of a cell in parametric coordinates.
+/// \remark Note that the parametric center is not always located at (0.5,0.5,0.5).
+/// \param[in]   tag      The cell tag.
+/// \param[out]  pcoords  The center of the cell in parametric coordinates.
+/// \return               A lcl::ErrorCode value indicating the status of the operation.
+///
+template<typename CoordType>
+LCL_EXEC inline lcl::ErrorCode parametricCenter(Cell tag, CoordType&& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  ErrorCode status = ErrorCode::SUCCESS;
+  switch (tag.shape())
+  {
+    lclGenericCellShapeMacro(status = parametricCenter(CellTag{tag}, std::forward<CoordType>(pcoords)));
+    default:
+      status = ErrorCode::INVALID_SHAPE_ID;
+      break;
+  }
+  return status;
+}
+
+/// \brief Return the parametric coordinates of a cell's point.
+/// \param[in]   tag      The cell tag.
+/// \param[in]   pointId  The point number.
+/// \param[out]  pcoords  The parametric coordinates of a cell's point.
+/// \return               A lcl::ErrorCode value indicating the status of the operation.
+///
+template<typename CoordType>
+LCL_EXEC inline lcl::ErrorCode parametricPoint(
+  Cell tag, IdComponent pointId, CoordType&& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  ErrorCode status = ErrorCode::SUCCESS;
+  switch (tag.shape())
+  {
+    lclGenericCellShapeMacro(status = parametricPoint(CellTag{tag}, pointId, std::forward<CoordType>(pcoords)));
+    default:
+      status = ErrorCode::INVALID_SHAPE_ID;
+      break;
+  }
+  return status;
+}
+
+/// \brief Return the parametric distance of a parametric coordinate to a cell.
+/// \param[in]  tag      The cell tag.
+/// \param[in]  pcoords  The parametric coordinates of the point.
+/// \return              The parametric distance of the point to the cell.
+///                      If point is inside the cell, 0 is returned.
+/// \pre tag should be a valid cell, otherwise the result is undefined.
+///
+template<typename CoordType>
+LCL_EXEC inline ComponentType<CoordType> parametricDistance(Cell tag, const CoordType& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  ComponentType<CoordType> dist{0};
+  switch (tag.shape())
+  {
+    lclGenericCellShapeMacro(dist = parametricDistance(CellTag{tag}, pcoords));
+    default:
+      break;
+  }
+  return dist;
+}
+
+/// \brief Check if the given parametric point lies inside a cell.
+/// \param[in]  tag      The cell tag.
+/// \param[in]  pcoords  The parametric coordinates of the point.
+/// \return              true if inside, false otherwise.
+/// \pre tag should be a valid cell, otherwise the result is undefined.
+///
+template<typename CoordType>
+LCL_EXEC inline bool cellInside(Cell tag, const CoordType& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  bool inside = false;
+  switch (tag.shape())
+  {
+    lclGenericCellShapeMacro(inside = cellInside(CellTag{tag}, pcoords));
+    default:
+      break;
+  }
+  return inside;
+}
+
+/// \brief Interpolate \c values at the paramteric coordinates \c pcoords
+/// \param[in]   tag      The cell tag.
+/// \param[in]   values   A \c FieldAccessor for values to interpolate
+/// \param[in]   pcoords  The parametric coordinates.
+/// \param[out]  result   The interpolation result
+/// \return               A lcl::ErrorCode value indicating the status of the operation.
+///
+template <typename Values, typename CoordType, typename Result>
+LCL_EXEC inline lcl::ErrorCode interpolate(
+  Cell tag,
+  const Values& values,
+  const CoordType& pcoords,
+  Result&& result) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  ErrorCode status = ErrorCode::SUCCESS;
+  switch (tag.shape())
+  {
+    lclGenericCellShapeMacro(status = interpolate(CellTag{tag}, values, pcoords, std::forward<Result>(result)));
+    default:
+      status = ErrorCode::INVALID_SHAPE_ID;
+  }
+  return status;
+}
+
+/// \brief Compute derivative of \c values at the paramteric coordinates \c pcoords
+/// \param[in]   tag      The cell tag.
+/// \param[in]   points   A \c FieldAccessor for points of the cell
+/// \param[in]   values   A \c FieldAccessor for the values to compute derivative of.
+/// \param[in]   pcoords  The parametric coordinates.
+/// \param[out]  dx       The derivative along X
+/// \param[out]  dy       The derivative along Y
+/// \param[out]  dz       The derivative along Z
+/// \return               A lcl::ErrorCode value indicating the status of the operation.
+///
+template <typename Points, typename Values, typename CoordType, typename Result>
+LCL_EXEC inline lcl::ErrorCode derivative(
+  Cell tag,
+  const Points& points,
+  const Values& values,
+  const CoordType& pcoords,
+  Result&& dx,
+  Result&& dy,
+  Result&& dz) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  ErrorCode status = ErrorCode::SUCCESS;
+  switch (tag.shape())
+  {
+    lclGenericCellShapeMacro(status = derivative(CellTag{tag},
+                                                points,
+                                                values,
+                                                pcoords,
+                                                std::forward<Result>(dx),
+                                                std::forward<Result>(dy),
+                                                std::forward<Result>(dz)));
+    default:
+      status = ErrorCode::INVALID_SHAPE_ID;
+  }
+  return status;
+}
+
+/// \brief Compute world coordinates from parametric coordinates
+/// \param[in]   tag      The cell tag.
+/// \param[in]   points   A \c FieldAccessor for points of the cell
+/// \param[in]   pcoords  The parametric coordinates.
+/// \param[out]  wcoords  The world coordinates.
+///
+template <typename Points, typename PCoordType, typename WCoordType>
+LCL_EXEC inline lcl::ErrorCode parametricToWorld(
+  Cell tag,
+  const Points& points,
+  const PCoordType& pcoords,
+  WCoordType&& wcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(PCoordType);
+
+  ErrorCode status = ErrorCode::SUCCESS;
+  switch (tag.shape())
+  {
+    lclGenericCellShapeMacro(status = parametricToWorld(CellTag{tag}, points, pcoords, std::forward<WCoordType>(wcoords)));
+    default:
+      status = ErrorCode::INVALID_SHAPE_ID;
+  }
+  return status;
+}
+
+/// \brief Compute parametric coordinates from world coordinates
+/// \param[in]   tag      The cell tag.
+/// \param[in]   points   A \c FieldAccessor for points of the cell
+/// \param[in]   wcoords  The world coordinates.
+/// \param[out]  pcoords  The parametric coordinates.
+///
+template <typename Points, typename WCoordType, typename PCoordType>
+LCL_EXEC inline lcl::ErrorCode worldToParametric(
+  Cell tag,
+  const Points& points,
+  const WCoordType& wcoords,
+  PCoordType&& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(PCoordType);
+
+  ErrorCode status = ErrorCode::SUCCESS;
+  switch (tag.shape())
+  {
+    lclGenericCellShapeMacro(status = worldToParametric(CellTag{tag}, points, wcoords, std::forward<PCoordType>(pcoords)));
+    default:
+      status = ErrorCode::INVALID_SHAPE_ID;
+  }
+  return status;
+}
+
+} //namespace lcl
+
+#endif //lcl_lcl_h