Merge branch 'jit-lagrange' into 'master'

WIP: Jit lagrange

See merge request vtk/vtk-m!2287

vtkm/CellShape.h

@@ -13,6 +13,7 @@
 #include <vtkm/StaticAssert.h>
 #include <vtkm/Types.h>
 
+#include <lcl/Lagrange_Hexahedron.h>
 #include <lcl/Polygon.h>
 #include <lcl/Shapes.h>
 
@@ -66,6 +67,7 @@ enum CellShapeIdEnum
   CELL_SHAPE_WEDGE = lcl::ShapeId::WEDGE,
   /// A pyramid with a quadrilateral base and four triangular faces.0
   CELL_SHAPE_PYRAMID = lcl::ShapeId::PYRAMID,
+  CELL_SHAPE_LAGRANGE_HEXAHEDRON = lcl::ShapeId::LAGRANGE_HEXAHEDRON,
 
   NUMBER_OF_CELL_SHAPES
 };
@@ -159,6 +161,7 @@ VTKM_DEFINE_CELL_TAG(Tetra, CELL_SHAPE_TETRA);
 VTKM_DEFINE_CELL_TAG(Hexahedron, CELL_SHAPE_HEXAHEDRON);
 VTKM_DEFINE_CELL_TAG(Wedge, CELL_SHAPE_WEDGE);
 VTKM_DEFINE_CELL_TAG(Pyramid, CELL_SHAPE_PYRAMID);
+VTKM_DEFINE_CELL_TAG(Lagrange_Hexahedron, CELL_SHAPE_LAGRANGE_HEXAHEDRON);
 
 #undef VTKM_DEFINE_CELL_TAG
 
@@ -200,6 +203,13 @@ inline lcl::Polygon make_LclCellShapeTag(const vtkm::CellShapeTagPolygon&,
   return lcl::Polygon(numPoints);
 }
 
+VTKM_EXEC_CONT
+inline lcl::Lagrange_Hexahedron make_LclCellShapeTag(const vtkm::CellShapeTagLagrange_Hexahedron&,
+                                                     vtkm::IdComponent numPoints = 0)
+{
+  return lcl::Lagrange_Hexahedron(numPoints);
+}
+
 VTKM_EXEC_CONT
 inline lcl::Cell make_LclCellShapeTag(const vtkm::CellShapeTagGeneric& tag,
                                       vtkm::IdComponent numPoints = 0)
@@ -258,7 +268,8 @@ inline lcl::Cell make_LclCellShapeTag(const vtkm::CellShapeTagGeneric& tag,
   vtkmGenericCellShapeMacroCase(CELL_SHAPE_TETRA, call);      \
   vtkmGenericCellShapeMacroCase(CELL_SHAPE_HEXAHEDRON, call); \
   vtkmGenericCellShapeMacroCase(CELL_SHAPE_WEDGE, call);      \
-  vtkmGenericCellShapeMacroCase(CELL_SHAPE_PYRAMID, call)
+  vtkmGenericCellShapeMacroCase(CELL_SHAPE_PYRAMID, call);    \
+  vtkmGenericCellShapeMacroCase(CELL_SHAPE_LAGRANGE_HEXAHEDRON, call);
 
 } // namespace vtkm
 

vtkm/CellTraits.h

@@ -115,6 +115,7 @@ VTKM_DEFINE_CELL_TRAITS(Tetra, 3, 4);
 VTKM_DEFINE_CELL_TRAITS(Hexahedron, 3, 8);
 VTKM_DEFINE_CELL_TRAITS(Wedge, 3, 6);
 VTKM_DEFINE_CELL_TRAITS(Pyramid, 3, 5);
+VTKM_DEFINE_CELL_TRAITS_VARIABLE(Lagrange_Hexahedron, 3);
 
 #undef VTKM_DEFINE_CELL_TRAITS
 

vtkm/VecFromPortalPermute.h

@@ -40,11 +40,20 @@ public:
     : Indices(indices)
     , Portal(portal)
   {
+    Limit = 0;
   }
 
+  void LimitNumberOfComponents(vtkm::IdComponent l) { Limit = l; }
+
   VTKM_SUPPRESS_EXEC_WARNINGS
   VTKM_EXEC_CONT
