diff --git a/src/MeshField.hpp b/src/MeshField.hpp
index 08d5942..2424710 100644
--- a/src/MeshField.hpp
+++ b/src/MeshField.hpp
@@ -299,7 +299,7 @@ class OmegahMeshField {
 
   template <typename DataType, size_t order, size_t numComp>
   // Ordering of field indexing changed to 'entity, node, component'
-  auto CreateLagrangeField() {
+  auto CreateLagrangeField() const {
     return MeshField::CreateLagrangeField<ExecutionSpace, Controller, DataType,
                                           order, dim, numComp>(meshInfo);
   }
@@ -307,7 +307,7 @@ class OmegahMeshField {
   auto getCoordField() { return coordField; }
 
   // FIXME support 2d and 3d and fields with order>1
-  template <typename Field> void writeVtk(Field &field) {
+  template <typename Field> void writeVtk(Field &field) const {
     using FieldDataType = typename decltype(field.vtxField)::BaseType;
     // HACK assumes there is a vertex field.. in the Field Mixin object
     auto field_view = field.vtxField.serialize();
@@ -318,7 +318,7 @@ class OmegahMeshField {
   }
 
   template <typename ViewType = Kokkos::View<MeshField::LO *>>
-  ViewType createOffsets(size_t numTri, size_t numPtsPerElem) {
+  ViewType createOffsets(size_t numTri, size_t numPtsPerElem) const {
     ViewType offsets("offsets", numTri + 1);
     Kokkos::parallel_for(
         "setOffsets", numTri,
@@ -330,8 +330,8 @@ class OmegahMeshField {
 
   // evaluate a field at the specified local coordinate for each triangle
   template <typename ViewType, typename ShapeField>
-  auto triangleLocalPointEval(ViewType localCoords, size_t NumPtsPerElem,
-                              ShapeField field) {
+  auto triangleLocalPointEval(const ViewType &localCoords, size_t NumPtsPerElem,
+                              const ShapeField &field) const {
     auto offsets = createOffsets(meshInfo.numTri, NumPtsPerElem);
     auto eval = triangleLocalPointEval<ViewType, ShapeField>(localCoords,
                                                              offsets, field);
@@ -340,8 +340,9 @@ class OmegahMeshField {
 
   // evaluate a field at the specified local coordinates for each triangle
   template <typename ViewType, typename ShapeField>
-  auto triangleLocalPointEval(ViewType localCoords, Kokkos::View<LO *> offsets,
-                              ShapeField field) {
+  auto triangleLocalPointEval(const ViewType &localCoords,
+                              Kokkos::View<LO *> offsets,
+                              const ShapeField &field) const {
     const auto MeshDim = 2;
     if (mesh.dim() != MeshDim) {
       MeshField::fail("input mesh must be 2d\n");
@@ -360,16 +361,18 @@ class OmegahMeshField {
   }
 
   template <typename ViewType, typename ShapeField>
-  auto tetrahedronLocalPointEval(ViewType localCoords, size_t NumPtsPerElem,
-                                 ShapeField field) {
+  auto tetrahedronLocalPointEval(const ViewType &localCoords,
+                                 size_t NumPtsPerElem,
+                                 const ShapeField &field) const {
     auto offsets = createOffsets(meshInfo.numTet, NumPtsPerElem);
     auto eval = tetrahedronLocalPointEval(localCoords, offsets, field);
     return eval;
   }
 
   template <typename ViewType, typename ShapeField>
-  auto tetrahedronLocalPointEval(ViewType localCoords,
-                                 Kokkos::View<LO *> offsets, ShapeField field) {
+  auto tetrahedronLocalPointEval(const ViewType &localCoords,
+                                 Kokkos::View<LO *> offsets,
+                                 const ShapeField &field) const {
     const auto MeshDim = 3;
     if (mesh.dim() != MeshDim) {
       MeshField::fail("input mesh must be 3d\n");
diff --git a/src/MeshField_Element.hpp b/src/MeshField_Element.hpp
index d7612bf..dfe302f 100644
--- a/src/MeshField_Element.hpp
+++ b/src/MeshField_Element.hpp
@@ -236,7 +236,7 @@ struct FieldElement {
    * heavily based on SCOREC/core @ 7cd76473 apf/apfVectorElement.cc
    */
   template <typename Matrices>
-  Kokkos::View<Real *> getJacobianDeterminants(Matrices const &J) {
+  Kokkos::View<Real *> getJacobianDeterminants(Matrices const &J) const {
     static_assert(has_static_rank<Matrices>::value,
                   "Matrices must have a static rank() method.");
     static_assert(has_extent_method<Matrices>::value,
@@ -321,7 +321,7 @@ struct FieldElement {
    * @return Kokkos::View containing the jacobian for all the mesh elements
    */
   Kokkos::View<Real ***> getJacobians(Kokkos::View<Real **> localCoords,
-                                      Kokkos::View<LO *> offsets) {
+                                      Kokkos::View<LO *> offsets) const {
     if (Debug) {
       // check input parametric coords are positive and sum to one
       LO numErrors = 0;