diff --git a/src/algebra/IntegrationTools.hpp b/src/algebra/IntegrationTools.hpp
index 99676b0480451ea883785bcac6d11c59c725f2e1..bdf6de6929f194696b0f467ea7af45ded1b52e62 100644
--- a/src/algebra/IntegrationTools.hpp
+++ b/src/algebra/IntegrationTools.hpp
@@ -9,7 +9,6 @@
 
 #include <cmath>
 #include <language/utils/EvaluateAtPoints.hpp>
-#include <language/utils/PugsFunctionAdapter.hpp>
 #include <utils/Types.hpp>
 
 enum class QuadratureType : int8_t
@@ -358,10 +357,10 @@ class IntegrationTools
   PUGS_INLINE T
   integrate(const FunctionSymbolId& function, const double& a, const double& b) const
   {
-    T result                       = 0;
-    Array<TinyVector<1>> positions = quadraturePoints(a, b);
-    Array<double> weights          = m_integration_method->weights();
-    Array<T> values                = EvaluateAtPoints<double(TinyVector<1>)>::evaluate(function, positions);
+    T result                             = 0;
+    Array<const TinyVector<1>> positions = quadraturePoints(a, b);
+    Array<double> weights                = m_integration_method->weights();
+    Array<T> values                      = EvaluateAtPoints<double(const TinyVector<1>)>::evaluate(function, positions);
     Assert(positions.size() == weights.size(), "Wrong number of quadrature points and/or weights in Gauss quadrature");
     for (size_t i = 0; i < values.size(); ++i) {
       result += weights[i] * values[i];
@@ -372,11 +371,11 @@ class IntegrationTools
   PUGS_INLINE Array<T>
   integrateFunction(const FunctionSymbolId& function, const Array<TinyVector<1>>& vertices) const
   {
-    Array<TinyVector<1>> positions = quadraturePoints(vertices);
+    Array<const TinyVector<1>> positions = quadraturePoints(vertices);
     Array<T> result{vertices.size() - 1};
     Array<double> interval_size{vertices.size() - 1};
     Array<double> weights = m_integration_method->weights();
-    Array<T> values       = EvaluateAtPoints<double(TinyVector<1>)>::evaluate(function, positions);
+    Array<T> values       = EvaluateAtPoints<double(const TinyVector<1>)>::evaluate(function, positions);
 
     for (size_t i = 0; i < interval_size.size(); ++i) {
       interval_size[i] = vertices[i + 1][0] - vertices[i][0];
diff --git a/src/language/modules/MeshModule.cpp b/src/language/modules/MeshModule.cpp
index 61767078aeaa0f8362393196a01ce9cde167a064..d06d378dd8eb5a79a80f41497cb592407db37849 100644
--- a/src/language/modules/MeshModule.cpp
+++ b/src/language/modules/MeshModule.cpp
@@ -4,7 +4,6 @@
 #include <language/node_processor/ExecutionPolicy.hpp>
 #include <language/utils/BuiltinFunctionEmbedder.hpp>
 #include <language/utils/FunctionTable.hpp>
-#include <language/utils/PugsFunctionAdapter.hpp>
 #include <language/utils/SymbolTable.hpp>
 #include <language/utils/TypeDescriptor.hpp>
 #include <mesh/CartesianMeshBuilder.hpp>
diff --git a/src/language/utils/EvaluateAtPoints.hpp b/src/language/utils/EvaluateAtPoints.hpp
index 2ba01887d79575e7b320b0d1187cb5995ba16350..d842dd44b96b02caf23515d4c35da0a0ee277865 100644
--- a/src/language/utils/EvaluateAtPoints.hpp
+++ b/src/language/utils/EvaluateAtPoints.hpp
@@ -2,31 +2,36 @@
 #define EVALUATE_AT_POINTS_HPP
 
 #include <language/utils/PugsFunctionAdapter.hpp>
+#include <utils/Array.hpp>
+
+class FunctionSymbolId;
 
 template <typename T>
 class EvaluateAtPoints;
 template <typename OutputType, typename InputType>
 class EvaluateAtPoints<OutputType(InputType)> : public PugsFunctionAdapter<OutputType(InputType)>
 {
-  //  static constexpr size_t Dimension = OutputType::Dimension;
   using Adapter = PugsFunctionAdapter<OutputType(InputType)>;
 
  public:
-  //  template <size_t vectorlength>
-  static inline Array<OutputType>
-  evaluate(const FunctionSymbolId& function_symbol_id, const Array<const InputType>& position)
+  template <typename InputArrayT, typename OutputArrayT>
+  static PUGS_INLINE void
+  evaluateTo(const FunctionSymbolId& function_symbol_id, const InputArrayT& position, OutputArrayT& value)
   {
-    //    static_assert(position.size()==vectorlength, "The length of the output must be equal to the template
-    //    argument");
+    static_assert(std::is_same_v<std::remove_const_t<typename InputArrayT::data_type>, InputType>,
+                  "invalid input data type");
+    static_assert(std::is_same_v<std::remove_const_t<typename OutputArrayT::data_type>, OutputType>,
+                  "invalid output data type");
+    Assert(size(value) == size(position));
+
     auto& expression    = Adapter::getFunctionExpression(function_symbol_id);
     auto convert_result = Adapter::getResultConverter(expression.m_data_type);
 
-    Array<ExecutionPolicy> context_list = Adapter::getContextList(expression);
+    auto context_list = Adapter::getContextList(expression);
 
     using execution_space = typename Kokkos::DefaultExecutionSpace::execution_space;
     Kokkos::Experimental::UniqueToken<execution_space, Kokkos::Experimental::UniqueTokenScope::Global> tokens;
 
-    Array<OutputType> value(position.size());
     if constexpr (std::is_arithmetic_v<OutputType>) {
       value.fill(0);
     } else if constexpr (is_tiny_vector_v<OutputType> or is_tiny_matrix_v<OutputType>) {
@@ -35,7 +40,7 @@ class EvaluateAtPoints<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
       static_assert(std::is_same_v<OutputType, double>, "unexpected output type");
     }
 
-    parallel_for(position.size(), [=, &expression, &tokens](size_t i) {
+    parallel_for(size(position), [=, &expression, &tokens](typename InputArrayT::index_type i) {
       const int32_t t = tokens.acquire();
 
       auto& execution_policy = context_list[t];
@@ -46,8 +51,18 @@ class EvaluateAtPoints<OutputType(InputType)> : public PugsFunctionAdapter<Outpu
 
       tokens.release(t);
     });
+  }
 
+  template <class InputArrayT>
+  static PUGS_INLINE Array<OutputType>
+  evaluate(const FunctionSymbolId& function_symbol_id, const InputArrayT& position)
+  {
+    static_assert(std::is_same_v<std::remove_const_t<typename InputArrayT::data_type>, InputType>,
+                  "invalid input data type");
+    Array<OutputType> value(size(position));
+    evaluateTo(function_symbol_id, position, value);
     return value;
   }
 };
-#endif   // INTERPOLATE_ITEM_VALUE_HPP
+
+#endif   // EVALUATE_AT_POINTS_HPP
diff --git a/src/language/utils/InterpolateItemArray.hpp b/src/language/utils/InterpolateItemArray.hpp
index c507c527ca9f1c59070410c09f802f524a3de083..17b72de6d62e0ac68062343010154257c43d3006 100644
--- a/src/language/utils/InterpolateItemArray.hpp
+++ b/src/language/utils/InterpolateItemArray.hpp
@@ -2,17 +2,15 @@
 #define INTERPOLATE_ITEM_ARRAY_HPP
 
 #include <language/utils/InterpolateItemValue.hpp>
-#include <language/utils/PugsFunctionAdapter.hpp>
 #include <mesh/ItemArray.hpp>
 #include <mesh/ItemType.hpp>
 
 template <typename T>
 class InterpolateItemArray;
 template <typename OutputType, typename InputType>
-class InterpolateItemArray<OutputType(InputType)> : public PugsFunctionAdapter<OutputType(InputType)>
+class InterpolateItemArray<OutputType(InputType)>
 {
   static constexpr size_t Dimension = OutputType::Dimension;
-  using Adapter                     = PugsFunctionAdapter<OutputType(InputType)>;
 
  public:
   template <ItemType item_type>
diff --git a/src/language/utils/InterpolateItemValue.hpp b/src/language/utils/InterpolateItemValue.hpp
index de12b563a19ded640f0b4bbae3184f383d90a2f3..dc3ef056b0693d4d32e695a33b7772aa11264b73 100644
--- a/src/language/utils/InterpolateItemValue.hpp
+++ b/src/language/utils/InterpolateItemValue.hpp
@@ -1,46 +1,23 @@
 #ifndef INTERPOLATE_ITEM_VALUE_HPP
 #define INTERPOLATE_ITEM_VALUE_HPP
 
-#include <language/utils/PugsFunctionAdapter.hpp>
+#include <language/utils/EvaluateAtPoints.hpp>
 #include <mesh/ItemType.hpp>
 #include <mesh/ItemValue.hpp>
 
 template <typename T>
 class InterpolateItemValue;
 template <typename OutputType, typename InputType>
-class InterpolateItemValue<OutputType(InputType)> : public PugsFunctionAdapter<OutputType(InputType)>
+class InterpolateItemValue<OutputType(InputType)>
 {
-  static constexpr size_t Dimension = OutputType::Dimension;
-  using Adapter                     = PugsFunctionAdapter<OutputType(InputType)>;
-
  public:
   template <ItemType item_type>
   PUGS_INLINE static ItemValue<OutputType, item_type>
   interpolate(const FunctionSymbolId& function_symbol_id, const ItemValue<const InputType, item_type>& position)
   {
-    auto& expression    = Adapter::getFunctionExpression(function_symbol_id);
-    auto convert_result = Adapter::getResultConverter(expression.m_data_type);
-
-    Array<ExecutionPolicy> context_list = Adapter::getContextList(expression);
-
-    using execution_space = typename Kokkos::DefaultExecutionSpace::execution_space;
-    Kokkos::Experimental::UniqueToken<execution_space, Kokkos::Experimental::UniqueTokenScope::Global> tokens;
     const IConnectivity& connectivity = *position.connectivity_ptr();
-
     ItemValue<OutputType, item_type> value(connectivity);
-    using ItemId = ItemIdT<item_type>;
-
-    parallel_for(connectivity.template numberOf<item_type>(), [=, &expression, &tokens](ItemId i) {
-      const int32_t t = tokens.acquire();
-
-      auto& execution_policy = context_list[t];
-
-      Adapter::convertArgs(execution_policy.currentContext(), position[i]);
-      auto result = expression.execute(execution_policy);
-      value[i]    = convert_result(std::move(result));
-
-      tokens.release(t);
-    });
+    EvaluateAtPoints<OutputType(const InputType)>::evaluateTo(function_symbol_id, position, value);
 
     return value;
   }
@@ -51,31 +28,15 @@ class InterpolateItemValue<OutputType(InputType)> : public PugsFunctionAdapter<O
               const ItemValue<const InputType, item_type>& position,
               const Array<const ItemIdT<item_type>>& list_of_items)
   {
-    auto& expression    = Adapter::getFunctionExpression(function_symbol_id);
-    auto convert_result = Adapter::getResultConverter(expression.m_data_type);
-
-    Array<ExecutionPolicy> context_list = Adapter::getContextList(expression);
-
-    using execution_space = typename Kokkos::DefaultExecutionSpace::execution_space;
-    Kokkos::Experimental::UniqueToken<execution_space, Kokkos::Experimental::UniqueTokenScope::Global> tokens;
-
-    Array<OutputType> value{list_of_items.size()};
+    Array<InputType> item_position{list_of_items.size()};
     using ItemId = ItemIdT<item_type>;
+    parallel_for(
+      list_of_items.size(), PUGS_LAMBDA(size_t i_item) {
+        ItemId item_id        = list_of_items[i_item];
+        item_position[i_item] = position[item_id];
+      });
 
-    parallel_for(list_of_items.size(), [=, &expression, &tokens](size_t i_item) {
-      ItemId item_id  = list_of_items[i_item];
-      const int32_t t = tokens.acquire();
-
-      auto& execution_policy = context_list[t];
-
-      Adapter::convertArgs(execution_policy.currentContext(), position[item_id]);
-      auto result   = expression.execute(execution_policy);
-      value[i_item] = convert_result(std::move(result));
-
-      tokens.release(t);
-    });
-
-    return value;
+    return EvaluateAtPoints<OutputType(const InputType)>::evaluate(function_symbol_id, item_position);
   }
 
   template <ItemType item_type>
diff --git a/src/mesh/MeshTransformer.cpp b/src/mesh/MeshTransformer.cpp
index 67b3688926dcfb765c44799cff4674955515eb8b..167c4f337cd5c48147594623f3ae436874a47a22 100644
--- a/src/mesh/MeshTransformer.cpp
+++ b/src/mesh/MeshTransformer.cpp
@@ -3,15 +3,12 @@
 #include <mesh/Connectivity.hpp>
 #include <mesh/Mesh.hpp>
 
-#include <language/utils/FunctionTable.hpp>
-#include <language/utils/PugsFunctionAdapter.hpp>
-#include <language/utils/SymbolTable.hpp>
+#include <language/utils/EvaluateAtPoints.hpp>
 
 template <typename OutputType, typename InputType>
-class MeshTransformer::MeshTransformation<OutputType(InputType)> : public PugsFunctionAdapter<OutputType(InputType)>
+class MeshTransformer::MeshTransformation<OutputType(InputType)>
 {
   static constexpr size_t Dimension = OutputType::Dimension;
-  using Adapter                     = PugsFunctionAdapter<OutputType(InputType)>;
 
  public:
   static inline std::shared_ptr<Mesh<Connectivity<Dimension>>>
@@ -20,28 +17,9 @@ class MeshTransformer::MeshTransformation<OutputType(InputType)> : public PugsFu
     using MeshType             = Mesh<Connectivity<Dimension>>;
     const MeshType& given_mesh = dynamic_cast<const MeshType&>(*p_mesh);
 
-    auto& expression    = Adapter::getFunctionExpression(function_symbol_id);
-    auto convert_result = Adapter::getResultConverter(expression.m_data_type);
-
-    Array<ExecutionPolicy> context_list = Adapter::getContextList(expression);
-
-    NodeValue<const InputType> given_xr = given_mesh.xr();
     NodeValue<OutputType> xr(given_mesh.connectivity());
-
-    using execution_space = typename Kokkos::DefaultExecutionSpace::execution_space;
-    Kokkos::Experimental::UniqueToken<execution_space, Kokkos::Experimental::UniqueTokenScope::Global> tokens;
-
-    parallel_for(given_mesh.numberOfNodes(), [=, &expression, &tokens](NodeId r) {
-      const int32_t t = tokens.acquire();
-
-      auto& execution_policy = context_list[t];
-
-      Adapter::convertArgs(execution_policy.currentContext(), given_xr[r]);
-      auto result = expression.execute(execution_policy);
-      xr[r]       = convert_result(std::move(result));
-
-      tokens.release(t);
-    });
+    NodeValue<const InputType> given_xr = given_mesh.xr();
+    EvaluateAtPoints<OutputType(InputType)>::evaluateTo(function_symbol_id, given_xr, xr);
 
     return std::make_shared<MeshType>(given_mesh.shared_connectivity(), xr);
   }
diff --git a/tests/test_PugsFunctionAdapter.cpp b/tests/test_PugsFunctionAdapter.cpp
index 76f3a06aa3f96c208e6c02c16f8006e22160a8f6..cb0dbc7c5e1b76dfeeb6390ddca70524e22d7dc1 100644
--- a/tests/test_PugsFunctionAdapter.cpp
+++ b/tests/test_PugsFunctionAdapter.cpp
@@ -33,7 +33,7 @@ class TestBinary<OutputType(InputType...)> : public PugsFunctionAdapter<OutputTy
     auto& expression    = Adapter::getFunctionExpression(function_symbol_id);
     auto convert_result = Adapter::getResultConverter(expression.m_data_type);
 
-    Array<ExecutionPolicy> context_list = Adapter::getContextList(expression);
+    auto context_list = Adapter::getContextList(expression);
 
     auto& execution_policy = context_list[0];
 
@@ -50,7 +50,7 @@ class TestBinary<OutputType(InputType...)> : public PugsFunctionAdapter<OutputTy
     auto& expression    = Adapter::getFunctionExpression(function_symbol_id);
     auto convert_result = Adapter::getResultConverter(expression.m_data_type);
 
-    Array<ExecutionPolicy> context_list = Adapter::getContextList(expression);
+    auto context_list = Adapter::getContextList(expression);
 
     auto& execution_policy = context_list[0];