-  vtkm::IdComponent GetNumberOfComponents() const { return this->Indices->GetNumberOfComponents(); }
+  vtkm::IdComponent GetNumberOfComponents() const
+  {
+    vtkm::IdComponent numComps = this->Indices->GetNumberOfComponents();
+    if (Limit > 0)
+      return (Limit < numComps ? Limit : numComps);
+    return numComps;
+  }
 
   VTKM_SUPPRESS_EXEC_WARNINGS
   template <vtkm::IdComponent DestSize>
@@ -67,6 +76,7 @@ public:
 private:
   const IndexVecType* const Indices;
   PortalType Portal;
+  vtkm::IdComponent Limit;
 };
 
 template <typename IndexVecType, typename PortalType>

vtkm/cont/testlib/MakeTestDataSet.cxx

@@ -233,33 +233,56 @@ vtkm::cont::DataSet MakeTestDataSet::Make3DUniformDataSet1()
   constexpr vtkm::Id nVerts = 125;
   constexpr vtkm::Id nCells = 64;
   constexpr vtkm::Float32 pointvar[nVerts] = {
-    0.0f,  0.0f, 0.0f, 0.0f,  0.0f,  0.0f,  0.0f,  0.0f,  0.0f,  0.0f, 0.0f, 0.0f,  0.0f,
-    0.0f,  0.0f, 0.0f, 0.0f,  0.0f,  0.0f,  0.0f,  0.0f,  0.0f,  0.0f, 0.0f, 0.0f,
+    0.0f, 0.0f,  0.0f,  0.0f,  0.0f, // 0-4, 0, 0
+    0.0f, 0.0f,  0.0f,  0.0f,  0.0f, // 0-4, 1, 0
+    0.0f, 0.0f,  0.0f,  0.0f,  0.0f, // 0-4, 2, 0
+    0.0f, 0.0f,  0.0f,  0.0f,  0.0f, // 0-4, 3, 0
+    0.0f, 0.0f,  0.0f,  0.0f,  0.0f, // 0-4, 4, 0
 
-    0.0f,  0.0f, 0.0f, 0.0f,  0.0f,  0.0f,  99.0f, 90.0f, 85.0f, 0.0f, 0.0f, 95.0f, 80.0f,
-    95.0f, 0.0f, 0.0f, 85.0f, 90.0f, 99.0f, 0.0f,  0.0f,  0.0f,  0.0f, 0.0f, 0.0f,
+    0.0f, 0.0f,  0.0f,  0.0f,  0.0f, // 0-4, 0, 1
+    0.0f, 99.0f, 90.0f, 85.0f, 0.0f, // 0-4, 1, 1
+    0.0f, 95.0f, 80.0f, 95.0f, 0.0f, // 0-4, 2, 1
+    0.0f, 85.0f, 90.0f, 99.0f, 0.0f, // 0-4, 3, 1
+    0.0f, 0.0f,  0.0f,  0.0f,  0.0f, // 0-4, 4, 1
 
-    0.0f,  0.0f, 0.0f, 0.0f,  0.0f,  0.0f,  75.0f, 50.0f, 65.0f, 0.0f, 0.0f, 55.0f, 15.0f,
-    45.0f, 0.0f, 0.0f, 60.0f, 40.0f, 70.0f, 0.0f,  0.0f,  0.0f,  0.0f, 0.0f, 0.0f,
+    0.0f, 0.0f,  0.0f,  0.0f,  0.0f, // 0-4, 0, 2
+    0.0f, 75.0f, 50.0f, 65.0f, 0.0f, // 0-4, 1, 2
+    0.0f, 55.0f, 15.0f, 45.0f, 0.0f, // 0-4, 2, 2
+    0.0f, 60.0f, 40.0f, 70.0f, 0.0f, // 0-4, 3, 2
+    0.0f, 0.0f,  0.0f,  0.0f,  0.0f, // 0-4, 4, 2
 
-    0.0f,  0.0f, 0.0f, 0.0f,  0.0f,  0.0f,  97.0f, 87.0f, 82.0f, 0.0f, 0.0f, 92.0f, 77.0f,
-    92.0f, 0.0f, 0.0f, 82.0f, 87.0f, 97.0f, 0.0f,  0.0f,  0.0f,  0.0f, 0.0f, 0.0f,
+    0.0f, 0.0f,  0.0f,  0.0f,  0.0f, // 0-4, 0, 3
+    0.0f, 97.0f, 87.0f, 82.0f, 0.0f, // 0-4, 1, 3
+    0.0f, 92.0f, 77.0f, 92.0f, 0.0f, // 0-4, 2, 3
+    0.0f, 82.0f, 87.0f, 97.0f, 0.0f, // 0-4, 3, 3
+    0.0f, 0.0f,  0.0f,  0.0f,  0.0f, // 0-4, 4, 3
 
-    0.0f,  0.0f, 0.0f, 0.0f,  0.0f,  0.0f,  0.0f,  0.0f,  0.0f,  0.0f, 0.0f, 0.0f,  0.0f,
-    0.0f,  0.0f, 0.0f, 0.0f,  0.0f,  0.0f,  0.0f,  0.0f,  0.0f,  0.0f, 0.0f, 0.0f
+    0.0f, 0.0f,  0.0f,  0.0f,  0.0f, // 0-4, 0, 4
+    0.0f, 0.0f,  0.0f,  0.0f,  0.0f, // 0-4, 1, 4
+    0.0f, 0.0f,  0.0f,  0.0f,  0.0f, // 0-4, 2, 4
+    0.0f, 0.0f,  0.0f,  0.0f,  0.0f, // 0-4, 3, 4
+    0.0f, 0.0f,  0.0f,  0.0f,  0.0f  // 0-4, 4, 4
   };
   constexpr vtkm::Float32 cellvar[nCells] = {
-    0.0f,  1.0f,  2.0f,  3.0f,  4.0f,  5.0f,  6.0f,  7.0f,
-    8.0f,  9.0f,  10.0f, 11.0f, 12.0f, 13.0f, 14.0f, 15.0f,
+    0.0f,  1.0f,  2.0f,  3.0f,  // 0-3, 0, 0
+    4.0f,  5.0f,  6.0f,  7.0f,  // 0-3, 1, 0
+    8.0f,  9.0f,  10.0f, 11.0f, // 0-3, 2, 0
+    12.0f, 13.0f, 14.0f, 15.0f, // 0-3, 3, 0
 
-    16.0f, 17.0f, 18.0f, 19.0f, 20.0f, 21.0f, 22.0f, 23.0f,
-    24.0f, 25.0f, 26.0f, 27.0f, 28.0f, 29.0f, 30.0f, 31.0f,
+    16.0f, 17.0f, 18.0f, 19.0f, // 0-3, 0, 1
+    20.0f, 21.0f, 22.0f, 23.0f, // 0-3, 1, 1
+    24.0f, 25.0f, 26.0f, 27.0f, // 0-3, 2, 1
+    28.0f, 29.0f, 30.0f, 31.0f, // 0-3, 3, 1
 
-    32.0f, 33.0f, 34.0f, 35.0f, 36.0f, 37.0f, 38.0f, 39.0f,
-    40.0f, 41.0f, 42.0f, 43.0f, 44.0f, 45.0f, 46.0f, 47.0f,
+    32.0f, 33.0f, 34.0f, 35.0f, // 0-3, 0, 2
+    36.0f, 37.0f, 38.0f, 39.0f, // 0-3, 1, 2
+    40.0f, 41.0f, 42.0f, 43.0f, // 0-3, 2, 2
+    44.0f, 45.0f, 46.0f, 47.0f, // 0-3, 3, 2
 
-    48.0f, 49.0f, 50.0f, 51.0f, 52.0f, 53.0f, 54.0f, 55.0f,
-    56.0f, 57.0f, 58.0f, 59.0f, 60.0f, 61.0f, 62.0f, 63.0f
+    48.0f, 49.0f, 50.0f, 51.0f, // 0-3, 0, 3
+    52.0f, 53.0f, 54.0f, 55.0f, // 0-3, 1, 3
+    56.0f, 57.0f, 58.0f, 59.0f, // 0-3, 2, 3
+    60.0f, 61.0f, 62.0f, 63.0f  // 0-3, 3, 3
   };
 
   dataSet.AddPointField("pointvar", pointvar, nVerts);

