diff --git a/NumLib/Function/Interpolation.h b/NumLib/Function/Interpolation.h
new file mode 100644
index 0000000000000000000000000000000000000000..627cca53ec68f076a1aa45793913953f77bb468e
--- /dev/null
+++ b/NumLib/Function/Interpolation.h
@@ -0,0 +1,53 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef NUMLIB_INTERPOLATION_H
+#define NUMLIB_INTERPOLATION_H
+
+#include<array>
+#include<cassert>
+
+namespace NumLib
+{
+
+/**
+ * Interpolates variables given at element nodes according to the given shape matrix.
+ *
+ * This function simply does the usual finite-element interpolation, i.e. multiplication
+ * of nodal values with the shape function.
+ *
+ * @param nodal_values   vector of nodal values, ordered by component
+ * @param shape_matrix_N shape matrix of the point to which will be interpolated
+ * @param interpolated_values array of addresses to which the interpolated values will be written
+ */
+template<typename NodalValues, typename ShapeMatrix, std::size_t NodalDOF>
+void shapeFunctionInterpolate(
+        const NodalValues& nodal_values,
+        const ShapeMatrix& shape_matrix_N,
+        std::array<double*, NodalDOF> interpolated_values
+        )
+{
+    auto const num_nodes = shape_matrix_N.size();
+
+    assert(num_nodes*NodalDOF == nodal_values.size());
+
+    for (auto d=decltype(NodalDOF){0}; d<NodalDOF; ++d)
+    {
+        *interpolated_values[d] = 0.0;
+
+        for (auto n=decltype(num_nodes){0}; n<num_nodes; ++n)
+        {
+            *interpolated_values[d] += nodal_values[d*num_nodes+n] * shape_matrix_N[n];
+        }
+    }
+}
+
+}
+
+#endif // NUMLIB_INTERPOLATION_H
diff --git a/Tests/NumLib/TestFunctionInterpolation.cpp b/Tests/NumLib/TestFunctionInterpolation.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..dd6b3a353acb7e57df75e39be0e366ab88bdfcac
--- /dev/null
+++ b/Tests/NumLib/TestFunctionInterpolation.cpp
@@ -0,0 +1,98 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/LICENSE.txt
+ */
+
+#include <gtest/gtest.h>
+
+#include <array>
+
+#include "NumLib/Function/Interpolation.h"
+
+#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+#include "NumLib/Fem/ShapeFunction/ShapeLine2.h"
+
+#include "NumLib/Fem/ShapeMatrixPolicy.h"
+
+#include "MeshLib/Elements/Element.h"
+
+#include "NumLib/Fem/Integration/GaussIntegrationPolicy.h"
+
+
+TEST(NumLibFunctionInterpolationTest, TwoVariablesTwoNodes)
+{
+    double variable1 = 0.0;
+    double variable2 = 0.0;
+    std::array<double*, 2> interpolated_values = { &variable1, &variable2 };
+
+    const std::array<double, 4> nodal_values = {
+        0.0, 1.0, // for variable1
+        -1.0, 1.0  // for variable2
+    };
+
+    const std::array<double, 2> shape_matrix = { 0.25, 0.75 };
+
+    NumLib::shapeFunctionInterpolate(nodal_values, shape_matrix,
+                                     interpolated_values);
+
+    ASSERT_EQ(0.75, variable1);
+    ASSERT_EQ(0.5,  variable2);
+}
+
+
+TEST(NumLibFunctionInterpolationTest, Linear1DElement)
+{
+    // typedefs
+    using ShapeFunction = NumLib::ShapeLine2;
+    using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, 1u>;
+
+    using FemType = NumLib::TemplateIsoparametric<ShapeFunction, ShapeMatricesType>;
+
+    using IntegrationMethod = NumLib::GaussIntegrationPolicy<
+        ShapeFunction::MeshElement>::IntegrationMethod;
+
+    // set up mesh element
+    auto pt_a = MeshLib::Node({ 0.0, 0.0, 0.0 }, 0);
+    auto pt_b = MeshLib::Node({ 1.0, 0.0, 0.0 }, 0);
+
+    auto const element = MeshLib::Line(std::array<MeshLib::Node*, 2>{
+                                           &pt_a, &pt_b
+                                       }, 0u, 0u);
+
+    // set up shape function
+    FemType finite_element(*static_cast<const ShapeFunction::MeshElement*>(&element));
+
+    const unsigned integration_order = 2;
+    IntegrationMethod integration_method(integration_order);
+
+    ShapeMatricesType::ShapeMatrices shape_matrix;
+
+    finite_element.computeShapeFunctions(
+            integration_method.getWeightedPoint(0).getCoords(),
+            shape_matrix);
+    ASSERT_EQ(2, shape_matrix.N.size());
+
+    // actual test
+    double variable1 = 0.0;
+    double variable2 = 0.0;
+    std::array<double*, 2> interpolated_values = { &variable1, &variable2 };
+
+    const std::array<double, 4> nodal_values = {
+        0.0, 1.0, // for variable1
+        -1.0, 1.0  // for variable2
+    };
+
+    NumLib::shapeFunctionInterpolate(nodal_values, shape_matrix.N,
+                                     interpolated_values);
+
+    const double n0 = shape_matrix.N[0];
+    const double n1 = shape_matrix.N[1];
+
+    ASSERT_EQ( 0.0 * n0 + 1.0 * n1, variable1);
+    ASSERT_EQ(-1.0 * n0 + 1.0 * n1, variable2);
+}
+
+