diff --git a/FileIO/MeshIO/GMSHAdaptiveMeshDensity.cpp b/FileIO/MeshIO/GMSHAdaptiveMeshDensity.cpp
index c9a903ea067cee58df86c4b3ac918c2d0b4a4290..5249c2cb6527de0a969a4ca1913e0f166c97aa6f 100644
--- a/FileIO/MeshIO/GMSHAdaptiveMeshDensity.cpp
+++ b/FileIO/MeshIO/GMSHAdaptiveMeshDensity.cpp
@@ -29,6 +29,9 @@
 #include "Polyline.h"
 #endif
 
+// MathLib
+#include "MathTools.h"
+
 namespace FileIO
 {
 GMSHAdaptiveMeshDensity::GMSHAdaptiveMeshDensity(double pnt_density, double station_density,
diff --git a/GeoLib/GeoObject.h b/GeoLib/GeoObject.h
index a7848b59e0891ce69bef4622d6c12c4419526366..88ab3853926b1e66f6c79b38b013c6ea6689a75e 100644
--- a/GeoLib/GeoObject.h
+++ b/GeoLib/GeoObject.h
@@ -2,7 +2,8 @@
  * \file
  * \author Thomas Fischer
  * \date   2010-08-27
- * \brief  Definition of the GeoObject class.
+ * \brief  Base class for classes Point, Polyline, Surface.
+ * \ingroup GeoLib
  *
  * \copyright
  * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
@@ -15,20 +16,14 @@
 #ifndef GEOOBJECT_H_
 #define GEOOBJECT_H_
 
-namespace GeoLib {
-
-/**
- * \ingroup GeoLib
- *
- * \brief Base class for classes Point, Polyline, Surface.
- */
-
-class GeoObject {
+namespace GeoLib
+{
+class GeoObject
+{
 public:
-	GeoObject() {};
-	virtual ~GeoObject() {};
+	GeoObject() {}
+	virtual ~GeoObject() {}
 };
-
 } // end namespace GeoLib
 
 #endif /* GEOOBJECT_H_ */
diff --git a/GeoLib/Grid.h b/GeoLib/Grid.h
index 55cae81f9808b6f8d0e705fff8ce0219cd883ac6..1396fe936923f28177ed696fa15057eb5f3ab5d6 100644
--- a/GeoLib/Grid.h
+++ b/GeoLib/Grid.h
@@ -29,8 +29,11 @@
 #include "StringTools.h"
 #endif
 
-namespace GeoLib {
+// MathLib
+#include "MathTools.h"
 
+namespace GeoLib
+{
 template <typename POINT>
 class Grid : public GeoLib::AABB<POINT>
 {
@@ -181,7 +184,7 @@ public:
 	 * @return a pointer to the point with the smallest distance within the grid cells that are
 	 * outlined above or NULL
 	 */
-	POINT* getNearestPoint(double const*const pnt) const
+	POINT* getNearestPoint(double const* const pnt) const
 	{
 		std::size_t coords[3];
 		getGridCoords(POINT(pnt), coords);
@@ -317,11 +320,13 @@ private:
 	 * ordered in the same sequence as above described
 	 * @param coords coordinates of the grid cell
 	 */
-	void getPointCellBorderDistances(double const*const pnt, double dists[6], std::size_t const* const coords) const;
+	void getPointCellBorderDistances(double const* const pnt,
+	                                 double dists[6],
+	                                 std::size_t const* const coords) const;
 
 	bool calcNearestPointInGridCell(double const* const pnt, std::size_t const* const coords,
-					double &sqr_min_dist,
-					POINT* &nearest_pnt) const
+	                                double &sqr_min_dist,
+	                                POINT* &nearest_pnt) const
 	{
 		const std::size_t grid_idx (coords[0] + coords[1] * _n_steps[0] + coords[2] * _n_steps[0] * _n_steps[1]);
 		std::vector<typename std::add_pointer<typename std::remove_pointer<POINT>::type>::type> const& pnts(_grid_cell_nodes_map[grid_idx]);
@@ -353,10 +358,10 @@ private:
 	std::vector<POINT*>* _grid_cell_nodes_map;
 };
 
-
 template<typename POINT>
 void Grid<POINT>::getPntVecsOfGridCellsIntersectingCube(POINT const& center,
-				double half_len, std::vector<std::vector<POINT*> const*>& pnts) const
+                                                        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
 	std::size_t min_coords[3];
@@ -382,7 +387,8 @@ void Grid<POINT>::getPntVecsOfGridCellsIntersectingCube(POINT const& center,
 
 template<typename POINT>
 void Grid<POINT>::getPntVecsOfGridCellsIntersectingCuboid(POINT const& min_pnt,
-				POINT const& max_pnt, std::vector<std::vector<POINT*> const*>& pnts) const
+                                                          POINT const& max_pnt,
+                                                          std::vector<std::vector<POINT*> const*>& pnts) const
 {
 	std::size_t min_coords[3];
 	getGridCoords(min_pnt, min_coords);
@@ -411,9 +417,9 @@ void Grid<POINT>::createGridGeometry(GeoLib::GEOObjects* geo_obj) const
 	GeoLib::Point const& llf (this->getMinPoint());
 	GeoLib::Point const& urb (this->getMaxPoint());
 
-	const double dx ((urb[0]-llf[0])/_n_steps[0]);
-	const double dy ((urb[1]-llf[1])/_n_steps[1]);
-	const double dz ((urb[2]-llf[2])/_n_steps[2]);
+	const double dx ((urb[0] - llf[0]) / _n_steps[0]);
+	const double dy ((urb[1] - llf[1]) / _n_steps[1]);
+	const double dz ((urb[2] - llf[2]) / _n_steps[2]);
 
 	// create grid names and grid boxes as geometry
 	for (std::size_t i(0); i<_n_steps[0]; i++) {
@@ -422,7 +428,7 @@ void Grid<POINT>::createGridGeometry(GeoLib::GEOObjects* geo_obj) const
 
 				std::string name("Grid-");
 				name += BaseLib::number2str<std::size_t>(i);
-				name +="-";
+				name += "-";
 				name += BaseLib::number2str<std::size_t>(j);
 				name += "-";
 				name += BaseLib::number2str<std::size_t> (k);
@@ -439,18 +445,17 @@ void Grid<POINT>::createGridGeometry(GeoLib::GEOObjects* geo_obj) const
 				points->push_back(new GeoLib::Point(llf[0]+(i+1)*dx, llf[1]+j*dy, llf[2]+(k+1)*dz));
 				geo_obj->addPointVec(points, grid_names[grid_names.size()-1], NULL);
 
-				std::vector<GeoLib::Polyline*>* plys (new std::vector<GeoLib::Polyline*>);
+				std::vector<GeoLib::Polyline*>* plys (
+				        new std::vector<GeoLib::Polyline*>);
 				GeoLib::Polyline* ply0 (new GeoLib::Polyline(*points));
-				for (std::size_t l(0); l<4; l++) {
+				for (std::size_t l(0); l < 4; l++)
 					ply0->addPoint(l);
-				}
 				ply0->addPoint(0);
 				plys->push_back(ply0);
 
 				GeoLib::Polyline* ply1 (new GeoLib::Polyline(*points));
-				for (std::size_t l(4); l<8; l++) {
+				for (std::size_t l(4); l < 8; l++)
 					ply1->addPoint(l);
-				}
 				ply1->addPoint(4);
 				plys->push_back(ply1);
 
@@ -474,7 +479,8 @@ void Grid<POINT>::createGridGeometry(GeoLib::GEOObjects* geo_obj) const
 				ply5->addPoint(7);
 				plys->push_back(ply5);
 
-				geo_obj->addPolylineVec(plys, grid_names[grid_names.size()-1], NULL);
+				geo_obj->addPolylineVec(plys,
+				                        grid_names[grid_names.size() - 1], NULL);
 			}
 		}
 	}
@@ -501,19 +507,19 @@ void Grid<POINT>::getGridCoords(POINT const& pnt, std::size_t* coords) const
 }
 
 template <typename POINT>
-void Grid<POINT>::getPointCellBorderDistances(double const*const pnt, double dists[6], std::size_t const* const coords) const
+void Grid<POINT>::getPointCellBorderDistances(double const* const pnt,
+                                              double dists[6],
+                                              std::size_t const* const coords) const
 {
-
-	dists[0] = (pnt[2] - this->_min_pnt[2] + coords[2]*_step_sizes[2]); // bottom
+	dists[0] = (pnt[2] - this->_min_pnt[2] + coords[2] * _step_sizes[2]); // bottom
 	dists[5] = (_step_sizes[2] - dists[0]); // top
 
-	dists[1] = (pnt[1] - this->_min_pnt[1] + coords[1]*_step_sizes[1]); // front
+	dists[1] = (pnt[1] - this->_min_pnt[1] + coords[1] * _step_sizes[1]); // front
 	dists[3] = (_step_sizes[1] - dists[1]); // back
 
-	dists[4] = (pnt[0] - this->_min_pnt[0] + coords[0]*_step_sizes[0]); // left
+	dists[4] = (pnt[0] - this->_min_pnt[0] + coords[0] * _step_sizes[0]); // left
 	dists[2] = (_step_sizes[0] - dists[4]); // right
 }
-
 } // end namespace GeoLib
 
 #endif /* MESHGRID_H_ */
diff --git a/GeoLib/Polygon.cpp b/GeoLib/Polygon.cpp
index fb1733cfd694aae5faee132217bda02799682c40..1bd58aecda5ed0cfae65ac829a1a5ae7e12a9fda 100644
--- a/GeoLib/Polygon.cpp
+++ b/GeoLib/Polygon.cpp
@@ -71,12 +71,12 @@ bool Polygon::isPntInPolygon (GeoLib::Point const & pnt) const
 	    max_aabb_pnt[1] < pnt[1])
 		return false;
 
-	size_t n_intersections (0);
+	std::size_t n_intersections (0);
 	GeoLib::Point s;
 
 	if (_simple_polygon_list.empty ()) {
-		const size_t n_nodes (getNumberOfPoints() - 1);
-		for (size_t k(0); k < n_nodes; k++) {
+		const std::size_t n_nodes (getNumberOfPoints() - 1);
+		for (std::size_t k(0); k < n_nodes; k++) {
 			if (((*(getPoint(k)))[1] <= pnt[1] && pnt[1] <= (*(getPoint(k + 1)))[1]) ||
 			    ((*(getPoint(k + 1)))[1] <= pnt[1] && pnt[1] <= (*(getPoint(k)))[1])) {
 				switch (getEdgeType(k, pnt))
@@ -115,9 +115,9 @@ bool Polygon::isPntInPolygon(double x, double y, double z) const
 
 bool Polygon::isPolylineInPolygon(const Polyline& ply) const
 {
-	size_t ply_size (ply.getNumberOfPoints()), cnt (0);
-	for (size_t k(0); k < ply_size; k++) {
-		if (isPntInPolygon (*(ply[k]))) {
+	std::size_t ply_size (ply.getNumberOfPoints()), cnt (0);
+	for (std::size_t k(0); k < ply_size; k++) {
+		if (isPntInPolygon (*(ply.getPoint(k)))) {
 			cnt++;
 		}
 	}
@@ -129,18 +129,18 @@ bool Polygon::isPolylineInPolygon(const Polyline& ply) const
 
 bool Polygon::isPartOfPolylineInPolygon(const Polyline& ply) const
 {
-	const size_t ply_size (ply.getNumberOfPoints());
+	const std::size_t ply_size (ply.getNumberOfPoints());
 	// check points
-	for (size_t k(0); k < ply_size; k++) {
-		if (isPntInPolygon (*(ply[k]))) {
+	for (std::size_t k(0); k < ply_size; k++) {
+		if (isPntInPolygon (*(ply.getPoint(k)))) {
 			return true;
 		}
 	}
 	// check segment intersections
 	GeoLib::Point* s (new GeoLib::Point (0,0,0));
-	const size_t n_nodes(getNumberOfPoints() - 1);
-	for (size_t k(0); k < ply_size - 1; k++) {
-		for (size_t j(0); j < n_nodes; j++) {
+	const std::size_t n_nodes(getNumberOfPoints() - 1);
+	for (std::size_t k(0); k < ply_size - 1; k++) {
+		for (std::size_t j(0); j < n_nodes; j++) {
 			if (MathLib::lineSegmentIntersect(*(getPoint(j)), *(getPoint(j + 1)),
 							*(ply.getPoint(k)), *(ply.getPoint(k + 1)), *s)) {
 				delete s;
@@ -155,11 +155,11 @@ bool Polygon::isPartOfPolylineInPolygon(const Polyline& ply) const
 
 bool Polygon::getNextIntersectionPointPolygonLine (GeoLib::Point const & a,
                 GeoLib::Point const & b, GeoLib::Point* intersection_pnt,
-                size_t& seg_num) const
+                std::size_t& seg_num) const
 {
 	if (_simple_polygon_list.empty()) {
-		const size_t n_nodes(getNumberOfPoints() - 1);
-		for (size_t k(seg_num); k < n_nodes; k++) {
+		const std::size_t n_nodes(getNumberOfPoints() - 1);
+		for (std::size_t k(seg_num); k < n_nodes; k++) {
 			if (MathLib::lineSegmentIntersect(*(getPoint(k)), *(getPoint(k + 1)), a, b, *intersection_pnt)) {
 				seg_num = k;
 				return true;
@@ -169,8 +169,8 @@ bool Polygon::getNextIntersectionPointPolygonLine (GeoLib::Point const & a,
 		for (std::list<Polygon*>::const_iterator it(_simple_polygon_list.begin()); it
 					!= _simple_polygon_list.end(); ++it) {
 			const Polygon* polygon(*it);
-			const size_t n_nodes_simple_polygon(polygon->getNumberOfPoints() - 1);
-			for (size_t k(0); k < n_nodes_simple_polygon; k++) {
+			const std::size_t n_nodes_simple_polygon(polygon->getNumberOfPoints() - 1);
+			for (std::size_t k(0); k < n_nodes_simple_polygon; k++) {
 				if (MathLib::lineSegmentIntersect(*(polygon->getPoint(k)), *(polygon->getPoint(k + 1)),
 								a, b, *intersection_pnt)) {
 					seg_num = k;
@@ -203,7 +203,7 @@ void Polygon::computeListOfSimplePolygons ()
 		(*it)->initialise ();
 }
 
-EdgeType::value Polygon::getEdgeType (size_t k, GeoLib::Point const & pnt) const
+EdgeType::value Polygon::getEdgeType (std::size_t k, GeoLib::Point const & pnt) const
 {
 	switch (getLocationOfPoint(k, pnt))
 	{
@@ -238,9 +238,9 @@ void Polygon::ensureCWOrientation ()
 {
 	// *** pre processing: rotate points to xy-plan
 	// *** copy points to vector - last point is identical to the first
-	size_t n_pnts (this->getNumberOfPoints() - 1);
+	std::size_t n_pnts (this->getNumberOfPoints() - 1);
 	std::vector<GeoLib::Point*> tmp_polygon_pnts;
-	for (size_t k(0); k < n_pnts; k++)
+	for (std::size_t k(0); k < n_pnts; k++)
 		tmp_polygon_pnts.push_back (new GeoLib::Point (*(this->getPoint(k))));
 
 	// *** calculate supporting plane (plane normal and
@@ -254,12 +254,12 @@ void Polygon::ensureCWOrientation ()
 		// rotate copied points into x-y-plane
 		MathLib::rotatePointsToXY(plane_normal, tmp_polygon_pnts);
 
-	for (size_t k(0); k < tmp_polygon_pnts.size(); k++)
+	for (std::size_t k(0); k < tmp_polygon_pnts.size(); k++)
 		(*(tmp_polygon_pnts[k]))[2] = 0.0; // should be -= d but there are numerical errors
 
 	// *** get the left most upper point
-	size_t min_x_max_y_idx (0); // for orientation check
-	for (size_t k(0); k < n_pnts; k++)
+	std::size_t min_x_max_y_idx (0); // for orientation check
+	for (std::size_t k(0); k < n_pnts; k++)
 		if ((*(tmp_polygon_pnts[k]))[0] <= (*(tmp_polygon_pnts[min_x_max_y_idx]))[0])
 		{
 			if ((*(tmp_polygon_pnts[k]))[0] < (*(tmp_polygon_pnts[min_x_max_y_idx]))[0])
@@ -293,19 +293,19 @@ void Polygon::ensureCWOrientation ()
 	if (orient == MathLib::CCW)
 	{
 		// switch orientation
-		size_t tmp_n_pnts (n_pnts);
+		std::size_t tmp_n_pnts (n_pnts);
 		tmp_n_pnts++; // include last point of polygon (which is identical to the first)
-		for (size_t k(0); k < tmp_n_pnts / 2; k++)
+		for (std::size_t k(0); k < tmp_n_pnts / 2; k++)
 			std::swap (_ply_pnt_ids[k], _ply_pnt_ids[tmp_n_pnts - 1 - k]);
 	}
 
-	for (size_t k(0); k < n_pnts; k++)
+	for (std::size_t k(0); k < n_pnts; k++)
 		delete tmp_polygon_pnts[k];
 }
 
 void Polygon::splitPolygonAtIntersection (std::list<Polygon*>::iterator polygon_it)
 {
-	size_t idx0 (0), idx1 (0);
+	std::size_t idx0 (0), idx1 (0);
 	while (polygon_it != _simple_polygon_list.end())
 	{
 		GeoLib::Point* intersection_pnt (new GeoLib::Point);
@@ -316,7 +316,7 @@ void Polygon::splitPolygonAtIntersection (std::list<Polygon*>::iterator polygon_
 		if (!is_simple)
 		{
 			// adding intersection point to pnt_vec
-			size_t intersection_pnt_id (_ply_pnts.size());
+			std::size_t intersection_pnt_id (_ply_pnts.size());
 			const_cast<std::vector<Point*>& >(_ply_pnts).push_back (intersection_pnt);
 
 			// split Polygon
@@ -325,10 +325,10 @@ void Polygon::splitPolygonAtIntersection (std::list<Polygon*>::iterator polygon_
 
 			GeoLib::Polygon* polygon0 (new GeoLib::Polygon(
 			                                   (*polygon_it)->getPointsVec(), false));
-			for (size_t k(0); k <= idx0; k++)
+			for (std::size_t k(0); k <= idx0; k++)
 				polygon0->addPoint ((*polygon_it)->getPointID (k));
 			polygon0->addPoint (intersection_pnt_id);
-			for (size_t k(idx1 + 1); k < (*polygon_it)->getNumberOfPoints(); k++)
+			for (std::size_t k(idx1 + 1); k < (*polygon_it)->getNumberOfPoints(); k++)
 				polygon0->addPoint ((*polygon_it)->getPointID (k));
 			if (!polygon0->initialise())
 			{
@@ -339,7 +339,7 @@ void Polygon::splitPolygonAtIntersection (std::list<Polygon*>::iterator polygon_
 			GeoLib::Polygon* polygon1 (new GeoLib::Polygon(
 			                                   (*polygon_it)->getPointsVec(), false));
 			polygon1->addPoint (intersection_pnt_id);
-			for (size_t k(idx0 + 1); k <= idx1; k++)
+			for (std::size_t k(idx0 + 1); k <= idx1; k++)
 				polygon1->addPoint ((*polygon_it)->getPointID (k));
 			polygon1->addPoint (intersection_pnt_id);
 			if (!polygon1->initialise())
@@ -365,9 +365,9 @@ void Polygon::splitPolygonAtIntersection (std::list<Polygon*>::iterator polygon_
 
 void Polygon::splitPolygonAtPoint (std::list<GeoLib::Polygon*>::iterator polygon_it)
 {
-	size_t n ((*polygon_it)->getNumberOfPoints() - 1), idx0 (0), idx1(0);
-	size_t* id_vec (new size_t[n]), *perm (new size_t[n]);
-	for (size_t k(0); k < n; k++)
+	std::size_t n ((*polygon_it)->getNumberOfPoints() - 1), idx0 (0), idx1(0);
+	std::size_t* id_vec (new std::size_t[n]), *perm (new std::size_t[n]);
+	for (std::size_t k(0); k < n; k++)
 	{
 		id_vec[k] = (*polygon_it)->getPointID (k);
 		perm[k] = k;
@@ -375,7 +375,7 @@ void Polygon::splitPolygonAtPoint (std::list<GeoLib::Polygon*>::iterator polygon
 
 	BaseLib::quicksort (id_vec, 0, n, perm);
 
-	for (size_t k(0); k < n - 1; k++)
+	for (std::size_t k(0); k < n - 1; k++)
 		if (id_vec[k] == id_vec[k + 1])
 		{
 			idx0 = perm[k];
@@ -388,14 +388,14 @@ void Polygon::splitPolygonAtPoint (std::list<GeoLib::Polygon*>::iterator polygon
 
 			// create two closed polylines
 			GeoLib::Polygon* polygon0 (new GeoLib::Polygon(*(*polygon_it)));
-			for (size_t k(0); k <= idx0; k++)
+			for (std::size_t k(0); k <= idx0; k++)
 				polygon0->addPoint ((*polygon_it)->getPointID (k));
-			for (size_t k(idx1 + 1); k < (*polygon_it)->getNumberOfPoints(); k++)
+			for (std::size_t k(idx1 + 1); k < (*polygon_it)->getNumberOfPoints(); k++)
 				polygon0->addPoint ((*polygon_it)->getPointID (k));
 			polygon0->initialise();
 
 			GeoLib::Polygon* polygon1 (new GeoLib::Polygon(*(*polygon_it)));
-			for (size_t k(idx0); k <= idx1; k++)
+			for (std::size_t k(idx0); k <= idx1; k++)
 				polygon1->addPoint ((*polygon_it)->getPointID (k));
 			polygon1->initialise();
 
@@ -416,12 +416,12 @@ void Polygon::splitPolygonAtPoint (std::list<GeoLib::Polygon*>::iterator polygon
 }
 
 GeoLib::Polygon* createPolygonFromCircle (GeoLib::Point const& middle_pnt, double radius,
-                                          std::vector<GeoLib::Point*> & pnts, size_t resolution)
+                                          std::vector<GeoLib::Point*> & pnts, std::size_t resolution)
 {
-	const size_t off_set (pnts.size());
+	const std::size_t off_set (pnts.size());
 	// create points
 	double angle (2.0 * M_PI / resolution);
-	for (size_t k(0); k < resolution; k++)
+	for (std::size_t k(0); k < resolution; k++)
 	{
 		GeoLib::Point* pnt (new GeoLib::Point(middle_pnt.getCoords()));
 		(*pnt)[0] += radius * cos (k * angle);
@@ -431,7 +431,7 @@ GeoLib::Polygon* createPolygonFromCircle (GeoLib::Point const& middle_pnt, doubl
 
 	// create polygon
 	GeoLib::Polygon* polygon (new GeoLib::Polygon (pnts, false));
-	for (size_t k(0); k < resolution; k++)
+	for (std::size_t k(0); k < resolution; k++)
 		polygon->addPoint (k + off_set);
 	polygon->addPoint (off_set);
 
@@ -443,12 +443,12 @@ bool operator==(Polygon const& lhs, Polygon const& rhs)
 	if (lhs.getNumberOfPoints() != rhs.getNumberOfPoints())
 		return false;
 
-	const size_t n(lhs.getNumberOfPoints());
-	const size_t start_pnt(lhs.getPointID(0));
+	const std::size_t n(lhs.getNumberOfPoints());
+	const std::size_t start_pnt(lhs.getPointID(0));
 
 	// search start point of first polygon in second polygon
 	bool nfound(true);
-	size_t k(0);
+	std::size_t k(0);
 	for (; k < n-1 && nfound; k++) {
 		if (start_pnt == rhs.getPointID(k)) {
 			nfound = false;
@@ -472,7 +472,7 @@ bool operator==(Polygon const& lhs, Polygon const& rhs)
 
 	// same direction - start point of first polygon at arbitrary position in second polygon
 	if (lhs.getPointID(1) == rhs.getPointID(k+1)) {
-		size_t j(k+2);
+		std::size_t j(k+2);
 		for (; j<n-1; j++) {
 			if (lhs.getPointID(j-k) != rhs.getPointID(j)) {
 				return false;
@@ -491,7 +491,7 @@ bool operator==(Polygon const& lhs, Polygon const& rhs)
 		WARN("operator==(Polygon const& lhs, Polygon const& rhs) - not tested case (implementation is probably buggy) - please contact thomas.fischer@ufz.de mentioning the problem.");
 		// in second polygon
 		if (lhs.getPointID(1) == rhs.getPointID(k-1)) {
-			size_t j(k-2);
+			std::size_t j(k-2);
 			for (; j>0; j--) {
 				if (lhs.getPointID(k-2-j) != rhs.getPointID(j)) {
 					return false;
diff --git a/GeoLib/Polyline.cpp b/GeoLib/Polyline.cpp
index fce2bae427b29b64e9ae9db6a92bd8249afdd119..2c7d732b11a25b853da6b7832631e74fe5071e86 100644
--- a/GeoLib/Polyline.cpp
+++ b/GeoLib/Polyline.cpp
@@ -12,6 +12,9 @@
  *
  */
 
+// STL
+#include <algorithm>
+
 // ThirdParty/logog
 #include "logog/include/logog.hpp"
 
@@ -32,33 +35,34 @@ Polyline::Polyline(const std::vector<Point*>& pnt_vec) :
 Polyline::Polyline(const Polyline& ply) :
 	GeoObject(), _ply_pnts (ply._ply_pnts)
 {
-	for (size_t k(0); k < ply.getNumberOfPoints(); ++k)
+	for (std::size_t k(0); k < ply.getNumberOfPoints(); ++k)
 		_ply_pnt_ids.push_back (ply.getPointID (k));
 
 	if (ply.getNumberOfPoints() > 0)
-		for (size_t k(0); k < ply.getNumberOfPoints(); ++k)
+		for (std::size_t k(0); k < ply.getNumberOfPoints(); ++k)
 			_length.push_back (ply.getLength (k));
-
 }
 
 void Polyline::write(std::ostream &os) const
 {
-	size_t size(_ply_pnt_ids.size());
-	for (size_t k(0); k < size; k++)
+	std::size_t size(_ply_pnt_ids.size());
+	for (std::size_t k(0); k < size; k++)
 		os << *(_ply_pnts[_ply_pnt_ids[k]]) << std::endl;
 }
 
-void Polyline::addPoint(size_t pnt_id)
+void Polyline::addPoint(std::size_t pnt_id)
 {
 	assert(pnt_id < _ply_pnts.size());
-	size_t n_pnts (_ply_pnt_ids.size());
+	std::size_t n_pnts (_ply_pnt_ids.size());
 
-	// don't insert point if ID if this would result in identical IDs for two adjacent points
-	if (n_pnts>0 && _ply_pnt_ids[n_pnts-1] == pnt_id) return;
+	// don't insert point if this would result in identical IDs for two adjacent points
+	if (n_pnts > 0 && _ply_pnt_ids[n_pnts - 1] == pnt_id)
+		return;
 
 	_ply_pnt_ids.push_back(pnt_id);
 
-	if (n_pnts > 0) {
+	if (n_pnts > 0)
+	{
 		double act_dist (sqrt(MathLib::sqrDist (_ply_pnts[_ply_pnt_ids[n_pnts - 1]],
 		                                        _ply_pnts[pnt_id])));
 		double dist_until_now (0.0);
@@ -69,49 +73,66 @@ void Polyline::addPoint(size_t pnt_id)
 	}
 }
 
-void Polyline::insertPoint(size_t pos, size_t pnt_id)
+void Polyline::insertPoint(std::size_t pos, std::size_t pnt_id)
 {
 	assert(pnt_id < _ply_pnts.size());
-	assert(pos < _ply_pnt_ids.size());
+	assert(pos <= _ply_pnt_ids.size());
+
+	if (pos == _ply_pnt_ids.size()) {
+		addPoint(pnt_id);
+		return;
+	}
 
 	// check if inserting pnt_id would result in two identical IDs for adjacent points
-	if (pos == 0 && pnt_id == _ply_pnt_ids[0]) return;
-	else if (pos == (_ply_pnt_ids.size()-1) && pnt_id == _ply_pnt_ids[pos]) return;
-	else if (pnt_id == _ply_pnt_ids[pos-1] || pnt_id == _ply_pnt_ids[pos]) return;
+	if (pos == 0 && pnt_id == _ply_pnt_ids[0])
+		return;
+	else if (pos == (_ply_pnt_ids.size() - 1) && pnt_id == _ply_pnt_ids[pos])
+		return;
+	else if (pnt_id == _ply_pnt_ids[pos - 1] || pnt_id == _ply_pnt_ids[pos])
+		return;
 
-	std::vector<size_t>::iterator it(_ply_pnt_ids.begin() + pos);
+	std::vector<std::size_t>::iterator it(_ply_pnt_ids.begin() + pos);
 	_ply_pnt_ids.insert(it, pnt_id);
 
 	if (_ply_pnt_ids.size() > 1) {
 		// update the _length vector
 		if (pos == 0) {
 			// insert at first position
-			double act_dist(sqrt(MathLib::sqrDist(_ply_pnts[_ply_pnt_ids[1]], _ply_pnts[pnt_id])));
-			_length.insert(_length.begin(), act_dist);
-			const size_t s(_length.size());
-			for (size_t k(1); k<s; k++) {
-				_length[k] += _length[0];
-			}
+			double act_dist(sqrt(MathLib::sqrDist(_ply_pnts[_ply_pnt_ids[1]],
+			                                      _ply_pnts[pnt_id])));
+			_length.insert(_length.begin() + 1, act_dist);
+			const std::size_t s(_length.size());
+			for (std::size_t k(2); k < s; k++)
+				_length[k] += _length[1];
 		} else {
-			if (pos == _ply_pnt_ids.size()-1) {
+			if (pos == _ply_pnt_ids.size() - 1) {
 				// insert at last position
-				double act_dist(sqrt(MathLib::sqrDist(_ply_pnts[_ply_pnt_ids[_ply_pnt_ids.size()-2]], _ply_pnts[pnt_id])));
+				double act_dist(sqrt(MathLib::sqrDist(
+				                             _ply_pnts[_ply_pnt_ids[_ply_pnt_ids.size() - 2]],
+				                             _ply_pnts[pnt_id])));
 				double dist_until_now (0.0);
 				if (_ply_pnt_ids.size() > 2)
 					dist_until_now = _length[_ply_pnt_ids.size() - 2];
 
-				_length.insert(_length.begin()+pos, dist_until_now + act_dist);
+				_length.insert(_length.begin() + pos, dist_until_now + act_dist);
 			} else {
 				// insert at arbitrary position within the vector
 				double dist_until_now (0.0);
-				if (pos > 1) {
-					dist_until_now = _length[pos-2];
-				}
-				double len_seg0(sqrt(MathLib::sqrDist(_ply_pnts[_ply_pnt_ids[pos-1]], _ply_pnts[pnt_id])));
-				double len_seg1(sqrt(MathLib::sqrDist(_ply_pnts[_ply_pnt_ids[pos+1]], _ply_pnts[pnt_id])));
-				_length[pos-1] = dist_until_now + len_seg0;
-				std::vector<double>::iterator it(_length.begin()+pos);
-				_length.insert(it, _length[pos-1]+len_seg1);
+				if (pos > 1)
+					dist_until_now = _length[pos - 1];
+				double len_seg0(sqrt(MathLib::sqrDist(
+				                             _ply_pnts[_ply_pnt_ids[pos - 1]],
+				                             _ply_pnts[pnt_id])));
+				double len_seg1(sqrt(MathLib::sqrDist(
+				                             _ply_pnts[_ply_pnt_ids[pos + 1]],
+				                             _ply_pnts[pnt_id])));
+				double update_dist(
+				        len_seg0 + len_seg1 - (_length[pos] - dist_until_now));
+				_length[pos] = dist_until_now + len_seg0;
+				std::vector<double>::iterator it(_length.begin() + pos + 1);
+				_length.insert(it, _length[pos] + len_seg1);
+				for (it = _length.begin() + pos + 2; it != _length.end(); ++it)
+					*it += update_dist;
 			}
 		}
 	}
@@ -122,79 +143,67 @@ void Polyline::removePoint(std::size_t pos)
 	if (pos >= _ply_pnt_ids.size())
 		return;
 
-	_ply_pnt_ids.erase(_ply_pnt_ids.begin()+pos);
+	_ply_pnt_ids.erase(_ply_pnt_ids.begin() + pos);
 
-	if (pos == _ply_pnt_ids.size()-1) {
-		_length.erase(_length.begin()+pos);
+	if (pos == _ply_pnt_ids.size())
+	{
+		_length.erase(_length.begin() + pos);
 		return;
 	}
 
-	const size_t n_ply_pnt_ids(_ply_pnt_ids.size());
+	const std::size_t n_ply_pnt_ids(_ply_pnt_ids.size());
 	if (pos == 0) {
 		double seg_length(_length[0]);
-		for (unsigned k(0); k<n_ply_pnt_ids; k++) {
-			_length[k] = _length[k+1] - seg_length;
-		}
+		for (std::size_t k(0); k < n_ply_pnt_ids; k++)
+			_length[k] = _length[k + 1] - seg_length;
 		_length.pop_back();
 	} else {
-		const double len_seg0(_length[pos] - _length[pos-1]);
-		const double len_seg1(_length[pos+1] - _length[pos]);
-		_length.erase(_length.begin()+pos);
-		const double len_new_seg(sqrt(MathLib::sqrDist(_ply_pnts[_ply_pnt_ids[pos-1]], _ply_pnts[pos])));
+		const double len_seg0(_length[pos] - _length[pos - 1]);
+		const double len_seg1(_length[pos + 1] - _length[pos]);
+		_length.erase(_length.begin() + pos);
+		const double len_new_seg(sqrt(MathLib::sqrDist(_ply_pnts[_ply_pnt_ids[pos - 1]],
+		                                               _ply_pnts[_ply_pnt_ids[pos]])));
 		double seg_length_diff(len_new_seg - len_seg0 - len_seg1);
 
-		for (unsigned k(pos); k<n_ply_pnt_ids; k++) {
+		for (std::size_t k(pos); k < n_ply_pnt_ids; k++)
 			_length[k] += seg_length_diff;
-		}
 	}
 }
 
-size_t Polyline::getNumberOfPoints() const
+std::size_t Polyline::getNumberOfPoints() const
 {
 	return _ply_pnt_ids.size();
 }
 
 bool Polyline::isClosed() const
 {
+	if (_ply_pnt_ids.size() < 3)
+		return false;
+
 	if (_ply_pnt_ids.front() == _ply_pnt_ids.back())
 		return true;
 	else
 		return false;
 }
 
-bool Polyline::isPointIDInPolyline(size_t pnt_id) const
+bool Polyline::isPointIDInPolyline(std::size_t pnt_id) const
 {
-	const size_t n_ply_pnt_ids(_ply_pnt_ids.size());
-	size_t k(0);
-	while (k<n_ply_pnt_ids && _ply_pnt_ids[k] != pnt_id) {
-		k++;
-	}
-
-	if (k == n_ply_pnt_ids) {
-		return false;
-	}
-	return true;
+	return std::find(_ply_pnt_ids.begin(), _ply_pnt_ids.end(), pnt_id) != _ply_pnt_ids.end();
 }
 
-size_t Polyline::getPointID(size_t i) const
+std::size_t Polyline::getPointID(std::size_t i) const
 {
 	assert(i < _ply_pnt_ids.size());
 	return _ply_pnt_ids[i];
 }
 
-void Polyline::setPointID(size_t idx, size_t id)
+void Polyline::setPointID(std::size_t idx, std::size_t id)
 {
 	assert(idx < _ply_pnt_ids.size());
 	_ply_pnt_ids[idx] = id;
 }
 
-const Point* Polyline::operator[](size_t i) const
-{
-	assert(i < _ply_pnt_ids.size());
-	return _ply_pnts[_ply_pnt_ids[i]];
-}
-
-const Point* Polyline::getPoint(size_t i) const
+const Point* Polyline::getPoint(std::size_t i) const
 {
 	assert(i < _ply_pnt_ids.size());
 	return _ply_pnts[_ply_pnt_ids[i]];
@@ -205,27 +214,22 @@ std::vector<Point*> const& Polyline::getPointsVec () const
 	return _ply_pnts;
 }
 
-double Polyline::getLength (size_t k) const
+double Polyline::getLength (std::size_t k) const
 {
 	assert(k < _length.size());
 	return _length[k];
 }
 
-const std::vector<double>& Polyline::getLengthVec () const
-{
-	return _length;
-}
-
 Polyline* Polyline::constructPolylineFromSegments(const std::vector<Polyline*> &ply_vec,
                                                   double prox)
 {
-	size_t nLines = ply_vec.size();
+	std::size_t nLines = ply_vec.size();
 
 	Polyline* new_ply = new Polyline(*ply_vec[0]);
 	std::vector<GeoLib::Point*> pnt_vec(new_ply->getPointsVec());
 
 	std::vector<Polyline*> local_ply_vec;
-	for (size_t i = 1; i < nLines; i++)
+	for (std::size_t i = 1; i < nLines; i++)
 		local_ply_vec.push_back(ply_vec[i]);
 
 	while (!local_ply_vec.empty())
@@ -237,18 +241,18 @@ Polyline* Polyline::constructPolylineFromSegments(const std::vector<Polyline*> &
 		{
 			if (pnt_vec == (*it)->getPointsVec())
 			{
-				size_t nPoints((*it)->getNumberOfPoints());
+				std::size_t nPoints((*it)->getNumberOfPoints());
 
 				//if (new_ply->getPointID(0) == (*it)->getPointID(0))
 				if (pointsAreIdentical(pnt_vec, new_ply->getPointID(0),
 				                       (*it)->getPointID(0), prox))
 				{
 					Polyline* tmp = new Polyline((*it)->getPointsVec());
-					for (size_t k = 0; k < nPoints; k++)
+					for (std::size_t k = 0; k < nPoints; k++)
 						tmp->addPoint((*it)->getPointID(nPoints - k - 1));
 
-					size_t new_ply_size(new_ply->getNumberOfPoints());
-					for (size_t k = 1; k < new_ply_size; k++)
+					std::size_t new_ply_size(new_ply->getNumberOfPoints());
+					for (std::size_t k = 1; k < new_ply_size; k++)
 						tmp->addPoint(new_ply->getPointID(k));
 					delete new_ply;
 					new_ply = tmp;
@@ -259,8 +263,8 @@ Polyline* Polyline::constructPolylineFromSegments(const std::vector<Polyline*> &
 				                            (*it)->getPointID(nPoints - 1), prox))
 				{
 					Polyline* tmp = new Polyline(**it);
-					size_t new_ply_size(new_ply->getNumberOfPoints());
-					for (size_t k = 1; k < new_ply_size; k++)
+					std::size_t new_ply_size(new_ply->getNumberOfPoints());
+					for (std::size_t k = 1; k < new_ply_size; k++)
 						tmp->addPoint(new_ply->getPointID(k));
 					delete new_ply;
 					new_ply = tmp;
@@ -273,7 +277,7 @@ Polyline* Polyline::constructPolylineFromSegments(const std::vector<Polyline*> &
 				                                                - 1),
 				                            (*it)->getPointID(0), prox))
 				{
-					for (size_t k = 1; k < nPoints; k++)
+					for (std::size_t k = 1; k < nPoints; k++)
 						new_ply->addPoint((*it)->getPointID(k));
 					ply_found = true;
 				}
@@ -284,7 +288,7 @@ Polyline* Polyline::constructPolylineFromSegments(const std::vector<Polyline*> &
 				                                                - 1),
 				                            (*it)->getPointID(nPoints - 1), prox))
 				{
-					for (size_t k = 1; k < nPoints; k++)
+					for (std::size_t k = 1; k < nPoints; k++)
 						new_ply->addPoint((*it)->getPointID(nPoints - k - 1));
 					ply_found = true;
 				}
@@ -308,31 +312,17 @@ Polyline* Polyline::constructPolylineFromSegments(const std::vector<Polyline*> &
 	return new_ply;
 }
 
-bool Polyline::pointsAreIdentical(const std::vector<Point*> &pnt_vec,
-                                  size_t i,
-                                  size_t j,
-                                  double prox)
-{
-	if (i == j)
-		return true;
-	return MathLib::checkDistance( *pnt_vec[i], *pnt_vec[j], prox );
-}
-
-Polyline* Polyline::closePolyline(const Polyline& ply)
+void Polyline::closePolyline()
 {
-	if (ply.getNumberOfPoints() > 2)
-	{
-		Polyline* new_ply = new Polyline(ply);
-		if (ply.isClosed())
-			return new_ply;
-		new_ply->addPoint(new_ply->getPointID(0));
-		return new_ply;
+	if (getNumberOfPoints() < 2) {
+		ERR("Polyline::closePolyline(): Input polyline needs to be composed of at least three points.");
+	}
+	if (! isClosed()) {
+		addPoint(getPointID(0));
 	}
-	ERR("Error in Polyline::closePolyline() - Input polyline needs to be composed of at least three points.");
-	return NULL;
 }
 
-Location::type Polyline::getLocationOfPoint (size_t k, GeoLib::Point const & pnt) const
+Location::type Polyline::getLocationOfPoint (std::size_t k, GeoLib::Point const & pnt) const
 {
 	assert (k < _ply_pnt_ids.size() - 1);
 
@@ -349,11 +339,14 @@ Location::type Polyline::getLocationOfPoint (size_t k, GeoLib::Point const & pnt
 		return Location::RIGHT;
 	if (a[0] * b[0] < 0.0 || a[1] * b[1] < 0.0)
 		return Location::BEHIND;
-	if (a[0]*a[0]+a[1]*a[1] < b[0]*b[0]+b[1]*b[1])
+	if (a[0] * a[0] + a[1] * a[1] < b[0] * b[0] + b[1] * b[1])
 		return Location::BEYOND;
-	if (MathLib::sqrDist (&pnt, _ply_pnts[_ply_pnt_ids[k]]) < sqrt(std::numeric_limits<double>::min()))
+	if (MathLib::sqrDist (&pnt,
+	                      _ply_pnts[_ply_pnt_ids[k]]) < sqrt(std::numeric_limits<double>::min()))
 		return Location::SOURCE;
-	if (MathLib::sqrDist (&pnt, _ply_pnts[_ply_pnt_ids[k + 1]]) < sqrt(std::numeric_limits<double>::min()))
+	if (MathLib::sqrDist (&pnt,
+	                      _ply_pnts[_ply_pnt_ids[k + 1]]) <
+	    sqrt(std::numeric_limits<double>::min()))
 		return Location::DESTINATION;
 	return Location::BETWEEN;
 }
@@ -364,7 +357,7 @@ std::ostream& operator<< (std::ostream &os, const Polyline &pl)
 	return os;
 }
 
-bool containsEdge (const Polyline& ply, size_t id0, size_t id1)
+bool containsEdge (const Polyline& ply, std::size_t id0, std::size_t id1)
 {
 	if (id0 == id1)
 	{
@@ -373,11 +366,11 @@ bool containsEdge (const Polyline& ply, size_t id0, size_t id1)
 	}
 	if (id0 > id1)
 		std::swap (id0,id1);
-	const size_t n (ply.getNumberOfPoints() - 1);
-	for (size_t k(0); k < n; k++)
+	const std::size_t n (ply.getNumberOfPoints() - 1);
+	for (std::size_t k(0); k < n; k++)
 	{
-		size_t ply_pnt0 (ply.getPointID (k));
-		size_t ply_pnt1 (ply.getPointID (k + 1));
+		std::size_t ply_pnt0 (ply.getPointID (k));
+		std::size_t ply_pnt1 (ply.getPointID (k + 1));
 		if (ply_pnt0 > ply_pnt1)
 			std::swap (ply_pnt0, ply_pnt1);
 		if (ply_pnt0 == id0 && ply_pnt1 == id1)
@@ -386,14 +379,17 @@ bool containsEdge (const Polyline& ply, size_t id0, size_t id1)
 	return false;
 }
 
-bool isLineSegmentIntersecting (const Polyline& ply, GeoLib::Point const& s0, GeoLib::Point const& s1)
+bool isLineSegmentIntersecting (const Polyline& ply,
+                                GeoLib::Point const& s0,
+                                GeoLib::Point const& s1)
 {
-	const size_t n (ply.getNumberOfPoints() - 1);
+	const std::size_t n (ply.getNumberOfPoints() - 1);
 	bool intersect(false);
 	GeoLib::Point intersection_pnt;
-	for (size_t k(0); k < n && !intersect; k++) {
-		intersect = MathLib::lineSegmentIntersect (*(ply.getPoint(k)), *(ply.getPoint(k+1)), s0, s1, intersection_pnt);
-	}
+	for (std::size_t k(0); k < n && !intersect; k++)
+		intersect = MathLib::lineSegmentIntersect (*(ply.getPoint(k)), *(ply.getPoint(
+		                                                                         k + 1)),
+		                                           s0, s1, intersection_pnt);
 	return intersect;
 }
 
@@ -402,13 +398,21 @@ bool operator==(Polyline const& lhs, Polyline const& rhs)
 	if (lhs.getNumberOfPoints() != rhs.getNumberOfPoints())
 		return false;
 
-	const size_t n(lhs.getNumberOfPoints());
-	for (size_t k(0); k<n; k++) {
+	const std::size_t n(lhs.getNumberOfPoints());
+	for (std::size_t k(0); k < n; k++)
 		if (lhs.getPointID(k) != rhs.getPointID(k))
 			return false;
-	}
 
 	return true;
 }
 
+bool pointsAreIdentical(const std::vector<Point*> &pnt_vec,
+                        std::size_t i,
+                        std::size_t j,
+                        double prox)
+{
+	if (i == j)
+		return true;
+	return MathLib::checkDistance( *pnt_vec[i], *pnt_vec[j], prox );
+}
 } // end namespace GeoLib
diff --git a/GeoLib/Polyline.h b/GeoLib/Polyline.h
index 8a68c07d990cf5092f840df1d0a18e9b17dc7d21..07a65e0a3584565cad96b081f68de65845985a79 100644
--- a/GeoLib/Polyline.h
+++ b/GeoLib/Polyline.h
@@ -15,16 +15,13 @@
 #ifndef POLYLINE_H_
 #define POLYLINE_H_
 
+#include <cmath>
+#include <vector>
+
 // GeoLib
 #include "GeoObject.h"
 #include "Point.h"
 
-// MathLib
-#include "MathTools.h"
-
-#include <cmath>
-#include <vector>
-
 namespace GeoLib
 {
 class Location
@@ -93,10 +90,8 @@ public:
 
 	/**
 	 * Closes a polyline by adding a line segment that connects start- and end-point.
-	 * \param ply A Polyline containing at least three points.
-	 * \return A polygon.
 	 */
-	static Polyline* closePolyline(const Polyline& ply);
+	void closePolyline();
 
 	/// Constructs one polyline from a vector of connected polylines.
 	/// All polylines in this vector need to reference the same point vector.
@@ -133,10 +128,6 @@ public:
 	 */
 	void setPointID(std::size_t idx, std::size_t id);
 
-	/** \brief const access operator for the access to the i-th point of the polyline.
-	 */
-	const Point* operator[](std::size_t i) const;
-
 	/**
 	 * \brief returns the i-th point contained in the polyline
 	 * */
@@ -151,12 +142,6 @@ public:
 	 */
 	double getLength (std::size_t k) const;
 
-	/**
-	 * get the complete length vector
-	 * @return the length vector of the polyline
-	 */
-	const std::vector<double>& getLengthVec () const;
-
 	friend bool operator==(Polyline const& lhs, Polyline const& rhs);
 protected:
 	/**
@@ -170,11 +155,6 @@ protected:
 	 */
 	Location::type getLocationOfPoint (std::size_t k, GeoLib::Point const & pnt) const;
 
-	static bool pointsAreIdentical(const std::vector<Point*> &pnt_vec,
-	                               std::size_t i,
-	                               std::size_t j,
-	                               double prox);
-
 	/** a reference to the vector of pointers to the geometric points */
 	const std::vector<Point*> &_ply_pnts;
 	/** position of pointers to the geometric points */
@@ -190,7 +170,9 @@ std::ostream& operator<< (std::ostream &os, Polyline const& pl);
 
 bool containsEdge (const Polyline& ply, std::size_t id0, std::size_t id1);
 
-bool isLineSegmentIntersecting (const Polyline& ply, GeoLib::Point const& s0, GeoLib::Point const& s1);
+bool isLineSegmentIntersecting (const Polyline& ply,
+                                GeoLib::Point const& s0,
+                                GeoLib::Point const& s1);
 
 /**
  * comparison operator
@@ -200,6 +182,8 @@ bool isLineSegmentIntersecting (const Polyline& ply, GeoLib::Point const& s0, Ge
  */
 bool operator==(Polyline const& lhs, Polyline const& rhs);
 
+bool pointsAreIdentical(const std::vector<Point*> &pnt_vec, std::size_t i, std::size_t j,
+                        double prox);
 } // end namespace
 
 #endif /* POLYLINE_H_ */
diff --git a/Gui/DataView/GEOModels.cpp b/Gui/DataView/GEOModels.cpp
index 3c0fa7667572b008c31f5157c69ffda794f49aa4..ad27aa4a9fd99674a0e9fb78730c5177d9f213b2 100644
--- a/Gui/DataView/GEOModels.cpp
+++ b/Gui/DataView/GEOModels.cpp
@@ -15,6 +15,9 @@
 // ** INCLUDES **
 #include "GEOModels.h"
 
+// ThirdParty/logog
+#include "logog/include/logog.hpp"
+
 #include "GeoTreeModel.h"
 #include "StationTreeModel.h"
 
@@ -65,7 +68,7 @@ void GEOModels::updateGeometry(const std::string &geo_name)
 		}
 	}
 	else
-		std::cout << "Error in GEOModels::updateGeometry() - Geometry \"" << geo_name << "\" not found." << std::endl;
+		ERR("GEOModels::updateGeometry() - Geometry \"%s\" not found.", geo_name.c_str());
 }
 
 void GEOModels::removeGeometry(std::string geo_name, GeoLib::GEOTYPE type)
@@ -80,7 +83,7 @@ void GEOModels::removeGeometry(std::string geo_name, GeoLib::GEOTYPE type)
 
 void GEOModels::addPointVec( std::vector<GeoLib::Point*>* points,
                              std::string &name,
-                             std::map<std::string, size_t>* name_pnt_id_map,
+                             std::map<std::string, std::size_t>* name_pnt_id_map,
                              double eps)
 {
 	GEOObjects::addPointVec(points, name, name_pnt_id_map, eps);
@@ -89,7 +92,7 @@ void GEOModels::addPointVec( std::vector<GeoLib::Point*>* points,
 }
 
 bool GEOModels::appendPointVec(const std::vector<GeoLib::Point*> &points,
-                               const std::string &name, std::vector<size_t>* ids)
+                               const std::string &name, std::vector<std::size_t>* ids)
 {
 	bool ret (GeoLib::GEOObjects::appendPointVec (points, name, ids));
 	// TODO import new points into geo-treeview
@@ -136,7 +139,7 @@ bool GEOModels::removeStationVec( const std::string &name )
 
 void GEOModels::addPolylineVec( std::vector<GeoLib::Polyline*>* lines,
                                 const std::string &name,
-                                std::map<std::string,size_t>* ply_names )
+                                std::map<std::string,std::size_t>* ply_names )
 {
 	GEOObjects::addPolylineVec(lines, name, ply_names);
 	if (lines->empty())
@@ -164,7 +167,7 @@ bool GEOModels::removePolylineVec( const std::string &name )
 
 void GEOModels::addSurfaceVec( std::vector<GeoLib::Surface*>* surfaces,
                                const std::string &name,
-                               std::map<std::string,size_t>* sfc_names )
+                               std::map<std::string,std::size_t>* sfc_names )
 {
 	GEOObjects::addSurfaceVec(surfaces, name, sfc_names);
 	if (surfaces->empty())
@@ -184,7 +187,7 @@ bool GEOModels::appendSurfaceVec(const std::vector<GeoLib::Surface*> &surfaces,
 	else
 	{
 		std::vector<GeoLib::Surface*>* sfc = new std::vector<GeoLib::Surface*>;
-		for (size_t i = 0; i < surfaces.size(); i++)
+		for (std::size_t i = 0; i < surfaces.size(); i++)
 			sfc->push_back(surfaces[i]);
 		this->addSurfaceVec(sfc, name);
 	}
@@ -200,7 +203,7 @@ bool GEOModels::removeSurfaceVec( const std::string &name )
 }
 
 void GEOModels::connectPolylineSegments(const std::string &geoName,
-                                        std::vector<size_t> indexlist,
+                                        std::vector<std::size_t> indexlist,
                                         double proximity,
                                         std::string ply_name,
                                         bool closePly,
@@ -212,7 +215,7 @@ void GEOModels::connectPolylineSegments(const std::string &geoName,
 	{
 		const std::vector<GeoLib::Polyline*>* polylines = plyVec->getVector();
 		std::vector<GeoLib::Polyline*> ply_list;
-		for (size_t i = 0; i < indexlist.size(); i++)
+		for (std::size_t i = 0; i < indexlist.size(); i++)
 			ply_list.push_back( (*polylines)[indexlist[i]] );
 
 		// connect polylines
@@ -227,7 +230,7 @@ void GEOModels::connectPolylineSegments(const std::string &geoName,
 
 			if (closePly)
 			{
-				new_line = GeoLib::Polyline::closePolyline(*new_line);
+				new_line->closePolyline();
 
 				if (triangulatePly)
 				{
@@ -249,8 +252,10 @@ void GEOModels::connectPolylineSegments(const std::string &geoName,
 		OGSError::box("Corresponding geometry not found.");
 }
 
-
-void GEOModels::addNameForElement(const std::string &geometry_name, const GeoLib::GEOTYPE object_type, size_t id, std::string new_name)
+void GEOModels::addNameForElement(const std::string &geometry_name,
+                                  const GeoLib::GEOTYPE object_type,
+                                  std::size_t id,
+                                  std::string new_name)
 {
 	if (object_type == GeoLib::POINT)
 		this->getPointVecObj(geometry_name)->setNameForElement(id, new_name);
@@ -259,33 +264,47 @@ void GEOModels::addNameForElement(const std::string &geometry_name, const GeoLib
 	else if (object_type == GeoLib::SURFACE)
 		this->getSurfaceVecObj(geometry_name)->setNameForElement(id, new_name);
 	else
-		std::cout << "Error in GEOModels::addNameForElement() - Unknown GEOTYPE..." << std::endl;
+		ERR("GEOModels::addNameForElement() - Unknown GEOTYPE %s.",
+		    GeoLib::convertGeoTypeToString(object_type).c_str());
 }
 
-void GEOModels::addNameForObjectPoints(const std::string &geometry_name, const GeoLib::GEOTYPE object_type, const std::string &geo_object_name, const std::string &new_name)
+void GEOModels::addNameForObjectPoints(const std::string &geometry_name,
+                                       const GeoLib::GEOTYPE object_type,
+                                       const std::string &geo_object_name,
+                                       const std::string &new_name)
 {
-	const GeoLib::GeoObject* obj = this->getGEOObject(geometry_name, object_type, geo_object_name);
+	const GeoLib::GeoObject* obj = this->getGEOObject(geometry_name,
+	                                                  object_type,
+	                                                  geo_object_name);
 	GeoLib::PointVec* pnt_vec = this->getPointVecObj(geometry_name);
 	if (object_type == GeoLib::POLYLINE)
 	{
 		const GeoLib::Polyline* ply = dynamic_cast<const GeoLib::Polyline*>(obj);
-		size_t nPoints = ply->getNumberOfPoints();
-		for (size_t i=0; i<nPoints; i++)
-			pnt_vec->setNameForElement(ply->getPointID(i), new_name + "_Point" + BaseLib::number2str(ply->getPointID(i)));
+		std::size_t nPoints = ply->getNumberOfPoints();
+		for (std::size_t i = 0; i < nPoints; i++)
+			pnt_vec->setNameForElement(ply->getPointID(
+			                                   i), new_name + "_Point" +
+			                           BaseLib::number2str(ply->getPointID(i)));
 	}
 	else if (object_type == GeoLib::SURFACE)
 	{
 		const GeoLib::Surface* sfc = dynamic_cast<const GeoLib::Surface*>(obj);
-		size_t nTriangles = sfc->getNTriangles();
-		for (size_t i=0; i<nTriangles; i++)
+		std::size_t nTriangles = sfc->getNTriangles();
+		for (std::size_t i = 0; i < nTriangles; i++)
 		{
 			const GeoLib::Triangle* tri = (*sfc)[i];
-			pnt_vec->setNameForElement((*tri)[0], new_name + "_Point" + BaseLib::number2str((*tri)[0]));
-			pnt_vec->setNameForElement((*tri)[1], new_name + "_Point" + BaseLib::number2str((*tri)[1]));
-			pnt_vec->setNameForElement((*tri)[2], new_name + "_Point" + BaseLib::number2str((*tri)[2]));
+			pnt_vec->setNameForElement((*tri)[0],
+			                           new_name + "_Point" +
+			                           BaseLib::number2str((*tri)[0]));
+			pnt_vec->setNameForElement((*tri)[1],
+			                           new_name + "_Point" +
+			                           BaseLib::number2str((*tri)[1]));
+			pnt_vec->setNameForElement((*tri)[2],
+			                           new_name + "_Point" +
+			                           BaseLib::number2str((*tri)[2]));
 		}
 	}
 	else
-		std::cout << "Error in GEOModels::addNameForElement() - Unknown GEOTYPE..." << std::endl;
+		ERR("GEOModels::addNameForObjectPoints() - Unknown GEOTYPE %s.",
+		    GeoLib::convertGeoTypeToString(object_type).c_str());
 }
-
diff --git a/Gui/DataView/GeoTreeModel.cpp b/Gui/DataView/GeoTreeModel.cpp
index 2b2f7e95c5dd24c3e3fe5dd0926fb7581ba5d35c..d1ed26b8a62b78b8ff861bfc25cdaee5500e0f7e 100644
--- a/Gui/DataView/GeoTreeModel.cpp
+++ b/Gui/DataView/GeoTreeModel.cpp
@@ -57,11 +57,13 @@ void GeoTreeModel::addPointList(QString geoName, const GeoLib::PointVec* pointVe
 		pointVec->getNameOfElementByID(j, pnt_name);
 		QList<QVariant> pnt_data;
 		pnt_data << static_cast<unsigned>(j)
-			     << QString::number(pnt[0], 'f')
-			     << QString::number(pnt[1], 'f')
-			     << QString::number(pnt[2], 'f')
-			     << QString::fromStdString(pnt_name);
-		GeoTreeItem* point(new GeoTreeItem(pnt_data, pointList, static_cast<const GeoLib::Point*>(&pnt)));
+		         << QString::number(pnt[0], 'f')
+		         << QString::number(pnt[1], 'f')
+		         << QString::number(pnt[2], 'f')
+		         << QString::fromStdString(pnt_name);
+		GeoTreeItem* point(new GeoTreeItem(pnt_data,
+		                                   pointList,
+		                                   static_cast<const GeoLib::Point*>(&pnt)));
 		pointList->appendChild(point);
 	}
 
@@ -143,12 +145,12 @@ void GeoTreeModel::addChildren(GeoObjectListItem* plyList,
 		int nPoints = static_cast<int>(lines[i]->getNumberOfPoints());
 		for (int j = 0; j < nPoints; j++)
 		{
-			const GeoLib::Point pnt(*(line[j]));
+			const GeoLib::Point pnt(*(line.getPoint(j)));
 			QList<QVariant> pnt_data;
 			pnt_data << static_cast<int>(line.getPointID(j))
-				     << QString::number(pnt[0], 'f')
-					 << QString::number(pnt[1], 'f')
-					 << QString::number(pnt[2], 'f');
+			         << QString::number(pnt[0], 'f')
+			         << QString::number(pnt[1], 'f')
+			         << QString::number(pnt[2], 'f');
 			TreeItem* child(new TreeItem(pnt_data, lineItem));
 			lineItem->appendChild(child);
 		}
@@ -168,8 +170,9 @@ void GeoTreeModel::addSurfaceList(QString geoName, const GeoLib::SurfaceVec* sur
 
 	if (geo == NULL)
 	{
-		std::cout << "GeoTreeModel::addSurfaceList() - Error: No corresponding geometry found..."
-				  << std::endl;
+		std::cout <<
+		"GeoTreeModel::addSurfaceList() - Error: No corresponding geometry found..."
+		          << std::endl;
 		return;
 	}
 
@@ -193,7 +196,8 @@ void GeoTreeModel::appendSurfaces(const std::string &name, GeoLib::SurfaceVec* s
 			int nChildren = _lists[i]->childCount();
 			for (int j = 0; j < nChildren; j++)
 			{
-				GeoObjectListItem* parent = static_cast<GeoObjectListItem*>(_lists[i]->child(j));
+				GeoObjectListItem* parent =
+				        static_cast<GeoObjectListItem*>(_lists[i]->child(j));
 				if (GeoLib::SURFACE == parent->getType())
 				{
 					this->addChildren(parent, surfaceVec,
@@ -222,7 +226,8 @@ void GeoTreeModel::addChildren(GeoObjectListItem* sfcList,
 		QList<QVariant> surface;
 		std::string sfc_name("");
 		surface_vec->getNameOfElementByID(i, sfc_name);
-		surface << "Surface " + QString::number(i) << QString::fromStdString(sfc_name) << "" << "";
+		surface << "Surface " + QString::number(i) << QString::fromStdString(sfc_name) <<
+		"" << "";
 
 		const GeoLib::Surface &sfc(*(*surfaces)[i]);
 		GeoTreeItem* surfaceItem(new GeoTreeItem(surface, sfcList, &sfc));
@@ -234,8 +239,8 @@ void GeoTreeModel::addChildren(GeoObjectListItem* sfcList,
 			QList<QVariant> elem;
 			const GeoLib::Triangle &triangle(*sfc[j]);
 			elem << j << static_cast<int>(triangle[0])
-				      << static_cast<int>(triangle[1])
-					  << static_cast<int>(triangle[2]);
+			     << static_cast<int>(triangle[1])
+			     << static_cast<int>(triangle[2]);
 			TreeItem* child(new TreeItem(elem, surfaceItem));
 			surfaceItem->appendChild(child);
 
@@ -244,9 +249,9 @@ void GeoTreeModel::addChildren(GeoObjectListItem* sfcList,
 				QList<QVariant> node;
 				const GeoLib::Point &pnt(*(nodesVec[triangle[k]]));
 				node << static_cast<int>(triangle[k])
-					 << QString::number(pnt[0], 'f')
-					 << QString::number(pnt[1], 'f')
-					 << QString::number(pnt[2], 'f');
+				     << QString::number(pnt[0], 'f')
+				     << QString::number(pnt[1], 'f')
+				     << QString::number(pnt[2], 'f');
 				TreeItem* nchild(new TreeItem(node, child));
 				child->appendChild(nchild);
 			}
@@ -297,38 +302,41 @@ vtkPolyDataAlgorithm* GeoTreeModel::vtkSource(const std::string &name, GeoLib::G
 	return NULL;
 }
 
-void GeoTreeModel::setNameForItem(const std::string &name, GeoLib::GEOTYPE type, size_t id, std::string item_name)
+void GeoTreeModel::setNameForItem(const std::string &name,
+                                  GeoLib::GEOTYPE type,
+                                  size_t id,
+                                  std::string item_name)
 {
 	int type_idx(0);
 	int col_idx(1);
 
 	switch(type)
 	{
-		case GeoLib::POINT:
-			type_idx = 0;
-			col_idx = 4; // for points the name is at a different position
-			break;
-		case GeoLib::POLYLINE:
-			type_idx = 1;
-			break;
-		case GeoLib::SURFACE:
-			type_idx = 2;
-			break;
-		case GeoLib::VOLUME:
-			type_idx = 3;
-			break;
-		default:
-			type_idx = -1;
+	case GeoLib::POINT:
+		type_idx = 0;
+		col_idx = 4; // for points the name is at a different position
+		break;
+	case GeoLib::POLYLINE:
+		type_idx = 1;
+		break;
+	case GeoLib::SURFACE:
+		type_idx = 2;
+		break;
+	case GeoLib::VOLUME:
+		type_idx = 3;
+		break;
+	default:
+		type_idx = -1;
 	}
 
-	for (size_t i=0; i<_lists.size(); i++)
+	for (size_t i = 0; i < _lists.size(); i++)
 	{
 		if ( name.compare( _lists[i]->data(0).toString().toStdString() ) == 0 )
 		{
 			TreeItem* object_list = _lists[i]->child(type_idx);
 //			for (int j=0; j<object_list->childCount(); j++)
 //			{
-				TreeItem* item = object_list->child(/*j*/id);
+			TreeItem* item = object_list->child(/*j*/ id);
 //				if (static_cast<size_t>(item->data(0).toInt()) == id)
 //				{
 					item->setData(col_idx, QString::fromStdString(item_name));
diff --git a/Gui/DataView/LinearEditDialog.cpp b/Gui/DataView/LinearEditDialog.cpp
index 5c5ae992a318f6965249589f755f4916f2d98baa..2fb07ed21dffa88e2a4a33b77decf9705d4223dc 100644
--- a/Gui/DataView/LinearEditDialog.cpp
+++ b/Gui/DataView/LinearEditDialog.cpp
@@ -14,20 +14,24 @@
 
 #include "LinearEditDialog.h"
 
-LinearEditDialog::LinearEditDialog(const GeoLib::Polyline &line, const std::vector<size_t> &dis_nodes, const std::vector<double> &dis_values, QDialog* parent)
+LinearEditDialog::LinearEditDialog(const GeoLib::Polyline &line,
+                                   const std::vector<size_t> &dis_nodes,
+                                   const std::vector<double> &dis_values,
+                                   QDialog* parent)
 	: QDialog(parent), _line(line)
 {
 	setupUi(this);
 	setupDialog(dis_nodes, dis_values);
 }
 
-void LinearEditDialog::setupDialog(const std::vector<size_t> &dis_nodes, const std::vector<double> &dis_values)
+void LinearEditDialog::setupDialog(const std::vector<size_t> &dis_nodes,
+                                   const std::vector<double> &dis_values)
 {
 	size_t nPoints(_line.getNumberOfPoints());
 	this->tableWidget->setRowCount(nPoints);
 	QList<QString> indexlist;
 
-	for (size_t i=0; i<nPoints; i++)
+	for (size_t i = 0; i < nPoints; i++)
 	{
 		indexlist.push_back(QString::number(i));
 		QTableWidgetItem *newItem = new QTableWidgetItem("");
@@ -37,7 +41,7 @@ void LinearEditDialog::setupDialog(const std::vector<size_t> &dis_nodes, const s
 	tableWidget->setVerticalHeaderLabels(vHeaders);
 
 	size_t nValues (dis_values.size());
-	for (size_t i=0; i<nValues; i++)
+	for (size_t i = 0; i < nValues; i++)
 		tableWidget->item(dis_nodes[i],0)->setText(QString::number(dis_values[i]));
 }
 
@@ -47,11 +51,11 @@ LinearEditDialog::~LinearEditDialog()
 
 void LinearEditDialog::on_comboBox_currentIndexChanged(int index)
 {
-	if (index>0) //elevation
+	if (index > 0) //elevation
 	{
 		size_t nRows = tableWidget->rowCount();
-		for (size_t i=0; i<nRows; i++)
-			tableWidget->item(i,0)->setText(QString::number(_line[i]->getCoords()[2]));
+		for (size_t i = 0; i < nRows; i++)
+			tableWidget->item(i,0)->setText(QString::number(_line.getPoint(i)->getCoords()[2]));
 	}
 }
 
@@ -60,10 +64,10 @@ void LinearEditDialog::accept()
 	std::vector< std::pair<size_t,double> > linear_values;
 
 	size_t nRows = tableWidget->rowCount();
-	for (size_t i=0; i<nRows; i++)
+	for (size_t i = 0; i < nRows; i++)
 	{
 		QString row_text (tableWidget->item(i,0)->text());
-		if (row_text.length()>0)
+		if (row_text.length() > 0)
 			linear_values.push_back( std::pair<size_t, double>(i, row_text.toDouble()) );
 	}
 
diff --git a/Gui/VtkVis/VtkPolylinesSource.cpp b/Gui/VtkVis/VtkPolylinesSource.cpp
index 9b03b974f802cb9553a071b57f40bbccfe35ba47..615c0296f75ca558148983f075acbe9fc6a39754 100644
--- a/Gui/VtkVis/VtkPolylinesSource.cpp
+++ b/Gui/VtkVis/VtkPolylinesSource.cpp
@@ -27,9 +27,9 @@
 #include <vtkPointData.h>
 #include <vtkPoints.h>
 #include <vtkPolyData.h>
+#include <vtkProperty.h>
 #include <vtkSmartPointer.h>
 #include <vtkStreamingDemandDrivenPipeline.h>
-#include <vtkProperty.h>
 
 vtkStandardNewMacro(VtkPolylinesSource);
 vtkCxxRevisionMacro(VtkPolylinesSource, "$Revision$");
@@ -62,7 +62,7 @@ void VtkPolylinesSource::PrintSelf( ostream& os, vtkIndent indent )
 		int numPoints = (*it)->getNumberOfPoints();
 		for (int i = 0; i < numPoints; i++)
 		{
-			const GeoLib::Point* point = (**it)[i];
+			const GeoLib::Point* point = (*it)->getPoint(i);
 			const double* coords = point->getCoords();
 			os << indent << "Point " << i << " (" << coords[0] << ", " << coords[1] <<
 			", " << coords[2] << ")\n";
@@ -115,7 +115,7 @@ int VtkPolylinesSource::RequestData( vtkInformation* request,
 		// Generate points
 		for (int i = 0; i < numPoints; i++)
 		{
-			const GeoLib::Point* point = (*(*_polylines)[j])[i];
+			const GeoLib::Point* point = (*_polylines)[j]->getPoint(i);
 			const double* coords = point->getCoords();
 			newPoints->InsertNextPoint(coords);
 		}
diff --git a/MathLib/AnalyticalGeometry.cpp b/MathLib/AnalyticalGeometry.cpp
index e708a4947e51f2b275855005912f3b2f8bd83c80..eb76424d28dee717bdc2b44bd53454855f71eb19 100644
--- a/MathLib/AnalyticalGeometry.cpp
+++ b/MathLib/AnalyticalGeometry.cpp
@@ -12,55 +12,56 @@
  *
  */
 
+#include "AnalyticalGeometry.h"
+
 #include <algorithm>
 #include <cmath>
 #include <cstdlib> // for exit
-#include <list>
-#include <limits>
 #include <fstream>
+#include <limits>
+#include <list>
 
-// Base
+// BaseLib
 #include "quicksort.h"
 
-// GEO
+// GeoLib
 #include "Polyline.h"
 #include "Triangle.h"
 
 // MathLib
-#include "MathTools.h"
-#include "AnalyticalGeometry.h"
-#include "LinAlg/Solvers/GaussAlgorithm.h"
 #include "LinAlg/Dense/Matrix.h" // for transformation matrix
+#include "LinAlg/Solvers/GaussAlgorithm.h"
+#include "MathTools.h"
 
-namespace MathLib {
-
-Orientation getOrientation (const double& p0_x, const double& p0_y,
-	const double& p1_x, const double& p1_y,
-	const double& p2_x, const double& p2_y)
+namespace MathLib
+{
+Orientation getOrientation(const double& p0_x, const double& p0_y, const double& p1_x,
+                           const double& p1_y, const double& p2_x, const double& p2_y)
 {
-	double h1 ((p1_x-p0_x)*(p2_y-p0_y));
-	double h2 ((p2_x-p0_x)*(p1_y-p0_y));
+	double h1((p1_x - p0_x) * (p2_y - p0_y));
+	double h2((p2_x - p0_x) * (p1_y - p0_y));
 
-	double tol (sqrt(std::numeric_limits<double>::min()));
-	if (fabs (h1-h2) <= tol * std::max (fabs(h1), fabs(h2)))
+	double tol(sqrt( std::numeric_limits<double>::min()));
+	if (fabs(h1 - h2) <= tol * std::max(fabs(h1), fabs(h2)))
 		return COLLINEAR;
-	if (h1-h2 > 0.0) return CCW;
+	if (h1 - h2 > 0.0)
+		return CCW;
 
 	return CW;
 }
 
-Orientation getOrientation (const GeoLib::Point* p0, const GeoLib::Point* p1, const GeoLib::Point* p2)
+Orientation getOrientation(const GeoLib::Point* p0, const GeoLib::Point* p1,
+                           const GeoLib::Point* p2)
 {
-		return getOrientation ((*p0)[0], (*p0)[1], (*p1)[0], (*p1)[1], (*p2)[0], (*p2)[1]);
+	return getOrientation((*p0)[0], (*p0)[1], (*p1)[0], (*p1)[1], (*p2)[0], (*p2)[1]);
 }
 
-bool lineSegmentIntersect (const GeoLib::Point& a, const GeoLib::Point& b,
-		const GeoLib::Point& c, const GeoLib::Point& d,
-		GeoLib::Point& s)
+bool lineSegmentIntersect(const GeoLib::Point& a, const GeoLib::Point& b, const GeoLib::Point& c,
+                          const GeoLib::Point& d, GeoLib::Point& s)
 {
-	Matrix<double> mat(2,2);
-	mat(0,0) = b[0] - a[0];
-	mat(1,0) = b[1] - a[1];
+	Matrix<double> mat(2, 2);
+	mat(0, 0) = b[0] - a[0];
+	mat(1, 0) = b[1] - a[1];
 	mat(0,1) = c[0] - d[0];
 	mat(1,1) = c[1] - d[1];
 
@@ -94,19 +95,22 @@ bool lineSegmentIntersect (const GeoLib::Point& a, const GeoLib::Point& b,
 		s[1] = a[1] + rhs[0] * (b[1] - a[1]);
 		s[2] = a[2] + rhs[0] * (b[2] - a[2]);
 		// check z component
-		double z0 (a[2] - d[2]), z1(rhs[0]*(b[2]-a[2]) + rhs[1]*(d[2]-c[2]));
+		double z0 (a[2] - d[2]), z1(rhs[0] * (b[2] - a[2]) + rhs[1] * (d[2] - c[2]));
 		delete [] rhs;
-		if (std::fabs (z0-z1) < eps)
+		if (std::fabs (z0 - z1) < eps)
 			return true;
 		else
 			return false;
-	} else delete [] rhs;
+	}
+	else
+		delete [] rhs;
 	return false;
 }
 
-bool lineSegmentsIntersect (const GeoLib::Polyline* ply, size_t &idx0, size_t &idx1, GeoLib::Point& intersection_pnt)
+bool lineSegmentsIntersect(const GeoLib::Polyline* ply, size_t &idx0, size_t &idx1,
+                           GeoLib::Point& intersection_pnt)
 {
-	size_t n_segs (ply->getNumberOfPoints() - 1);
+	size_t n_segs(ply->getNumberOfPoints() - 1);
 	/**
 	 * computing the intersections of all possible pairs of line segments of the given polyline
 	 * as follows:
@@ -114,10 +118,12 @@ bool lineSegmentsIntersect (const GeoLib::Polyline* ply, size_t &idx0, size_t &i
 	 * of the polyline and segment \f$s_2 = (C,B)\f$ defined by \f$j\f$-th and
 	 * \f$j+1\f$-st point of the polyline, \f$j>k+1\f$
 	 */
-	for (size_t k(0); k<n_segs-2; k++) {
-		for (size_t j(k+2); j<n_segs; j++) {
-			if (k!=0 || j<n_segs-1) {
-				if (lineSegmentIntersect (*(*ply)[k], *(*ply)[k+1], *(*ply)[j], *(*ply)[j+1], intersection_pnt)) {
+	for (size_t k(0); k < n_segs - 2; k++) {
+		for (size_t j(k + 2); j < n_segs; j++) {
+			if (k != 0 || j < n_segs - 1) {
+				if (lineSegmentIntersect(*(ply->getPoint(k)), *(ply->getPoint(k + 1)),
+				                         *(ply->getPoint(j)), *(ply->getPoint(j + 1)),
+				                         intersection_pnt)) {
 					idx0 = k;
 					idx1 = j;
 					return true;
@@ -128,54 +134,56 @@ bool lineSegmentsIntersect (const GeoLib::Polyline* ply, size_t &idx0, size_t &i
 	return false;
 }
 
-bool isPointInTriangle (const double p[3], const double a[3], const double b[3], const double c[3])
+bool isPointInTriangle(const double p[3], const double a[3], const double b[3], const double c[3])
 {
 	// criterion: p-b = u0 * (b - a) + u1 * (b - c); 0 <= u0, u1 <= 1, u0+u1 <= 1
-	MathLib::Matrix<double> mat (2,2);
-	mat(0,0) = a[0] - b[0];
-	mat(0,1) = c[0] - b[0];
-	mat(1,0) = a[1] - b[1];
-	mat(1,1) = c[1] - b[1];
-	double rhs[2] = {p[0]-b[0], p[1]-b[1]};
+	MathLib::Matrix<double> mat(2, 2);
+	mat(0, 0) = a[0] - b[0];
+	mat(0, 1) = c[0] - b[0];
+	mat(1, 0) = a[1] - b[1];
+	mat(1, 1) = c[1] - b[1];
+	double rhs[2] = { p[0] - b[0], p[1] - b[1] };
 
-	MathLib::GaussAlgorithm gauss (mat);
-	gauss.execute (rhs);
+	MathLib::GaussAlgorithm gauss(mat);
+	gauss.execute(rhs);
 
-	if (0 <= rhs[0] && rhs[0] <= 1 && 0 <= rhs[1] && rhs[1] <= 1 && rhs[0] + rhs[1] <= 1) return true;
+	if (0 <= rhs[0] && rhs[0] <= 1 && 0 <= rhs[1] && rhs[1] <= 1 && rhs[0] + rhs[1] <= 1)
+		return true;
 	return false;
 }
 
-bool isPointInTriangle (const GeoLib::Point* p,
-		const GeoLib::Point* a, const GeoLib::Point* b, const GeoLib::Point* c)
+bool isPointInTriangle(const GeoLib::Point* p, const GeoLib::Point* a, const GeoLib::Point* b,
+                       const GeoLib::Point* c)
 {
-	return isPointInTriangle (p->getCoords(), a->getCoords(), b->getCoords(), c->getCoords());
+	return isPointInTriangle(p->getCoords(), a->getCoords(), b->getCoords(), c->getCoords());
 }
 
-double getOrientedTriArea(GeoLib::Point const& a, GeoLib::Point const& b,
-				GeoLib::Point const& c)
+double getOrientedTriArea(GeoLib::Point const& a, GeoLib::Point const& b, GeoLib::Point const& c)
 {
-	const double u[3] = {c[0]-a[0], c[1]-a[1], c[2]-a[2]};
-	const double v[3] = {b[0]-a[0], b[1]-a[1], b[2]-a[2]};
+	const double u[3] = { c[0] - a[0], c[1] - a[1], c[2] - a[2] };
+	const double v[3] = { b[0] - a[0], b[1] - a[1], b[2] - a[2] };
 	double w[3];
 	MathLib::crossProd(u, v, w);
-	return 0.5 * sqrt(MathLib::scpr<double,3>(w,w));
+	return 0.5 * sqrt(MathLib::scpr<double, 3>(w, w));
 }
 
-bool isPointInTriangle(GeoLib::Point const& p, GeoLib::Point const& a,
-				GeoLib::Point const& b, GeoLib::Point const& c,
-				double eps)
+bool isPointInTriangle(GeoLib::Point const& p, GeoLib::Point const& a, GeoLib::Point const& b,
+                       GeoLib::Point const& c, double eps)
 {
 	const unsigned dim(3);
-	MathLib::Matrix<double> m (dim,dim);
-	for (unsigned i(0); i<dim; i++) m(i,0) = b[i] - a[i];
-	for (unsigned i(0); i<dim; i++) m(i,1) = c[i] - a[i];
-	for (unsigned i(0); i<dim; i++) m(i,2) = p[i] - a[i];
+	MathLib::Matrix<double> m(dim, dim);
+	for (unsigned i(0); i < dim; i++)
+		m(i, 0) = b[i] - a[i];
+	for (unsigned i(0); i < dim; i++)
+		m(i, 1) = c[i] - a[i];
+	for (unsigned i(0); i < dim; i++)
+		m(i, 2) = p[i] - a[i];
 
 	// point p is in the same plane as the triangle if and only if
 	// the following determinate of the 3x3 matrix equals zero (up to an eps)
-	double det3x3(m(0,0) * (m(1,1) * m(2,2) - m(2,1) * m(1,2))
-					- m(1,0) * (m(2,1) * m(0,2) - m(0,1) * m(2,2))
-					+ m(2,0) * (m(0,1) * m(1,2) - m(1,1) * m(0,2)));
+	double det3x3(m(0, 0) * (m(1, 1) * m(2, 2) - m(2, 1) * m(1, 2))
+	              - m(1, 0) * (m(2, 1) * m(0, 2) - m(0, 1) * m(2, 2))
+	              + m(2, 0) * (m(0, 1) * m(1, 2) - m(1, 1) * m(0, 2)));
 	if (fabs(det3x3) > eps)
 		return false;
 
@@ -190,19 +198,18 @@ bool isPointInTriangle(GeoLib::Point const& p, GeoLib::Point const& a,
 }
 
 // NewellPlane from book Real-Time Collision detection p. 494
-void getNewellPlane(const std::vector<GeoLib::Point*>& pnts, Vector &plane_normal,
-		double& d)
+void getNewellPlane(const std::vector<GeoLib::Point*>& pnts, Vector &plane_normal, double& d)
 {
 	d = 0;
 	Vector centroid;
-	size_t n_pnts (pnts.size());
+	size_t n_pnts(pnts.size());
 	for (size_t i(n_pnts - 1), j(0); j < n_pnts; i = j, j++) {
 		plane_normal[0] += ((*(pnts[i]))[1] - (*(pnts[j]))[1])
-				* ((*(pnts[i]))[2] + (*(pnts[j]))[2]); // projection on yz
+		                   * ((*(pnts[i]))[2] + (*(pnts[j]))[2]); // projection on yz
 		plane_normal[1] += ((*(pnts[i]))[2] - (*(pnts[j]))[2])
-				* ((*(pnts[i]))[0] + (*(pnts[j]))[0]); // projection on xz
+		                   * ((*(pnts[i]))[0] + (*(pnts[j]))[0]); // projection on xz
 		plane_normal[2] += ((*(pnts[i]))[0] - (*(pnts[j]))[0])
-				* ((*(pnts[i]))[1] + (*(pnts[j]))[1]); // projection on xy
+		                   * ((*(pnts[i]))[1] + (*(pnts[j]))[1]); // projection on xy
 
 		centroid += *(pnts[j]);
 	}
@@ -211,10 +218,9 @@ void getNewellPlane(const std::vector<GeoLib::Point*>& pnts, Vector &plane_norma
 	d = centroid.Dot(plane_normal) / n_pnts;
 }
 
-void rotatePointsToXY(Vector &plane_normal,
-                      std::vector<GeoLib::Point*> &pnts)
+void rotatePointsToXY(Vector &plane_normal, std::vector<GeoLib::Point*> &pnts)
 {
-	double small_value (sqrt (std::numeric_limits<double>::min()));
+	double small_value(sqrt( std::numeric_limits<double>::min()));
 	if (fabs(plane_normal[0]) < small_value && fabs(plane_normal[1]) < small_value)
 		return;
 
@@ -222,16 +228,16 @@ void rotatePointsToXY(Vector &plane_normal,
 	computeRotationMatrixToXY(plane_normal, rot_mat);
 	rotatePoints(rot_mat, pnts);
 
-	double* tmp (rot_mat * plane_normal.getCoords());
+	double* tmp(rot_mat * plane_normal.getCoords());
 	for (std::size_t j(0); j < 3; j++)
 		plane_normal[j] = tmp[j];
 
-	delete [] tmp;
+	delete[] tmp;
 }
 
 void rotatePointsToXZ(Vector &n, std::vector<GeoLib::Point*> &pnts)
 {
-	double small_value (sqrt (std::numeric_limits<double>::min()));
+	double small_value(sqrt( std::numeric_limits<double>::min()));
 	if (fabs(n[0]) < small_value && fabs(n[1]) < small_value)
 		return;
 
@@ -246,22 +252,22 @@ void rotatePointsToXZ(Vector &n, std::vector<GeoLib::Point*> &pnts)
 	Matrix<double> rot_mat(3, 3);
 	// calc rotation matrix
 	rot_mat(0, 0) = n[1] * h1;
-	rot_mat(0, 1) = - n[0] * h1;
+	rot_mat(0, 1) = -n[0] * h1;
 	rot_mat(0, 2) = 0.0;
 	rot_mat(1, 0) = n[0] * h2;
 	rot_mat(1, 1) = n[1] * h2;
 	rot_mat(1, 2) = n[2] * h2;
 	rot_mat(2, 0) = n[0] * n[2] * h1 * h2;
 	rot_mat(2, 1) = n[1] * n[2] * h1 * h2;
-	rot_mat(2, 2) = - sqrt(h0) * h2;
+	rot_mat(2, 2) = -sqrt(h0) * h2;
 
 	rotatePoints(rot_mat, pnts);
 
-	double *tmp (rot_mat * n.getCoords());
+	double *tmp(rot_mat * n.getCoords());
 	for (std::size_t j(0); j < 3; j++)
 		n[j] = tmp[j];
 
-	delete [] tmp;
+	delete[] tmp;
 }
 
 void computeRotationMatrixToXY(Vector const& plane_normal, Matrix<double> & rot_mat)
@@ -298,5 +304,4 @@ void rotatePoints(Matrix<double> const& rot_mat, std::vector<GeoLib::Point*> &pn
 		delete [] tmp;
 	}
 }
-
 } // end namespace MathLib
diff --git a/Tests/GeoLib/TestPolyline.cpp b/Tests/GeoLib/TestPolyline.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..91ac9d87b994896a5802d5d007e3d34db7f37728
--- /dev/null
+++ b/Tests/GeoLib/TestPolyline.cpp
@@ -0,0 +1,108 @@
+/**
+ * @file TestPolyline.cpp
+ * @author Thomas Fischer
+ * @date Jan 7, 2013
+ *
+ * @copyright
+ * Copyright (c) 2013, 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 <ctime>
+
+#include "Polyline.h"
+
+using namespace GeoLib;
+
+TEST(GeoLib, PolylineTest)
+{
+	std::vector<Point*> ply_pnts;
+	Polyline ply(ply_pnts);
+
+	// checking properties of empty polyline
+	ASSERT_EQ(ply.getNumberOfPoints(), 0);
+	ASSERT_FALSE(ply.isClosed());
+	ASSERT_FALSE(ply.isPointIDInPolyline(0));
+	ASSERT_EQ(ply.getLength(0), 0.0);
+
+	ply_pnts.push_back(new Point(0.0, 0.0, 0.0));
+	ply.addPoint(0);
+	// checking properties of polyline with one point
+	ASSERT_EQ(ply.getNumberOfPoints(), 1);
+	ASSERT_FALSE(ply.isClosed());
+	ASSERT_TRUE(ply.isPointIDInPolyline(0));
+	ASSERT_EQ(ply.getLength(0), 0.0);
+
+	ply_pnts.push_back(new Point(1.0, 0.0, 0.0));
+	ply.addPoint(1);
+	// checking properties of polyline with two points
+	ASSERT_EQ(ply.getNumberOfPoints(), 2);
+	ASSERT_FALSE(ply.isClosed());
+	ASSERT_TRUE(ply.isPointIDInPolyline(1));
+	ASSERT_EQ(ply.getLength(1), 1);
+
+	// checking properties of polyline with two points
+	ply_pnts.push_back(new Point(0.5, 0.5, 0.0));
+	ply.addPoint(2);
+	ASSERT_EQ(ply.getNumberOfPoints(), 3);
+	ASSERT_FALSE(ply.isClosed());
+	ASSERT_TRUE(ply.isPointIDInPolyline(2));
+	ASSERT_TRUE(fabs(ply.getLength(2) - (1.0 + sqrt(0.5))) < std::numeric_limits<double>::min());
+
+	// checking remove
+	ply.removePoint(1);
+	ASSERT_EQ(ply.getNumberOfPoints(), 2);
+	ASSERT_FALSE(ply.isClosed());
+	ASSERT_FALSE(ply.isPointIDInPolyline(1));
+	ASSERT_TRUE(fabs(ply.getLength(1) - sqrt(0.5)) < std::numeric_limits<double>::min());
+
+	// inserting point in the middle
+	ply.insertPoint(1,1);
+	ASSERT_EQ(ply.getNumberOfPoints(), 3);
+	ASSERT_FALSE(ply.isClosed());
+	ASSERT_TRUE(ply.isPointIDInPolyline(2));
+	ASSERT_TRUE(fabs(ply.getLength(2) - (1.0 + sqrt(0.5))) < std::numeric_limits<double>::epsilon());
+
+	// inserting point at the end
+	ply_pnts.push_back(new GeoLib::Point(1.0, 0.5, 0.0));
+	ply.insertPoint(3,3);
+	ASSERT_EQ(ply.getNumberOfPoints(), 4);
+	ASSERT_FALSE(ply.isClosed());
+	ASSERT_TRUE(ply.isPointIDInPolyline(3));
+	ASSERT_TRUE(fabs(ply.getLength(3) - (1.0 + sqrt(0.5) + 0.5)) < std::numeric_limits<double>::epsilon());
+
+	// inserting point at the beginning
+	ply_pnts.push_back(new GeoLib::Point(-1.0, 0.0, 0.0));
+	ply.insertPoint(0,4);
+	ASSERT_EQ(ply.getNumberOfPoints(), 5);
+	ASSERT_FALSE(ply.isClosed());
+	ASSERT_TRUE(ply.isPointIDInPolyline(4));
+	ASSERT_TRUE(fabs(ply.getLength(4) - (1.0 + 1.0 + sqrt(0.5) + 0.5)) <
+	            std::numeric_limits<double>::epsilon());
+
+	// inserting point in the middle
+	ply_pnts.push_back(new GeoLib::Point(0.0, 0.5, 0.0));
+	ply.insertPoint(2,5);
+	ASSERT_EQ(ply.getNumberOfPoints(), 6);
+	ASSERT_FALSE(ply.isClosed());
+	ASSERT_TRUE(ply.isPointIDInPolyline(5));
+	ASSERT_TRUE(fabs(ply.getLength(5) - (1.0 + 0.5 + sqrt(1.25) + sqrt(0.5) + 0.5)) < std::numeric_limits<double>::epsilon());
+
+	// close polyline
+	ply.closePolyline();
+	ASSERT_EQ(ply.getNumberOfPoints(), 7);
+	ASSERT_TRUE(ply.isClosed());
+
+	// remove last point -> polyline is not closed!
+	ply.removePoint(6);
+	ASSERT_EQ(ply.getNumberOfPoints(), 6);
+	ASSERT_FALSE(ply.isClosed());
+	ply.closePolyline();
+	ASSERT_TRUE(ply.isClosed());
+
+	for (std::size_t k(0); k < ply_pnts.size(); ++k)
+		delete ply_pnts[k];
+}
diff --git a/Utils/SimpleMeshCreation/createMeshElemPropertiesFromASCRaster.cpp b/Utils/SimpleMeshCreation/createMeshElemPropertiesFromASCRaster.cpp
index d4eb23402531d0e3e32e555a9297b1f8b667b1bc..85ee4aedb4bd13ba117b3470863245d60a3e2911 100644
--- a/Utils/SimpleMeshCreation/createMeshElemPropertiesFromASCRaster.cpp
+++ b/Utils/SimpleMeshCreation/createMeshElemPropertiesFromASCRaster.cpp
@@ -12,10 +12,15 @@
  *
  */
 
+// stl
+#include <numeric>
+
 // BaseLib
 #include "tclap/CmdLine.h"
+
 // ThirdParty/logog
 #include "logog/include/logog.hpp"
+
 // BaseLib
 #include "LogogSimpleFormatter.h"
 #include "StringTools.h"
@@ -32,6 +37,9 @@
 // Gui/VtkVis
 #include "VtkMeshConverter.h"
 
+// MathLib
+#include "MathTools.h"
+
 // MeshLib
 #include "ConvertRasterToMesh.h"
 #include "Elements/Element.h"
@@ -111,8 +119,10 @@ int main (int argc, char* argv[])
 		raster->refineRaster(refinement_arg.getValue());
 		if (refinement_raster_output_arg.getValue()) {
 			// write new asc file
-			std::string new_raster_fname (BaseLib::dropFileExtension(raster_arg.getValue()));
-			new_raster_fname += "-" + BaseLib::number2str(raster->getNRows()) + "x" + BaseLib::number2str(raster->getNCols()) + ".asc";
+			std::string new_raster_fname (BaseLib::dropFileExtension(
+			                                      raster_arg.getValue()));
+			new_raster_fname += "-" + BaseLib::number2str(raster->getNRows()) + "x" +
+			                    BaseLib::number2str(raster->getNCols()) + ".asc";
 			std::ofstream out(new_raster_fname);
 			raster->writeRasterAsASC(out);
 			out.close();
@@ -121,110 +131,104 @@ int main (int argc, char* argv[])
 
 	// put raster data in a std::vector
 	GeoLib::Raster::const_iterator raster_it(raster->begin());
-	unsigned n_cols(raster->getNCols()), n_rows(raster->getNRows());
-	std::vector<double> src_properties(n_cols * n_rows);
+	std::size_t n_cols(raster->getNCols()), n_rows(raster->getNRows());
+	std::size_t size(n_cols * n_rows);
+	std::vector<double> src_properties(size);
 	for (unsigned row(0); row<n_rows; row++) {
 		for (unsigned col(0); col<n_cols; col++) {
-			src_properties[row*n_cols+col] = *raster_it;
+			src_properties[row * n_cols + col] = *raster_it;
 			++raster_it;
 		}
 	}
 
 	{
-		double src_mean_value(src_properties[0]);
-		for (size_t k(1); k < n_cols * n_rows; k++)
-			src_mean_value += src_properties[k];
-		src_mean_value /= n_cols*n_rows;
-		std::cout << "mean value of source: " << src_mean_value << std::endl;
-
-		double src_varianz(MathLib::fastpow(src_properties[0] - src_mean_value, 2));
-		for (size_t k(1); k<n_cols*n_rows; k++) {
-			src_varianz += MathLib::fastpow(src_properties[k] - src_mean_value, 2);
+		const double mu(std::accumulate(src_properties.begin(), src_properties.end(), 0) / size);
+		INFO("Mean value of source: %f.", mu);
+
+		double src_variance(MathLib::fastpow(src_properties[0] - mu, 2));
+		for (std::size_t k(1); k<size; k++) {
+			src_variance += MathLib::fastpow(src_properties[k] - mu, 2);
 		}
-		src_varianz /= n_cols * n_rows;
-		std::cout << "variance of source: " << src_varianz << std::endl;
+		src_variance /= size;
+		INFO("Variance of source: %f.", src_variance);
 	}
 
 	MeshLib::Mesh* src_mesh(MeshLib::ConvertRasterToMesh(*raster, MshElemType::QUAD,
 					MeshLib::UseIntensityAs::MATERIAL).execute());
 
-	std::vector<size_t> src_perm(n_cols * n_rows);
-	for (size_t k(0); k < n_cols * n_rows; k++)
-		src_perm[k] = k;
-	BaseLib::Quicksort<double, std::size_t>(src_properties, 0, n_cols * n_rows, src_perm);
+	std::vector<std::size_t> src_perm(size);
+	std::iota(src_perm.begin(), src_perm.end(), 0);
+	BaseLib::Quicksort<double>(src_properties, 0, size, src_perm);
 
 	// compress the property data structure
-	const size_t mat_map_size(src_properties.size());
-	std::vector<size_t> mat_map(mat_map_size);
+	const std::size_t mat_map_size(src_properties.size());
+	std::vector<std::size_t> mat_map(mat_map_size);
 	mat_map[0] = 0;
-	size_t n_mat(1);
-	for (size_t k(1); k<mat_map_size; ++k) {
+	std::size_t n_mat(1);
+	for (std::size_t k(1); k<mat_map_size; ++k) {
 		if (std::fabs(src_properties[k] - src_properties[k-1]) > std::numeric_limits<double>::epsilon()) {
-			mat_map[k] = mat_map[k-1]+1;
+			mat_map[k] = mat_map[k - 1] + 1;
 			n_mat++;
 		} else
-			mat_map[k] = mat_map[k-1];
+			mat_map[k] = mat_map[k - 1];
 	}
 	std::vector<double> compressed_src_properties(n_mat);
 	compressed_src_properties[0] = src_properties[0];
-	for (size_t k(1), id(1); k<mat_map_size; ++k) {
+	for (std::size_t k(1), id(1); k<mat_map_size; ++k) {
 		if (std::fabs(src_properties[k] - src_properties[k-1]) > std::numeric_limits<double>::epsilon()) {
 			compressed_src_properties[id] = src_properties[k];
 			id++;
 		}
 	}
-	compressed_src_properties[n_mat-1] = src_properties[mat_map_size-1];
+	compressed_src_properties[n_mat - 1] = src_properties[mat_map_size - 1];
 
 	// reset materials in source mesh
-	const size_t n_mesh_elements(src_mesh->getNElements());
-	for (size_t k(0); k<n_mesh_elements; k++) {
+	const std::size_t n_mesh_elements(src_mesh->getNElements());
+	for (std::size_t k(0); k<n_mesh_elements; k++) {
 		const_cast<MeshLib::Element*>(src_mesh->getElement(src_perm[k]))->setValue(mat_map[k]);
 	}
 
 	// do the interpolation
-	MeshLib::Mesh2MeshPropertyInterpolation mesh_interpolation(src_mesh, &compressed_src_properties);
+	MeshLib::Mesh2MeshPropertyInterpolation mesh_interpolation(src_mesh,
+	                                                           &compressed_src_properties);
 	std::vector<double> dest_properties(dest_mesh->getNElements());
-	mesh_interpolation.setPropertiesForMesh(const_cast<MeshLib::Mesh*>(dest_mesh), dest_properties);
+	mesh_interpolation.setPropertiesForMesh(const_cast<MeshLib::Mesh*>(dest_mesh),
+	                                        dest_properties);
 
-	const size_t n_dest_mesh_elements(dest_mesh->getNElements());
+	const std::size_t n_dest_mesh_elements(dest_mesh->getNElements());
 
 	{ // write property file
 		std::string property_fname(mapping_arg.getValue());
 		std::ofstream property_out(property_fname.c_str());
-		if (! property_out) {
-			std::cerr << "could not open file " << property_fname << "PropertyMapping" << std::endl;
+		if (!property_out)
+		{
+			ERR("Could not open file %s for writing the mapping.", property_fname.c_str());
 			return -1;
 		}
 
-		for (size_t k(0); k<n_dest_mesh_elements; k++) {
+		for (std::size_t k(0); k < n_dest_mesh_elements; k++)
 			property_out << k << " " << dest_properties[k] << "\n";
-		}
 		property_out.close();
 	}
 
 	{
-		double mu(dest_properties[0]);
-		for (size_t k(1); k<n_dest_mesh_elements; k++) {
-			mu += dest_properties[k];
-		}
-		mu /= n_dest_mesh_elements;
-		std::cout << "mean value of destination: " << mu << std::endl;
+		const double mu(std::accumulate(dest_properties.begin(), dest_properties.end(), 0.0) / n_dest_mesh_elements);
+		INFO("Mean value of destination: %f.", mu);
 
 		double sigma_q(MathLib::fastpow(dest_properties[0] - mu, 2));
-		for (size_t k(1); k<n_dest_mesh_elements; k++) {
+		for (std::size_t k(1); k < n_dest_mesh_elements; k++)
 			sigma_q += MathLib::fastpow(dest_properties[k] - mu, 2);
-		}
 		sigma_q /= n_dest_mesh_elements;
-		std::cout << "variance of destination: " << sigma_q << std::endl;
+		INFO("Variance of destination: %f.", sigma_q);
 	}
 
 	if (! out_mesh_arg.getValue().empty()) {
-		std::vector<size_t> dest_perm(n_dest_mesh_elements);
-		for (size_t k(0); k<n_dest_mesh_elements; k++) dest_perm[k] = k;
-		BaseLib::Quicksort<double, std::size_t>(dest_properties, 0, n_dest_mesh_elements, dest_perm);
+		std::vector<std::size_t> dest_perm(n_dest_mesh_elements);
+		std::iota(dest_perm.begin(), dest_perm.end(), 0);
+		BaseLib::Quicksort<double>(dest_properties, 0, n_dest_mesh_elements, dest_perm);
 
 		// reset materials in destination mesh
-		for (size_t k(0); k<n_dest_mesh_elements; k++) {
+		for (std::size_t k(0); k<n_dest_mesh_elements; k++) {
 			const_cast<MeshLib::Element*>(dest_mesh->getElement(dest_perm[k]))->setValue(k);
 		}