vtkm/exec/CellEdge.h

@@ -160,6 +160,19 @@ static inline VTKM_EXEC vtkm::ErrorCode CellEdgeNumberOfEdges(vtkm::IdComponent
 /// @param[in]  shape A tag of type `CellShapeTag*` to identify the shape of the cell.
 ///     This method is overloaded for different shape types.
 /// @param[out] numEdges A reference to return the number of edges.
+static inline VTKM_EXEC vtkm::ErrorCode CellEdgeNumberOfEdges(vtkm::IdComponent numPoints,
+                                                              vtkm::CellShapeTagLagrange_Hexahedron,
+                                                              vtkm::IdComponent& numEdges)
+{
+  if (numPoints <= 0)
+  {
+    numEdges = -1;
+    return vtkm::ErrorCode::InvalidNumberOfPoints;
+  }
+  numEdges = numPoints; // HC: this is what SW did, but possibly it was a shortcut
+  return vtkm::ErrorCode::Success;
+}
+
 static inline VTKM_EXEC vtkm::ErrorCode CellEdgeNumberOfEdges(vtkm::IdComponent numPoints,
                                                               vtkm::CellShapeTagGeneric shape,
                                                               vtkm::IdComponent& numEdges)

vtkm/exec/arg/FetchTagArrayTopologyMapIn.h

@@ -19,6 +19,7 @@
 #include <vtkm/internal/ArrayPortalUniformPointCoordinates.h>
 
 #include <vtkm/VecAxisAlignedPointCoordinates.h>
+#include <vtkm/exec/ConnectivityExplicit.h>
 #include <vtkm/exec/ConnectivityExtrude.h>
 #include <vtkm/exec/ConnectivityStructured.h>
 
@@ -154,6 +155,67 @@ struct FetchArrayTopologyMapInImplementation<
   }
 };
 
