diff --git a/NumLib/CMakeLists.txt b/NumLib/CMakeLists.txt
index 1f4c8a0c78f8ca256bce04de27a98ae15fcd231a..18a164b0696b70ae0d6464c9e45c41c81907a4e6 100644
--- a/NumLib/CMakeLists.txt
+++ b/NumLib/CMakeLists.txt
@@ -19,6 +19,8 @@ SET ( SOURCES ${SOURCES} ${SOURCES_TIMESTEP} ${SOURCES_TIMESTEP_ALGORITHMS})
 
 GET_SOURCE_FILES(SOURCES_FUNCTION Function)
 SET ( SOURCES ${SOURCES} ${SOURCES_FUNCTION})
+GET_SOURCE_FILES(SOURCES_DISTRIBUTION Distribution)
+SET ( SOURCES ${SOURCES} ${SOURCES_DISTRIBUTION})
 
 
 INCLUDE_DIRECTORIES (
diff --git a/NumLib/Distribution/Distribution.cpp b/NumLib/Distribution/Distribution.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..bce94991f1f988ed90acfdb606340b0c9d88895c
--- /dev/null
+++ b/NumLib/Distribution/Distribution.cpp
@@ -0,0 +1,37 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2014, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+
+#include "Distribution.h"
+
+#include <algorithm>
+
+#include "MeshLib/Mesh.h"
+#include "MeshLib/Node.h"
+#include "NumLib/Function/ISpatialFunction.h"
+
+namespace NumLib
+{
+
+std::vector<double> generateNodeValueDistribution(
+	const NumLib::ISpatialFunction &func,
+	const MeshLib::Mesh &msh,
+	const std::vector<std::size_t> &vec_node_ids)
+{
+	// evaluate a given function with nodal coordinates
+	std::vector<double> vec_values(vec_node_ids.size());
+	std::transform(vec_node_ids.begin(), vec_node_ids.end(), vec_values.begin(),
+		[&func, &msh](std::size_t node_id)
+		{
+			return func(*msh.getNode(node_id));
+		});
+	return vec_values;
+}
+
+} //NumLib
diff --git a/NumLib/Distribution/Distribution.h b/NumLib/Distribution/Distribution.h
new file mode 100644
index 0000000000000000000000000000000000000000..ce64fb0107990e5803cc40774ef5d41242e5b332
--- /dev/null
+++ b/NumLib/Distribution/Distribution.h
@@ -0,0 +1,40 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2014, 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 DISTRIBUTION_H_
+#define DISTRIBUTION_H_
+
+#include <vector>
+
+namespace MeshLib
+{
+class Mesh;
+}
+
+namespace NumLib
+{
+class ISpatialFunction;
+
+/**
+ * Generate distributed node values from a given function
+ *
+ * @param func            a spatial function object
+ * @param msh             a mesh object
+ * @param vec_node_ids    a vector of mesh node ids where the function is evaluated
+ * @return a vector of nodal values
+ */
+std::vector<double> generateNodeValueDistribution(
+    const NumLib::ISpatialFunction &func,
+    const MeshLib::Mesh &msh,
+    const std::vector<std::size_t> &vec_node_ids);
+
+} // NumLib
+
+#endif // DISTRIBUTION_H_
diff --git a/Tests/NumLib/TestDistribution.cpp b/Tests/NumLib/TestDistribution.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..fbb037164b3d75b8af93a4f863118cc6ab4f606e
--- /dev/null
+++ b/Tests/NumLib/TestDistribution.cpp
@@ -0,0 +1,225 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2014, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include <gtest/gtest.h>
+
+#include <vector>
+#include <memory>
+
+#include "GEOObjects.h"
+
+#include "MeshLib/Mesh.h"
+#include "MeshLib/MeshGenerators/MeshGenerator.h"
+#include "MeshGeoToolsLib/MeshNodeSearcher.h"
+#include "MeshGeoToolsLib/MeshNodesAlongPolyline.h"
+#include "MeshGeoToolsLib/MeshNodesAlongSurface.h"
+#include "NumLib/Function/LinearInterpolationAlongPolyline.h"
+#include "NumLib/Function/LinearInterpolationOnSurface.h"
+#include "NumLib/Function/SpatialFunctionLinear.h"
+#include "NumLib/Distribution/Distribution.h"
+
+#include "../TestTools.h"
+
+class NumLibDistributionQuad : public testing::Test
+{
+public:
+	NumLibDistributionQuad() :
+		_geometric_size(10.0), _number_of_subdivisions_per_direction(10),
+		_msh(MeshLib::MeshGenerator::generateRegularQuadMesh(_geometric_size, _number_of_subdivisions_per_direction)),
+		_project_name("test"), _mshNodesSearcher(*_msh), _ply0(nullptr)
+	{
+		// create geometry
+		std::vector<GeoLib::Point*>* pnts (new std::vector<GeoLib::Point*>);
+		pnts->push_back(new GeoLib::Point(0.0, 0.0, 0.0));
+		pnts->push_back(new GeoLib::Point(_geometric_size, 0.0, 0.0));
+		pnts->push_back(new GeoLib::Point(_geometric_size, _geometric_size, 0.0));
+		pnts->push_back(new GeoLib::Point(0.0, _geometric_size, 0.0));
+
+		std::vector<GeoLib::Polyline*>* plys (new std::vector<GeoLib::Polyline*>);
+		_ply0 = new GeoLib::Polyline(*pnts);
+		_ply0->addPoint(0);
+		_ply0->addPoint(1);
+		plys->push_back(_ply0);
+
+		GeoLib::Polyline* ply1 = new GeoLib::Polyline(*pnts);
+		ply1->addPoint(0);
+		ply1->addPoint(1);
+		ply1->addPoint(2);
+		ply1->addPoint(3);
+		ply1->addPoint(0);
+		plys->push_back(ply1);
+
+		std::vector<GeoLib::Surface*>* sfcs (new std::vector<GeoLib::Surface*>);
+		_sfc1 = GeoLib::Surface::createSurface(*ply1);
+		sfcs->push_back(_sfc1);
+
+		_geo_objs.addPointVec(pnts,_project_name);
+		_geo_objs.addPolylineVec(plys, _project_name);
+		_geo_objs.addSurfaceVec(sfcs, _project_name);
+	}
+
+protected:
+	const double _geometric_size;
+	const std::size_t _number_of_subdivisions_per_direction;
+	std::unique_ptr<MeshLib::Mesh> _msh;
+	GeoLib::GEOObjects _geo_objs;
+	std::string _project_name;
+	MeshGeoToolsLib::MeshNodeSearcher _mshNodesSearcher;
+	GeoLib::Polyline* _ply0;
+	GeoLib::Surface* _sfc1;
+};
+
+class NumLibDistributionHex : public testing::Test
+{
+public:
+	NumLibDistributionHex() :
+		_geometric_size(10.0), _number_of_subdivisions_per_direction(10),
+		_msh(MeshLib::MeshGenerator::generateRegularHexMesh(_geometric_size, _number_of_subdivisions_per_direction)),
+		_project_name("test"), _mshNodesSearcher(*_msh), _ply0(nullptr)
+	{
+		// create geometry
+		std::vector<GeoLib::Point*>* pnts (new std::vector<GeoLib::Point*>);
+		pnts->push_back(new GeoLib::Point(0.0, 0.0, 0.0));
+		pnts->push_back(new GeoLib::Point(_geometric_size, 0.0, 0.0));
+		pnts->push_back(new GeoLib::Point(_geometric_size, _geometric_size, 0.0));
+		pnts->push_back(new GeoLib::Point(0.0, _geometric_size, 0.0));
+		pnts->push_back(new GeoLib::Point(0.0, 0.0, _geometric_size));
+		pnts->push_back(new GeoLib::Point(_geometric_size, 0.0, _geometric_size));
+		pnts->push_back(new GeoLib::Point(_geometric_size, _geometric_size, _geometric_size));
+		pnts->push_back(new GeoLib::Point(0.0, _geometric_size, _geometric_size));
+
+		std::vector<GeoLib::Polyline*>* plys (new std::vector<GeoLib::Polyline*>);
+		_ply0 = new GeoLib::Polyline(*pnts); // vertical polyline
+		_ply0->addPoint(0);
+		_ply0->addPoint(4);
+		plys->push_back(_ply0);
+		GeoLib::Polyline* ply1 = new GeoLib::Polyline(*pnts); // polygon for left surface
+		ply1->addPoint(0);
+		ply1->addPoint(3);
+		ply1->addPoint(7);
+		ply1->addPoint(4);
+		ply1->addPoint(0);
+		plys->push_back(ply1);
+
+		std::vector<GeoLib::Surface*>* sfcs (new std::vector<GeoLib::Surface*>);
+		_sfc1 = GeoLib::Surface::createSurface(*ply1);
+		sfcs->push_back(_sfc1);
+
+		_geo_objs.addPointVec(pnts,_project_name);
+		_geo_objs.addPolylineVec(plys, _project_name);
+		_geo_objs.addSurfaceVec(sfcs, _project_name);
+	}
+
+protected:
+	const double _geometric_size;
+	const std::size_t _number_of_subdivisions_per_direction;
+	std::unique_ptr<MeshLib::Mesh> _msh;
+	GeoLib::GEOObjects _geo_objs;
+	std::string _project_name;
+	MeshGeoToolsLib::MeshNodeSearcher _mshNodesSearcher;
+	GeoLib::Polyline* _ply0;
+	GeoLib::Surface* _sfc1;
+};
+
+TEST_F(NumLibDistributionQuad, Linear)
+{
+	// f(x,y,z) = 1 + 2x + 3y + 4z
+	std::array<double,4> f_coeff = {{1, 2, 3, 4}};
+	MathLib::LinearFunction<double,3> linear_f(f_coeff);
+	std::vector<double> expected = {{1, 3, 5, 7, 9, 11, 13, 15, 17, 19, 21}};
+
+	NumLib::SpatialFunctionLinear f(linear_f);
+	const std::vector<std::size_t>& vec_node_ids = _mshNodesSearcher.getMeshNodesAlongPolyline(*_ply0).getNodeIDs();
+	std::vector<double> interpolated_values = NumLib::generateNodeValueDistribution(f, *_msh, vec_node_ids);
+	ASSERT_ARRAY_NEAR(expected, interpolated_values, expected.size(), std::numeric_limits<double>::epsilon());
+}
+
+TEST_F(NumLibDistributionHex, Linear)
+{
+	// f(x,y,z) = 1 + 2x + 3y + 4z
+	std::array<double,4> f_coeff = {{1, 2, 3, 4}};
+	MathLib::LinearFunction<double,3> linear_f(f_coeff);
+	std::vector<double> expected(std::pow(_number_of_subdivisions_per_direction+1, 2));
+	const double dL = _geometric_size / _number_of_subdivisions_per_direction;
+	for (std::size_t i=0; i<expected.size(); i++) {
+		double x = 0;
+		double y = (i%(_number_of_subdivisions_per_direction+1)) * dL;
+		double z = (i/(_number_of_subdivisions_per_direction+1)) * dL;
+		expected[i] = f_coeff[0] + f_coeff[1]*x + f_coeff[2]*y + f_coeff[3]*z;
+	}
+
+	NumLib::SpatialFunctionLinear f(linear_f);
+	const std::vector<std::size_t>& vec_node_ids = _mshNodesSearcher.getMeshNodesAlongSurface(*_sfc1).getNodeIDs();
+	std::vector<double> interpolated_values = NumLib::generateNodeValueDistribution(f, *_msh, vec_node_ids);
+
+	ASSERT_ARRAY_NEAR(expected, interpolated_values, expected.size(), std::numeric_limits<double>::epsilon());
+}
+
+TEST_F(NumLibDistributionQuad, InterpolationPolyline)
+{
+	const std::vector<std::size_t> vec_point_ids = {{0, 1}};
+	const std::vector<double> vec_point_values = {{0., 100.}};
+	std::vector<double> expected = {{0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100}};
+
+	NumLib::LinearInterpolationAlongPolyline interpolate(
+			*_ply0, vec_point_ids, vec_point_values, std::numeric_limits<double>::epsilon(), std::numeric_limits<double>::max());
+	const std::vector<std::size_t>& vec_node_ids = _mshNodesSearcher.getMeshNodesAlongPolyline(*_ply0).getNodeIDs();
+	std::vector<double> interpolated_values = NumLib::generateNodeValueDistribution(interpolate, *_msh, vec_node_ids);
+
+	ASSERT_ARRAY_NEAR(expected, interpolated_values, expected.size(), std::numeric_limits<double>::epsilon());
+}
+
+TEST_F(NumLibDistributionHex, InterpolationPolyline)
+{
+	const std::vector<std::size_t> vec_point_ids = {{0, 4}};
+	const std::vector<double> vec_point_values = {{0., 100.}};
+	std::vector<double> expected = {{0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100}};
+
+	NumLib::LinearInterpolationAlongPolyline interpolate(
+			*_ply0, vec_point_ids, vec_point_values, std::numeric_limits<double>::epsilon(), std::numeric_limits<double>::max());
+	const std::vector<std::size_t>& vec_node_ids = _mshNodesSearcher.getMeshNodesAlongPolyline(*_ply0).getNodeIDs();
+	std::vector<double> interpolated_values = NumLib::generateNodeValueDistribution(interpolate, *_msh, vec_node_ids);
+
+	ASSERT_ARRAY_NEAR(expected, interpolated_values, expected.size(), std::numeric_limits<double>::epsilon());
+}
+
+TEST_F(NumLibDistributionQuad, InterpolationSurface)
+{
+	const std::vector<std::size_t> vec_point_ids = {{0, 1, 2, 3}};
+	const std::vector<double> vec_point_values = {{0., 100., 100., 0.}};
+	std::vector<double> expected(_msh->getNNodes());
+	for (std::size_t i=0; i<_msh->getNNodes(); i++) {
+		expected[i] = (i%(_number_of_subdivisions_per_direction+1)) * 10;
+	}
+
+	NumLib::LinearInterpolationOnSurface interpolate(*_sfc1, vec_point_ids, vec_point_values, std::numeric_limits<double>::max());
+	const std::vector<std::size_t>& vec_node_ids = _mshNodesSearcher.getMeshNodesAlongSurface(*_sfc1).getNodeIDs();
+	std::vector<double> interpolated_values = NumLib::generateNodeValueDistribution(interpolate, *_msh, vec_node_ids);
+
+	// the machine epsilon for double is too small for this test
+	ASSERT_ARRAY_NEAR(expected, interpolated_values, expected.size(), std::numeric_limits<float>::epsilon());
+}
+
+TEST_F(NumLibDistributionHex, InterpolationSurface)
+{
+	const std::vector<std::size_t> vec_point_ids = {{0, 3, 7, 4}};
+	const std::vector<double> vec_point_values = {{0., 100., 100., 0.}};
+	std::vector<double> expected(std::pow(_number_of_subdivisions_per_direction+1, 2));
+	for (std::size_t i=0; i<expected.size(); i++) {
+		expected[i] = (i%(_number_of_subdivisions_per_direction+1)) * 10;
+	}
+
+	NumLib::LinearInterpolationOnSurface interpolate(*_sfc1, vec_point_ids, vec_point_values, std::numeric_limits<double>::max());
+	const std::vector<std::size_t>& vec_node_ids = _mshNodesSearcher.getMeshNodesAlongSurface(*_sfc1).getNodeIDs();
+	std::vector<double> interpolated_values = NumLib::generateNodeValueDistribution(interpolate, *_msh, vec_node_ids);
+
+	ASSERT_ARRAY_NEAR(expected, interpolated_values, expected.size(), std::numeric_limits<float>::epsilon());
+}
+
+
diff --git a/Tests/NumLib/TestSpatialFunction.cpp b/Tests/NumLib/TestSpatialFunction.cpp
index 790b26b5e9e063cb37da36ae293753617ee45b7f..dec0aef74e609b2da4f49f461aa15f54dcbb36b5 100644
--- a/Tests/NumLib/TestSpatialFunction.cpp
+++ b/Tests/NumLib/TestSpatialFunction.cpp
@@ -16,239 +16,101 @@
 
 #include "MeshLib/Mesh.h"
 #include "MeshLib/MeshGenerators/MeshGenerator.h"
-#include "MeshGeoToolsLib/MeshNodeSearcher.h"
-#include "MeshGeoToolsLib/MeshNodesAlongPolyline.h"
-#include "MeshGeoToolsLib/MeshNodesAlongSurface.h"
 #include "NumLib/Function/LinearInterpolationAlongPolyline.h"
 #include "NumLib/Function/LinearInterpolationOnSurface.h"
 #include "NumLib/Function/SpatialFunctionLinear.h"
 
 #include "../TestTools.h"
 
-class NumLibSpatialFunctionQuad : public testing::Test
-{
-public:
-	NumLibSpatialFunctionQuad() :
-		_geometric_size(10.0), _number_of_subdivisions_per_direction(10),
-		_msh(MeshLib::MeshGenerator::generateRegularQuadMesh(_geometric_size, _number_of_subdivisions_per_direction)),
-		_project_name("test"), _mshNodesSearcher(*_msh), _ply0(nullptr)
-	{
-		// create geometry
-		std::vector<GeoLib::Point*>* pnts (new std::vector<GeoLib::Point*>);
-		pnts->push_back(new GeoLib::Point(0.0, 0.0, 0.0));
-		pnts->push_back(new GeoLib::Point(_geometric_size, 0.0, 0.0));
-		pnts->push_back(new GeoLib::Point(_geometric_size, _geometric_size, 0.0));
-		pnts->push_back(new GeoLib::Point(0.0, _geometric_size, 0.0));
-
-		std::vector<GeoLib::Polyline*>* plys (new std::vector<GeoLib::Polyline*>);
-		_ply0 = new GeoLib::Polyline(*pnts);
-		_ply0->addPoint(0);
-		_ply0->addPoint(1);
-		plys->push_back(_ply0);
-
-		GeoLib::Polyline* ply1 = new GeoLib::Polyline(*pnts);
-		ply1->addPoint(0);
-		ply1->addPoint(1);
-		ply1->addPoint(2);
-		ply1->addPoint(3);
-		ply1->addPoint(0);
-		plys->push_back(ply1);
-
-		std::vector<GeoLib::Surface*>* sfcs (new std::vector<GeoLib::Surface*>);
-		_sfc1 = GeoLib::Surface::createSurface(*ply1);
-		sfcs->push_back(_sfc1);
-
-		_geo_objs.addPointVec(pnts,_project_name);
-		_geo_objs.addPolylineVec(plys, _project_name);
-		_geo_objs.addSurfaceVec(sfcs, _project_name);
-	}
-
-protected:
-	const double _geometric_size;
-	const std::size_t _number_of_subdivisions_per_direction;
-	std::unique_ptr<MeshLib::Mesh> _msh;
-	GeoLib::GEOObjects _geo_objs;
-	std::string _project_name;
-	MeshGeoToolsLib::MeshNodeSearcher _mshNodesSearcher;
-	GeoLib::Polyline* _ply0;
-	GeoLib::Surface* _sfc1;
-};
-
-class NumLibSpatialFunctionHex : public testing::Test
-{
-public:
-	NumLibSpatialFunctionHex() :
-		_geometric_size(10.0), _number_of_subdivisions_per_direction(10),
-		_msh(MeshLib::MeshGenerator::generateRegularHexMesh(_geometric_size, _number_of_subdivisions_per_direction)),
-		_project_name("test"), _mshNodesSearcher(*_msh), _ply0(nullptr)
-	{
-		// create geometry
-		std::vector<GeoLib::Point*>* pnts (new std::vector<GeoLib::Point*>);
-		pnts->push_back(new GeoLib::Point(0.0, 0.0, 0.0));
-		pnts->push_back(new GeoLib::Point(_geometric_size, 0.0, 0.0));
-		pnts->push_back(new GeoLib::Point(_geometric_size, _geometric_size, 0.0));
-		pnts->push_back(new GeoLib::Point(0.0, _geometric_size, 0.0));
-		pnts->push_back(new GeoLib::Point(0.0, 0.0, _geometric_size));
-		pnts->push_back(new GeoLib::Point(_geometric_size, 0.0, _geometric_size));
-		pnts->push_back(new GeoLib::Point(_geometric_size, _geometric_size, _geometric_size));
-		pnts->push_back(new GeoLib::Point(0.0, _geometric_size, _geometric_size));
-
-		std::vector<GeoLib::Polyline*>* plys (new std::vector<GeoLib::Polyline*>);
-		_ply0 = new GeoLib::Polyline(*pnts); // vertical polyline
-		_ply0->addPoint(0);
-		_ply0->addPoint(4);
-		plys->push_back(_ply0);
-		GeoLib::Polyline* ply1 = new GeoLib::Polyline(*pnts); // polygon for left surface
-		ply1->addPoint(0);
-		ply1->addPoint(3);
-		ply1->addPoint(7);
-		ply1->addPoint(4);
-		ply1->addPoint(0);
-		plys->push_back(ply1);
-
-		std::vector<GeoLib::Surface*>* sfcs (new std::vector<GeoLib::Surface*>);
-		_sfc1 = GeoLib::Surface::createSurface(*ply1);
-		sfcs->push_back(_sfc1);
-
-		_geo_objs.addPointVec(pnts,_project_name);
-		_geo_objs.addPolylineVec(plys, _project_name);
-		_geo_objs.addSurfaceVec(sfcs, _project_name);
-	}
-
-protected:
-	const double _geometric_size;
-	const std::size_t _number_of_subdivisions_per_direction;
-	std::unique_ptr<MeshLib::Mesh> _msh;
-	GeoLib::GEOObjects _geo_objs;
-	std::string _project_name;
-	MeshGeoToolsLib::MeshNodeSearcher _mshNodesSearcher;
-	GeoLib::Polyline* _ply0;
-	GeoLib::Surface* _sfc1;
-};
-
-template <class T_TASK>
-class NodeIDtoNodeObject
-{
-public:
-	NodeIDtoNodeObject(const MeshLib::Mesh& msh, const T_TASK& task)
-	: _msh(msh), _task(task) {}
-
-	double operator()(std::size_t node_id)
-	{
-		return _task(*_msh.getNode(node_id));
-	}
-
-private:
-	const MeshLib::Mesh& _msh;
-	const T_TASK& _task;
-};
-
-TEST_F(NumLibSpatialFunctionQuad, Linear)
+TEST(NumLib, SpatialFunctionLinear)
 {
 	// f(x,y,z) = 1 + 2x + 3y + 4z
 	std::array<double,4> f_coeff = {{1, 2, 3, 4}};
 	MathLib::LinearFunction<double,3> linear_f(f_coeff);
-	std::vector<double> expected = {{1, 3, 5, 7, 9, 11, 13, 15, 17, 19, 21}};
 
 	NumLib::SpatialFunctionLinear f(linear_f);
-	const std::vector<std::size_t>& vec_node_ids = _mshNodesSearcher.getMeshNodesAlongPolyline(*_ply0).getNodeIDs();
-	std::vector<double> interpolated_values(vec_node_ids.size());
-	NodeIDtoNodeObject<NumLib::SpatialFunctionLinear> task(*_msh, f);
-	std::transform(vec_node_ids.begin(), vec_node_ids.end(), interpolated_values.begin(), task);
-
-	ASSERT_ARRAY_NEAR(expected, interpolated_values, expected.size(), std::numeric_limits<double>::epsilon());
-}
 
-TEST_F(NumLibSpatialFunctionHex, Linear)
-{
-	// f(x,y,z) = 1 + 2x + 3y + 4z
-	std::array<double,4> f_coeff = {{1, 2, 3, 4}};
-	MathLib::LinearFunction<double,3> linear_f(f_coeff);
-	std::vector<double> expected(std::pow(_number_of_subdivisions_per_direction+1, 2));
-	const double dL = _geometric_size / _number_of_subdivisions_per_direction;
-	for (std::size_t i=0; i<expected.size(); i++) {
-		double x = 0;
-		double y = (i%(_number_of_subdivisions_per_direction+1)) * dL;
-		double z = (i/(_number_of_subdivisions_per_direction+1)) * dL;
-		expected[i] = f_coeff[0] + f_coeff[1]*x + f_coeff[2]*y + f_coeff[3]*z;
+	ASSERT_DOUBLE_EQ(1., f(GeoLib::Point(0,0,0)));
+	ASSERT_DOUBLE_EQ(10., f(GeoLib::Point(1,1,1)));
+	ASSERT_DOUBLE_EQ(-8, f(GeoLib::Point(-1,-1,-1)));
+	for (std::size_t k(0); k < 5; ++k)
+	{
+		GeoLib::Point pt(0, 0, 0);
+		for (unsigned i = 0; i < 3; ++i)
+			pt[i] = (double) rand() - (RAND_MAX / 2.0);
+		double expected = 1+2*pt[0]+3*pt[1]+4*pt[2];
+		ASSERT_DOUBLE_EQ(expected, f(pt));
 	}
-
-	NumLib::SpatialFunctionLinear f(linear_f);
-	const std::vector<std::size_t>& vec_node_ids = _mshNodesSearcher.getMeshNodesAlongSurface(*_sfc1).getNodeIDs();
-	std::vector<double> interpolated_values(vec_node_ids.size());
-	NodeIDtoNodeObject<NumLib::SpatialFunctionLinear> task(*_msh, f);
-	std::transform(vec_node_ids.begin(), vec_node_ids.end(), interpolated_values.begin(), task);
-
-	ASSERT_ARRAY_NEAR(expected, interpolated_values, expected.size(), std::numeric_limits<double>::epsilon());
 }
 
-TEST_F(NumLibSpatialFunctionQuad, InterpolationPolyline)
+TEST(NumLib, SpatialFunctionInterpolationPolyline)
 {
+	// create geometry
+	GeoLib::Point pt1(0.0, 0.0, 0.0);
+	GeoLib::Point pt2(10.0, 0.0, 0.0);
+	std::vector<GeoLib::Point*> pnts = {&pt1, &pt2};
+	GeoLib::Polyline ply0(pnts);
+	ply0.addPoint(0);
+	ply0.addPoint(1);
+
+	// define a function
 	const std::vector<std::size_t> vec_point_ids = {{0, 1}};
 	const std::vector<double> vec_point_values = {{0., 100.}};
-	std::vector<double> expected = {{0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100}};
-
-	NumLib::LinearInterpolationAlongPolyline interpolate(
-			*_ply0, vec_point_ids, vec_point_values, std::numeric_limits<double>::epsilon(), std::numeric_limits<double>::max());
-	const std::vector<std::size_t>& vec_node_ids = _mshNodesSearcher.getMeshNodesAlongPolyline(*_ply0).getNodeIDs();
-	std::vector<double> interpolated_values(vec_node_ids.size());
-	NodeIDtoNodeObject<NumLib::LinearInterpolationAlongPolyline> task(*_msh, interpolate);
-	std::transform(vec_node_ids.begin(), vec_node_ids.end(), interpolated_values.begin(), task);
-
-	ASSERT_ARRAY_NEAR(expected, interpolated_values, expected.size(), std::numeric_limits<double>::epsilon());
-}
-
-TEST_F(NumLibSpatialFunctionHex, InterpolationPolyline)
-{
-	const std::vector<std::size_t> vec_point_ids = {{0, 4}};
-	const std::vector<double> vec_point_values = {{0., 100.}};
-	std::vector<double> expected = {{0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100}};
-
 	NumLib::LinearInterpolationAlongPolyline interpolate(
-			*_ply0, vec_point_ids, vec_point_values, std::numeric_limits<double>::epsilon(), std::numeric_limits<double>::max());
-	const std::vector<std::size_t>& vec_node_ids = _mshNodesSearcher.getMeshNodesAlongPolyline(*_ply0).getNodeIDs();
-	std::vector<double> interpolated_values(vec_node_ids.size());
-	NodeIDtoNodeObject<NumLib::LinearInterpolationAlongPolyline> task(*_msh, interpolate);
-	std::transform(vec_node_ids.begin(), vec_node_ids.end(), interpolated_values.begin(), task);
-
-	ASSERT_ARRAY_NEAR(expected, interpolated_values, expected.size(), std::numeric_limits<double>::epsilon());
+			ply0, vec_point_ids, vec_point_values, std::numeric_limits<double>::epsilon(), std::numeric_limits<double>::max());
+
+	// normal
+	for (unsigned k=0; k<10; k++)
+		ASSERT_DOUBLE_EQ(k*10., interpolate(GeoLib::Point(k,0,0)));
+
+	// failure
+	// x
+	ASSERT_DOUBLE_EQ(std::numeric_limits<double>::max(), interpolate(GeoLib::Point(-1,0,0)));
+	ASSERT_DOUBLE_EQ(std::numeric_limits<double>::max(), interpolate(GeoLib::Point(11,0,0)));
+	// y
+	ASSERT_DOUBLE_EQ(std::numeric_limits<double>::max(), interpolate(GeoLib::Point(0,1,0)));
+	ASSERT_DOUBLE_EQ(std::numeric_limits<double>::max(), interpolate(GeoLib::Point(0,-1,0)));
+	// z
+	ASSERT_DOUBLE_EQ(std::numeric_limits<double>::max(), interpolate(GeoLib::Point(0,0,1)));
+	ASSERT_DOUBLE_EQ(std::numeric_limits<double>::max(), interpolate(GeoLib::Point(0,0,-1)));
 }
 
-TEST_F(NumLibSpatialFunctionQuad, InterpolationSurface)
+TEST(NumLib, SpatialFunctionInterpolationSurface)
 {
+	// create geometry
+	GeoLib::Point pt1(0.0, 0.0, 0.0);
+	GeoLib::Point pt2(10.0, 0.0, 0.0);
+	GeoLib::Point pt3(10.0, 10.0, 0.0);
+	GeoLib::Point pt4(0.0, 10.0, 0.0);
+	std::vector<GeoLib::Point*> pnts = {&pt1, &pt2, &pt3, &pt4};
+	GeoLib::Polyline ply0(pnts);
+	ply0.addPoint(0);
+	ply0.addPoint(1);
+	ply0.addPoint(2);
+	ply0.addPoint(3);
+	ply0.addPoint(0);
+	std::unique_ptr<GeoLib::Surface>  sfc1(GeoLib::Surface::createSurface(ply0));
+
+	// define a function
 	const std::vector<std::size_t> vec_point_ids = {{0, 1, 2, 3}};
 	const std::vector<double> vec_point_values = {{0., 100., 100., 0.}};
-	std::vector<double> expected(_msh->getNNodes());
-	for (std::size_t i=0; i<_msh->getNNodes(); i++) {
-		expected[i] = (i%(_number_of_subdivisions_per_direction+1)) * 10;
-	}
-
-	NumLib::LinearInterpolationOnSurface interpolate(*_sfc1, vec_point_ids, vec_point_values, std::numeric_limits<double>::max());
-	const std::vector<std::size_t>& vec_node_ids = _mshNodesSearcher.getMeshNodesAlongSurface(*_sfc1).getNodeIDs();
-	std::vector<double> interpolated_values(vec_node_ids.size());
-	NodeIDtoNodeObject<NumLib::LinearInterpolationOnSurface> task(*_msh, interpolate);
-	std::transform(vec_node_ids.begin(), vec_node_ids.end(), interpolated_values.begin(), task);
-
-	// the machine epsilon for double is too small for this test
-	ASSERT_ARRAY_NEAR(expected, interpolated_values, expected.size(), std::numeric_limits<float>::epsilon());
-}
+	NumLib::LinearInterpolationOnSurface interpolate(*sfc1, vec_point_ids, vec_point_values, std::numeric_limits<double>::max());
 
-TEST_F(NumLibSpatialFunctionHex, InterpolationSurface)
-{
-	const std::vector<std::size_t> vec_point_ids = {{0, 3, 7, 4}};
-	const std::vector<double> vec_point_values = {{0., 100., 100., 0.}};
-	std::vector<double> expected(std::pow(_number_of_subdivisions_per_direction+1, 2));
-	for (std::size_t i=0; i<expected.size(); i++) {
-		expected[i] = (i%(_number_of_subdivisions_per_direction+1)) * 10;
+	// normal
+	for (unsigned k=0; k<10; k++) {
+		ASSERT_DOUBLE_EQ(k*10., interpolate(GeoLib::Point(k,0,0)));
+		ASSERT_DOUBLE_EQ(k*10., interpolate(GeoLib::Point(k,5,0)));
+		ASSERT_DOUBLE_EQ(k*10., interpolate(GeoLib::Point(k,10,0)));
 	}
 
-	NumLib::LinearInterpolationOnSurface interpolate(*_sfc1, vec_point_ids, vec_point_values, std::numeric_limits<double>::max());
-	const std::vector<std::size_t>& vec_node_ids = _mshNodesSearcher.getMeshNodesAlongSurface(*_sfc1).getNodeIDs();
-	std::vector<double> interpolated_values(vec_node_ids.size());
-	NodeIDtoNodeObject<NumLib::LinearInterpolationOnSurface> task(*_msh, interpolate);
-	std::transform(vec_node_ids.begin(), vec_node_ids.end(), interpolated_values.begin(), task);
-
-	ASSERT_ARRAY_NEAR(expected, interpolated_values, expected.size(), std::numeric_limits<float>::epsilon());
+	// failure
+	// x, y
+	ASSERT_DOUBLE_EQ(std::numeric_limits<double>::max(), interpolate(GeoLib::Point(-1,-1,0)));
+	ASSERT_DOUBLE_EQ(std::numeric_limits<double>::max(), interpolate(GeoLib::Point(11,-1,0)));
+	ASSERT_DOUBLE_EQ(std::numeric_limits<double>::max(), interpolate(GeoLib::Point(11,11,0)));
+	ASSERT_DOUBLE_EQ(std::numeric_limits<double>::max(), interpolate(GeoLib::Point(-1,11,0)));
+	// z
+	ASSERT_DOUBLE_EQ(std::numeric_limits<double>::max(), interpolate(GeoLib::Point(0,0,1)));
+	ASSERT_DOUBLE_EQ(std::numeric_limits<double>::max(), interpolate(GeoLib::Point(0,0,-1)));
 }
 
-