diff --git a/GeoLib/Grid.h b/GeoLib/Grid.h
index ed9f2a2e5487c6eadca8073db8afdef61c6f1645..61715d5c1107acd1c1728cb2695215487b080b0b 100644
--- a/GeoLib/Grid.h
+++ b/GeoLib/Grid.h
@@ -131,7 +131,8 @@ public:
 	 * @return a pointer to the point with the smallest distance within the grid cells that are
 	 * outlined above or nullptr
 	 */
-	POINT* getNearestPoint(POINT const& pnt) const
+	template <typename P>
+	POINT* getNearestPoint(P const& pnt) const
 	{
 		std::size_t coords[3];
 		getGridCoords(pnt, coords);
@@ -201,7 +202,7 @@ public:
 			} // end while
 		} // end else
 
-		double len (sqrt(MathLib::sqrDist(pnt.getCoords(), nearest_pnt->getCoords())));
+		double len(sqrt(MathLib::sqrDist(pnt, *nearest_pnt)));
 		// search all other grid cells within the cube with the edge nodes
 		std::vector<std::vector<POINT*> const*> vecs_of_pnts;
 		getPntVecsOfGridCellsIntersectingCube(pnt, len, vecs_of_pnts);
@@ -211,7 +212,7 @@ public:
 			std::vector<POINT*> const& pnts(*(vecs_of_pnts[j]));
 			const std::size_t n_pnts(pnts.size());
 			for (std::size_t k(0); k<n_pnts; k++) {
-				const double sqr_dist (MathLib::sqrDist(pnt.getCoords(), pnts[k]->getCoords()));
+				const double sqr_dist(MathLib::sqrDist(pnt, *pnts[k]));
 				if (sqr_dist < sqr_min_dist) {
 					sqr_min_dist = sqr_dist;
 					nearest_pnt = pnts[k];
@@ -232,7 +233,9 @@ public:
 	 * @param pnts (output) vector of vectors of points within grid cells that intersects
 	 * the axis aligned cube
 	 */
-	void getPntVecsOfGridCellsIntersectingCube(POINT const& center, double half_len, std::vector<std::vector<POINT*> const*>& pnts) const;
+	template <typename P>
+	void getPntVecsOfGridCellsIntersectingCube(P const& center, double half_len,
+		std::vector<std::vector<POINT*> const*>& pnts) const;
 
 	void getPntVecsOfGridCellsIntersectingCuboid(
 		MathLib::Point3d const& min_pnt,
@@ -354,11 +357,13 @@ private:
 	 * ordered in the same sequence as above described
 	 * @param coords coordinates of the grid cell
 	 */
-	void getPointCellBorderDistances(POINT const& pnt,
+	template <typename P>
+	void getPointCellBorderDistances(P const& pnt,
 	                                 double dists[6],
 	                                 std::size_t const* const coords) const;
 
-	bool calcNearestPointInGridCell(POINT const& pnt, std::size_t const* const coords,
+	template <typename P>
+	bool calcNearestPointInGridCell(P const& pnt, std::size_t const* const coords,
 	                                double &sqr_min_dist,
 	                                POINT* &nearest_pnt) const
 	{
@@ -367,7 +372,7 @@ private:
 		if (pnts.empty()) return false;
 
 		const std::size_t n_pnts(pnts.size());
-		sqr_min_dist = MathLib::sqrDist(pnts[0]->getCoords(), pnt.getCoords());
+		sqr_min_dist = MathLib::sqrDist(*pnts[0], pnt);
 		nearest_pnt = pnts[0];
 		for (std::size_t i(1); i < n_pnts; i++) {
 			const double sqr_dist(MathLib::sqrDist(pnts[i]->getCoords(), pnt.getCoords()));
@@ -393,11 +398,13 @@ private:
 };
 
 template<typename POINT>
-void Grid<POINT>::getPntVecsOfGridCellsIntersectingCube(POINT const& center,
+template <typename P>
+void Grid<POINT>::getPntVecsOfGridCellsIntersectingCube(P const& center,
                                                         double half_len,
                                                         std::vector<std::vector<POINT*> const*>& pnts) const
 {
-	double tmp_pnt[3] = { center[0] - half_len, center[1] - half_len, center[2] - half_len }; // min
+	MathLib::Point3d tmp_pnt{
+		{{center[0]-half_len, center[1]-half_len, center[2]-half_len}}}; // min
 	std::size_t min_coords[3];
 	getGridCoords(tmp_pnt, min_coords);
 
@@ -543,7 +550,8 @@ void Grid<POINT>::getGridCoords(T const& pnt, std::size_t* coords) const
 }
 
 template <typename POINT>
-void Grid<POINT>::getPointCellBorderDistances(POINT const& pnt,
+template <typename P>
+void Grid<POINT>::getPointCellBorderDistances(P const& pnt,
                                               double dists[6],
                                               std::size_t const* const coords) const
 {
diff --git a/GeoLib/Point.h b/GeoLib/Point.h
index 941f80b94521a3a7ce07c42eee1ece0dc23c9165..231050fd7cb2e7892571f814396c4923513b9c41 100644
--- a/GeoLib/Point.h
+++ b/GeoLib/Point.h
@@ -43,11 +43,6 @@ public:
 		MathLib::Point3dWithID(), GeoLib::GeoObject()
 	{}
 
-	Point(double const* x) :
-		MathLib::Point3dWithID(std::array<double,3>({{x[0], x[1], x[2]}}), 0),
-		GeoLib::GeoObject()
-	{}
-
 	Point(MathLib::Point3d const& x, std::size_t id) :
 		MathLib::Point3dWithID(x, id), GeoLib::GeoObject()
 	{}
diff --git a/GeoLib/Surface.cpp b/GeoLib/Surface.cpp
index a68de11567f2d89fe03360e22a43ec50e815f848..df6e075b5650445d031cf76219607fc44e6bdcce 100644
--- a/GeoLib/Surface.cpp
+++ b/GeoLib/Surface.cpp
@@ -115,17 +115,17 @@ const Triangle* Surface::operator[] (std::size_t i) const
 	return _sfc_triangles[i];
 }
 
-bool Surface::isPntInBoundingVolume(Point const& pnt) const
+bool Surface::isPntInBoundingVolume(MathLib::Point3d const& pnt) const
 {
 	return _bounding_volume->containsPoint (pnt);
 }
 
-bool Surface::isPntInSfc (Point const& pnt) const
+bool Surface::isPntInSfc(MathLib::Point3d const& pnt) const
 {
 	return (findTriangle(pnt)!=nullptr);
 }
 
-const Triangle* Surface::findTriangle (Point const& pnt) const
+const Triangle* Surface::findTriangle (MathLib::Point3d const& pnt) const
 {
 	for (std::size_t k(0); k<_sfc_triangles.size(); k++) {
 		if (_sfc_triangles[k]->containsPoint (pnt)) {
diff --git a/GeoLib/Surface.h b/GeoLib/Surface.h
index 7468b4294108d083b314de01640c9587e16474b6..75f64b228efdf3bd1999c01f1635d0ff4631bb00 100644
--- a/GeoLib/Surface.h
+++ b/GeoLib/Surface.h
@@ -61,14 +61,14 @@ public:
 	/**
 	 * is the given point in the bounding volume of the surface
 	 */
-	bool isPntInBoundingVolume(Point const& pnt) const;
+	bool isPntInBoundingVolume(MathLib::Point3d const& pnt) const;
 
 	/**
 	 * is the given point pnt located in the surface
 	 * @param pnt the point
 	 * @return true if the point is contained in the surface
 	 */
-	bool isPntInSfc (Point const& pnt) const;
+	bool isPntInSfc(MathLib::Point3d const& pnt) const;
 
 	/**
 	 * find a triangle in which the given point is located
@@ -76,7 +76,7 @@ public:
 	 * @return a pointer to a triangle. nullptr is returned if the point is not
 	 * contained in the surface
 	 */
-	const Triangle* findTriangle (Point const& pnt) const;
+	const Triangle* findTriangle(MathLib::Point3d const& pnt) const;
 
 	const std::vector<Point*> *getPointVec() const { return &_sfc_pnts; }
 
diff --git a/GeoLib/Triangle.cpp b/GeoLib/Triangle.cpp
index 214719674934a1cc22561c53b06d414b51328119..ea7137191526cc136a1b8b5e2b15649c87ff1142 100644
--- a/GeoLib/Triangle.cpp
+++ b/GeoLib/Triangle.cpp
@@ -63,7 +63,7 @@ void Triangle::setTriangle (std::size_t pnt_a, std::size_t pnt_b, std::size_t pn
 	_longest_edge = sqrt (_longest_edge);
 }
 
-bool Triangle::containsPoint(Point const& q, double eps) const
+bool Triangle::containsPoint(MathLib::Point3d const& q, double eps) const
 {
 	GeoLib::Point const& a(*(_pnts[_pnt_ids[0]]));
 	GeoLib::Point const& b(*(_pnts[_pnt_ids[1]]));
diff --git a/GeoLib/Triangle.h b/GeoLib/Triangle.h
index 6668a215d467d32714123e33037f33e6a7a2fc7c..cbc22cc00e01f309e95625774b37004ae6ab9b87 100644
--- a/GeoLib/Triangle.h
+++ b/GeoLib/Triangle.h
@@ -68,7 +68,7 @@ public:
 	 * @param eps Checks the 'epsilon'-neighbourhood
 	 * @return true, if point is in triangle, else false
 	 */
-	bool containsPoint(Point const& q, double eps = std::numeric_limits<float>::epsilon()) const;
+	bool containsPoint(MathLib::Point3d const& q, double eps = std::numeric_limits<float>::epsilon()) const;
 
 	/**
 	 * projects the triangle points to the x-y-plane and
diff --git a/MathLib/MathTools.h b/MathLib/MathTools.h
index f037e46526a1a58fdf73a965a39175d91b46a325..cd9eeb80a28c4ab3f4db2426542a94210911635b 100644
--- a/MathLib/MathTools.h
+++ b/MathLib/MathTools.h
@@ -23,7 +23,7 @@
 #include <omp.h>
 #endif
 
-#include "TemplatePoint.h"
+#include "Point3d.h"
 
 namespace MathLib
 {
@@ -107,13 +107,6 @@ void crossProd (const double u[3], const double v[3], double r[3]);
 double calcProjPntToLineAndDists(const double p[3], const double a[3],
                                  const double b[3], double &lambda, double &d0);
 
-template <typename POINT_T>
-typename POINT_T::FP_T sqrDist(POINT_T const& p0, POINT_T const& p1)
-{
-	typename POINT_T::FP_T const v[3] = {p1[0]-p0[0], p1[1]-p0[1], p1[2]-p0[2]};
-	return MathLib::scalarProduct<typename POINT_T::FP_T,3>(v,v);
-}
-
 template <typename T, std::size_t DIM>
 bool operator==(TemplatePoint<T,DIM> const& a, TemplatePoint<T,DIM> const& b)
 {
@@ -121,6 +114,14 @@ bool operator==(TemplatePoint<T,DIM> const& a, TemplatePoint<T,DIM> const& b)
 	return (sqr_dist < pow(std::numeric_limits<T>::epsilon(),2));
 }
 
+/// Computes the squared dist between the two points p0 and p1.
+inline
+double sqrDist(MathLib::Point3d const& p0, MathLib::Point3d const& p1)
+{
+	const double v[3] = {p1[0] - p0[0], p1[1] - p0[1], p1[2] - p0[2]};
+	return scalarProduct<double,3>(v,v);
+}
+
 /** squared dist between double arrays p0 and p1 (size of arrays is 3) */
 inline
 double sqrDist(const double* p0, const double* p1)
diff --git a/MeshGeoToolsLib/GeoMapper.cpp b/MeshGeoToolsLib/GeoMapper.cpp
index e381b7464aad113f2aa8f08c95199c4464539b79..b4023e16335e5ec9091282712dffaf04c2f984f1 100644
--- a/MeshGeoToolsLib/GeoMapper.cpp
+++ b/MeshGeoToolsLib/GeoMapper.cpp
@@ -153,8 +153,7 @@ float GeoMapper::getDemElevation(GeoLib::Point const& pnt) const
 
 double GeoMapper::getMeshElevation(double x, double y, double min_val, double max_val) const
 {
-	double coords[3] = {x,y,0};
-	const MeshLib::Node* pnt = _grid->getNearestPoint(coords);
+	const MeshLib::Node* pnt = _grid->getNearestPoint(MathLib::Point3d{{{x,y,0}}});
 	const std::vector<MeshLib::Element*> elements (_mesh->getNode(pnt->getID())->getElements());
 	GeoLib::Point* intersection (nullptr);
 
diff --git a/MeshGeoToolsLib/MeshNodeSearcher.cpp b/MeshGeoToolsLib/MeshNodeSearcher.cpp
index 2136909d4e68fbc9f9c01de63c9634c6b298b769..69677009cd3c82410cfed77a6c90df35fadd2a85 100644
--- a/MeshGeoToolsLib/MeshNodeSearcher.cpp
+++ b/MeshGeoToolsLib/MeshNodeSearcher.cpp
@@ -74,7 +74,7 @@ std::vector<std::size_t> MeshNodeSearcher::getMeshNodeIDs(GeoLib::GeoObject cons
 
 boost::optional<std::size_t> MeshNodeSearcher::getMeshNodeIDForPoint(GeoLib::Point const& pnt) const
 {
-	const MeshLib::Node* found = _mesh_grid.getNearestPoint(pnt.getCoords());
+	const MeshLib::Node* found = _mesh_grid.getNearestPoint(pnt);
 	if (found)
 		return found->getID();
 	else
diff --git a/MeshGeoToolsLib/MeshNodesAlongSurface.cpp b/MeshGeoToolsLib/MeshNodesAlongSurface.cpp
index bf25042943627827e00acb23fec9e1b55023e0c1..9fc50ec6b6ed2d440c746475491b6b97886744b0 100644
--- a/MeshGeoToolsLib/MeshNodesAlongSurface.cpp
+++ b/MeshGeoToolsLib/MeshNodesAlongSurface.cpp
@@ -32,9 +32,9 @@ MeshNodesAlongSurface::MeshNodesAlongSurface(
 	// loop over all nodes
 	for (std::size_t i = 0; i < n_nodes; i++) {
 		auto* node = mesh_nodes[i];
-		if (!sfc.isPntInBoundingVolume(node->getCoords()))
+		if (!sfc.isPntInBoundingVolume(*node))
 			continue;
-		if (sfc.isPntInSfc(node->getCoords())) {
+		if (sfc.isPntInSfc(*node)) {
 			_msh_node_ids.push_back(node->getID());
 		}
 	}
diff --git a/MeshLib/MeshEditing/MeshRevision.cpp b/MeshLib/MeshEditing/MeshRevision.cpp
index d20b7eab051c18df39fb23bc8f7d835d50a965b6..be60481c3df572e4ddcf8f10a46efb76742c6b93 100644
--- a/MeshLib/MeshEditing/MeshRevision.cpp
+++ b/MeshLib/MeshEditing/MeshRevision.cpp
@@ -200,7 +200,7 @@ std::vector<std::size_t> MeshRevision::collapseNodeIndices(double eps) const
 		MeshLib::Node const*const node(nodes[k]);
 		if (node->getID() != k)
 			continue;
-		grid.getPntVecsOfGridCellsIntersectingCube(node->getCoords(), half_eps, node_vectors);
+		grid.getPntVecsOfGridCellsIntersectingCube(*node, half_eps, node_vectors);
 
 		const std::size_t nVectors(node_vectors.size());
 		for (std::size_t i = 0; i < nVectors; ++i)
diff --git a/NumLib/Function/LinearInterpolationOnSurface.cpp b/NumLib/Function/LinearInterpolationOnSurface.cpp
index f5fd90ee7cacf80a8227a22a262dc8a1fc6beaf6..6f7ef7702a809f3c0122ead2ed83a813dada29c4 100644
--- a/NumLib/Function/LinearInterpolationOnSurface.cpp
+++ b/NumLib/Function/LinearInterpolationOnSurface.cpp
@@ -41,10 +41,9 @@ LinearInterpolationOnSurface::LinearInterpolationOnSurface(
 
 double LinearInterpolationOnSurface::operator()(const MathLib::Point3d& pnt) const
 {
-	const double* coords = pnt.getCoords();
-	if (!_sfc.isPntInBoundingVolume(coords))
+	if (!_sfc.isPntInBoundingVolume(pnt))
 		return _default_value;
-	auto* tri = _sfc.findTriangle(coords);
+	auto* tri = _sfc.findTriangle(pnt);
 	if (tri == nullptr)
 		return _default_value;
 
@@ -58,11 +57,14 @@ double LinearInterpolationOnSurface::operator()(const MathLib::Point3d& pnt) con
 			pnt_values[j] = _default_value;
 		}
 	}
-	double val = interpolateInTri(*tri, pnt_values.data(), coords);
+	double val = interpolateInTri(*tri, pnt_values.data(), pnt);
 	return val;
 }
 
-double LinearInterpolationOnSurface::interpolateInTri(const GeoLib::Triangle &tri, double const* const vertex_values, double const* const pnt) const
+double LinearInterpolationOnSurface::interpolateInTri(
+	const GeoLib::Triangle &tri,
+	double const* const vertex_values,
+	MathLib::Point3d const& pnt) const
 {
 	std::vector<GeoLib::Point> pnts;
 	for (unsigned i=0; i<3; i++)
diff --git a/NumLib/Function/LinearInterpolationOnSurface.h b/NumLib/Function/LinearInterpolationOnSurface.h
index f40f0e1a3a6506a5134d31b22edaa31fd7e6e023..b77ca7c7bed081000b9e7b4c4ea55a24e568e867 100644
--- a/NumLib/Function/LinearInterpolationOnSurface.h
+++ b/NumLib/Function/LinearInterpolationOnSurface.h
@@ -81,7 +81,9 @@ private:
 	 * @param pnt
 	 * @return
 	 */
-	double interpolateInTri(const GeoLib::Triangle &tri, double const* const vertex_values, double const* const pnt) const;
+	double interpolateInTri(const GeoLib::Triangle &tri,
+		double const* const vertex_values,
+		MathLib::Point3d const& pnt) const;
 
 	/// a surface object
 	const GeoLib::Surface& _sfc;
diff --git a/SimpleTests/MeshTests/MeshSearchTest.cpp b/SimpleTests/MeshTests/MeshSearchTest.cpp
index f7af931ad442af992f214157cb6c4994f4b74143..6514489d17d9b8793403ff23a6addd661da70fef 100644
--- a/SimpleTests/MeshTests/MeshSearchTest.cpp
+++ b/SimpleTests/MeshTests/MeshSearchTest.cpp
@@ -59,7 +59,7 @@ void testMeshGridAlgorithm(MeshLib::Mesh const*const mesh,
 		INFO ("[MeshGridAlgorithm] searching %d points ...", pnts_for_search.size());
 		clock_t start = clock();
 		for (size_t k(0); k<n_pnts_for_search; k++) {
-			MeshLib::Node const* node(mesh_grid.getNearestPoint(pnts_for_search[k]->getCoords()));
+			MeshLib::Node const* node(mesh_grid.getNearestPoint(*pnts_for_search[k]));
 			idx_found_nodes.push_back(node->getID());
 		}
 		clock_t stop = clock();
@@ -83,7 +83,7 @@ void testMeshGridAlgorithm(MeshLib::Mesh const*const mesh,
 		INFO ("[MeshGridAlgorithm] searching %d points ...", pnts_for_search.size());
 		clock_t start = clock();
 		for (size_t k(0); k<n_pnts_for_search; k++) {
-			MeshLib::Node const* node(mesh_grid.getNearestPoint(pnts_for_search[k]->getCoords()));
+			MeshLib::Node const* node(mesh_grid.getNearestPoint(pnts_for_search[k]));
 			idx_found_nodes.push_back(node->getID());
 		}
 		clock_t stop = clock();
@@ -138,7 +138,7 @@ int main(int argc, char *argv[])
 	std::vector<GeoLib::Point*> pnts_for_search;
 	unsigned n(std::min(static_cast<unsigned>(nodes.size()), number_arg.getValue()));
 	for (size_t k(0); k<n; k++) {
-		pnts_for_search.push_back(new GeoLib::Point(nodes[k]->getCoords()));
+		pnts_for_search.push_back(new GeoLib::Point(nodes[k]));
 	}
 
 	std::vector<size_t> idx_found_nodes;