+template <typename PortalType>
+VTKM_EXEC inline void PerformHigherOrderToLowerOrderCellSubstitution(
+  PortalType pt,
+  vtkm::CellShapeTagLagrange_Hexahedron)
+{
+  pt.LimitNumberOfComponents(8);
+}
+
+template <typename PortalType>
+VTKM_EXEC inline void PerformHigherOrderToLowerOrderCellSubstitution(
+  PortalType pt,
+  vtkm::CellShapeTagGeneric shape)
+{
+  if (shape.Id == vtkm::CELL_SHAPE_LAGRANGE_HEXAHEDRON)
+    pt.LimitNumberOfComponents(8);
+}
+
+
+template <typename PortalType, typename CellShapeTag>
+VTKM_EXEC inline void PerformHigherOrderToLowerOrderCellSubstitution(PortalType, CellShapeTag)
+{
+  return;
+}
+
+template <typename PortalType, typename CellShapeTag>
+VTKM_EXEC inline void PerformInlineShapeSubstitution(PortalType pt, CellShapeTag shape)
+{
+  PerformHigherOrderToLowerOrderCellSubstitution(pt, shape);
+}
+
+template <typename FieldExecObjectType, typename X, typename Y, typename Z, typename W>
+struct FetchArrayTopologyMapInImplementation<
+  vtkm::exec::ConnectivityExplicit<X, Y, Z>,
+  FieldExecObjectType,
+  ThreadIndicesTopologyMap<vtkm::exec::ConnectivityExplicit<X, Y, Z>, W>>
+{
+  using ThreadIndicesTopologyMapType =
+    const ThreadIndicesTopologyMap<vtkm::exec::ConnectivityExplicit<X, Y, Z>, W>;
+
+  // stored in a Vec-like object.
+  using IndexVecType = typename ThreadIndicesTopologyMapType::IndicesIncidentType;
+
+  // The FieldExecObjectType is expected to behave like an ArrayPortal.
+  using PortalType = FieldExecObjectType;
+
+  using ValueType = vtkm::VecFromPortalPermute<IndexVecType, PortalType>;
+
+  VTKM_SUPPRESS_EXEC_WARNINGS
+  VTKM_EXEC
+  static ValueType Load(ThreadIndicesTopologyMapType& indices, const FieldExecObjectType& field)
+  {
+    std::cerr << "LOAD5" << std::endl;
+    ValueType rv(indices.GetIndicesIncidentPointer(), field);
+    PerformInlineShapeSubstitution(rv, indices.GetCellShape());
+    //if (indices.GetCellShape() == vtkm::CellShapeTagLagrange_Hexahedron())
+    //rv.LimitNumberOfComponents(8);
+    return rv;
+  }
+};
+
+
 template <typename PermutationPortal, vtkm::IdComponent NumDimensions, typename ThreadIndicesType>
 struct FetchArrayTopologyMapInImplementation<
   vtkm::exec::ConnectivityPermutedVisitCellsWithPoints<

vtkm/exec/testing/UnitTestCellEdgeFace.cxx

@@ -221,6 +221,16 @@ struct TestCellFacesFunctor
       this->TryShapeWithNumPoints(numPoints, vtkm::CellShapeTagPolygon());
     }
   }
+
+  void operator()(vtkm::CellShapeTagLagrange_Hexahedron) const
+  {
+    // n^3
+    for (vtkm::IdComponent numPoints = 8; numPoints < 28;
+         numPoints = static_cast<vtkm::IdComponent>(vtkm::Pow(numPoints, 3.)))
+    {
+      this->TryShapeWithNumPoints(numPoints, vtkm::CellShapeTagLagrange_Hexahedron());
+    }
+  }
 };
 
 void TestAllShapes()

vtkm/filter/entity_extraction/ExtractStructured.h

@@ -85,17 +85,23 @@ public:
   VTKM_CONT
   void SetSampleRate(vtkm::Id3 sampleRate) { this->SampleRate = sampleRate; }
 
-  // Get if we should include the outer boundary on a subsample
-  // (NOTE: This method is deprecated because it is not clear from either
-  // the documentation or the source code what the intention of this feature
-  // is. If your are using this method somewhere and think it should remain,
-  // please open a merge request to "de-deprecate" it and add a test and
-  // documentation of the expected behavior.)
-  VTKM_DEPRECATED(2.2)
+  /// @brief Specifies if the outer boundary should always be included.
+  ///
+  /// When a subsample rate is specified, it is possible that some of the
+  /// boundary samples will not be included in the sampling. If this is the
+  /// case and `IncludeBoundary` is set to true, then an extra sample is
+  /// set in the output and the values on the boundary are included. For example,
+  /// say the input has resolution (5, 5, 1) (and the VOI matches), and the sample
+  /// rate is set to (3, 3, 1). If `IncludeBoundary` is false, then the output will
+  /// have the 4 points that correspond to the 3D indices (0, 0, 0), (3, 0, 0),
+  /// (0, 3, 0), and (3, 3, 0) of the input. This misses the outer boundary at
+  /// index 4 in the x and y directions. If `IncludeBoundary is set to true, then
+  /// the output will have the 9 points that correspond to the 3D indices (0, 0, 0),
+  /// (3, 0, 0), (4, 0, 0), (0, 3, 0), (3, 3, 0), (4, 3, 0), (0, 4, 0), (3, 4, 0),
+  /// and (4, 4, 0) to capture this outer boundary.
   VTKM_CONT
   bool GetIncludeBoundary() const { return this->IncludeBoundary; }
-  // Set if we should include the outer boundary on a subsample
-  VTKM_DEPRECATED(2.2)
+  /// @copydoc GetIncludeBoundary
   VTKM_CONT
   void SetIncludeBoundary(bool value) { this->IncludeBoundary = value; }
 

vtkm/filter/entity_extraction/testing/UnitTestExtractStructuredFilter.cxx

@@ -341,7 +341,7 @@ public:
 
   static void TestUniform3D7()
   {
-    std::cout << "Testing extract structured uniform" << std::endl;
+    std::cout << "Testing extract structured uniform, exclude boundary" << std::endl;
     vtkm::cont::DataSet dataset = MakeTestDataSet().Make3DUniformDataSet1();
 
     vtkm::filter::entity_extraction::ExtractStructured extract;
@@ -351,6 +351,7 @@ public:
     vtkm::Id3 sample(3, 3, 2);
     extract.SetVOI(range);
     extract.SetSampleRate(sample);
+    extract.SetIncludeBoundary(false); // default
 
     extract.SetFieldsToPass({ "pointvar", "cellvar" });
     vtkm::cont::DataSet output = extract.Execute(dataset);
@@ -371,16 +372,16 @@ public:
                      "Data/Geometry mismatch for ExtractStructured filter");
 
     VTKM_TEST_ASSERT(outPointData.ReadPortal().Get(0) == 0.0f, "Wrong point field data");
+    VTKM_TEST_ASSERT(outPointData.ReadPortal().Get(3) == 99.0f, "Wrong point field data");
+    VTKM_TEST_ASSERT(outPointData.ReadPortal().Get(4) == 0.0f, "Wrong point field data");
     VTKM_TEST_ASSERT(outPointData.ReadPortal().Get(7) == 97.0f, "Wrong point field data");
 
     VTKM_TEST_ASSERT(outCellData.ReadPortal().Get(0) == 16.0f, "Wrong cell field data");
   }
 
-  VTKM_DEPRECATED_SUPPRESS_BEGIN
-  // This test repeates TestUniform3D7 but uses deprecated behavior.
   static void TestUniform3D8()
   {
-    std::cout << "Testing extract structured uniform" << std::endl;
+    std::cout << "Testing extract structured uniform, include boundary" << std::endl;
     vtkm::cont::DataSet dataset = MakeTestDataSet().Make3DUniformDataSet1();
     vtkm::filter::entity_extraction::ExtractStructured extract;
 
@@ -411,12 +412,15 @@ public:
 
     VTKM_TEST_ASSERT(outPointData.ReadPortal().Get(0) == 0.0f, "Wrong point field data");
     VTKM_TEST_ASSERT(outPointData.ReadPortal().Get(4) == 99.0f, "Wrong point field data");
+    VTKM_TEST_ASSERT(outPointData.ReadPortal().Get(5) == 0.0f, "Wrong point field data");
+    VTKM_TEST_ASSERT(outPointData.ReadPortal().Get(7) == 0.0f, "Wrong point field data");
     VTKM_TEST_ASSERT(outPointData.ReadPortal().Get(13) == 97.0f, "Wrong point field data");
 
     VTKM_TEST_ASSERT(outCellData.ReadPortal().Get(0) == 16.0f, "Wrong cell field data");
+    VTKM_TEST_ASSERT(outCellData.ReadPortal().Get(1) == 19.0f, "Wrong cell field data");
+    VTKM_TEST_ASSERT(outCellData.ReadPortal().Get(2) == 28.0f, "Wrong cell field data");
     VTKM_TEST_ASSERT(outCellData.ReadPortal().Get(3) == 31.0f, "Wrong cell field data");
   }
-  VTKM_DEPRECATED_SUPPRESS_END
 
   static void TestRectilinear2D()
   {

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

@@ -0,0 +1,274 @@
+//============================================================================
+//  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_Lagrange_Hexahedron_h
+#define lcl_Lagrange_Hexahedron_h
+
+#include <vtkm/Math.h>
+#include <lcl/ErrorCode.h>
+#include <lcl/Shapes.h>
+
+#include <lcl/internal/Common.h>
+
+namespace lcl
+{
+
+class Lagrange_Hexahedron : public Cell
+{
+public:
+  constexpr LCL_EXEC Lagrange_Hexahedron() : Cell(ShapeId::LAGRANGE_HEXAHEDRON, 8) {}
+  constexpr LCL_EXEC explicit Lagrange_Hexahedron(lcl::IdComponent numPoints)
+      : Cell(ShapeId::LAGRANGE_HEXAHEDRON, numPoints)
+  {
+  }
+  constexpr LCL_EXEC explicit Lagrange_Hexahedron(const Cell& cell) : Cell(cell) {}
+};
+
+LCL_EXEC inline lcl::ErrorCode validate(Lagrange_Hexahedron tag) noexcept
+{
+  if (tag.shape() != ShapeId::LAGRANGE_HEXAHEDRON && tag.shape() != ShapeId::VOXEL)
+  {
+    return ErrorCode::WRONG_SHAPE_ID_FOR_TAG_TYPE;
+  }
+  if (static_cast<lcl::IdComponent>(vtkm::Cbrt(tag.numberOfPoints()))%3 != 0) // 0th order Lagrange is a Hex
+  {
+    return ErrorCode::INVALID_NUMBER_OF_POINTS;
+  }
+
+  return ErrorCode::SUCCESS;
+}
+
+template<typename CoordType>
+LCL_EXEC inline lcl::ErrorCode parametricCenter(Lagrange_Hexahedron, CoordType&& pcoords) noexcept
+{
+  LCL_STATIC_ASSERT_PCOORDS_IS_FLOAT_TYPE(CoordType);
+
+  // All Lagrange Unit Hexes have the same parametric center
+  component(pcoords, 0) = 0.5f;
+  component(pcoords, 1) = 0.5f;
+  component(pcoords, 2) = 0.5f;
+  return ErrorCode::SUCCESS;
+}
+
+/*
+ * Ordering inspired by vtk quadratic hex
+ * First 8 points are the normal unit hex
+ * Next 8 points are the midpoints for the unit hex (starts <0.5,0,0>)
+ * Next 4 points is unit square at z=0.5 (starts <0,0,0.5>)
+ * Next 4 points are the midpoints of unit square at z=0.5 (starts <0.5,0,0.5>)
+ * Final 3 points are the vertical central line (<0.5,0.5,0> -> <0.5,0.5,1>)
+ */
+template<typename CoordType>
+LCL_EXEC inline lcl::ErrorCode parametricPoint(
+  Lagrange_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(Lagrange_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(Lagrange_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);
+}
+
+// TODO: How does this work with Lagrange?
+// Implementation should ignore midpoints and approximate as linear cell
+template <typename Values, typename CoordType, typename Result>
+LCL_EXEC inline lcl::ErrorCode interpolate(
+  Lagrange_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
+{
+
+// TODO: Fix for 27
+// Current is Linear?
+template <typename Values, typename CoordType, typename Result>
+LCL_EXEC inline void parametricDerivative(Lagrange_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(
+  Lagrange_Hexahedron,
+  const Points& points,
+  const Values& values,
+  const CoordType& pcoords,
+  Result&& dx,
+  Result&& dy,
+  Result&& dz) noexcept
+{
+  return internal::derivative3D(Lagrange_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(
+  Lagrange_Hexahedron,
+  const Points& points,
+  const PCoordType& pcoords,
+  WCoordType&& wcoords) noexcept
+{
+  return interpolate(Lagrange_Hexahedron{}, points, pcoords, std::forward<WCoordType>(wcoords));
+}
+
+template <typename Points, typename WCoordType, typename PCoordType>
+LCL_EXEC inline lcl::ErrorCode worldToParametric(
+  Lagrange_Hexahedron,
+  const Points& points,
+  const WCoordType& wcoords,
+  PCoordType&& pcoords) noexcept
+{
+  return internal::worldToParametric3D(
+    Lagrange_Hexahedron{}, points, wcoords, std::forward<PCoordType>(pcoords));
+}
+
+} // lcl
+
+#endif // lcl_Lagrange_Hexahedron_h

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

@@ -35,7 +35,7 @@ enum ShapeId : IdShape
   HEXAHEDRON       = 12,
   WEDGE            = 13,
   PYRAMID          = 14,
-
+  LAGRANGE_HEXAHEDRON = 72, 
   NUMBER_OF_CELL_SHAPES
 };
 
@@ -87,6 +87,7 @@ inline LCL_EXEC int dimension(IdShape shapeId)
     case HEXAHEDRON:
     case WEDGE:
     case PYRAMID:
+    case LAGRANGE_HEXAHEDRON:
       return 3;
     case EMPTY:
     default:
@@ -115,6 +116,7 @@ class Voxel;
 class Hexahedron;
 class Wedge;
 class Pyramid;
+class Lagrange_Hexahedron;
 
 } //namespace lcl
 
@@ -137,6 +139,7 @@ class Pyramid;
   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)
+  lclGenericCellShapeMacroCase(lcl::ShapeId::PYRAMID,    lcl::Pyramid,    call);                \
+  lclGenericCellShapeMacroCase(lcl::ShapeId::LAGRANGE_HEXAHEDRON,    lcl::Lagrange_Hexahedron,    call)
 
 #endif //lcl_Shapes_h

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

@@ -66,6 +66,7 @@ FORWARD_DECLAR_PARAMETRIC_DERIVATIVE(Tetra);
 FORWARD_DECLAR_PARAMETRIC_DERIVATIVE(Hexahedron);
 FORWARD_DECLAR_PARAMETRIC_DERIVATIVE(Wedge);
 FORWARD_DECLAR_PARAMETRIC_DERIVATIVE(Pyramid);
+FORWARD_DECLAR_PARAMETRIC_DERIVATIVE(Lagrange_Hexahedron);
 
 #undef FORWARD_DECLAR_PARAMETRIC_DERIVATIVE
 

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

@@ -21,6 +21,7 @@
 #include <lcl/Vertex.h>
 #include <lcl/Voxel.h>
 #include <lcl/Wedge.h>
+#include <lcl/Lagrange_Hexahedron.h>
 
 #include <utility>