diff --git a/FileIO/MeshIO/GMSHAdaptiveMeshDensity.cpp b/FileIO/MeshIO/GMSHAdaptiveMeshDensity.cpp
index 0061f3cad45cf14907d3ceb31cfc239577a7a36c..13b28ceacb3adc1605a8e12a52d042fe41a1c22c 100644
--- a/FileIO/MeshIO/GMSHAdaptiveMeshDensity.cpp
+++ b/FileIO/MeshIO/GMSHAdaptiveMeshDensity.cpp
@@ -12,14 +12,17 @@
 
 #include <list>
 
+// ThirdParty/logog
+#include "logog/include/logog.hpp"
+
 // FileIO
 #include "MeshIO/GMSHAdaptiveMeshDensity.h"
 
 // GeoLib
 #include "Polygon.h"
 
-namespace FileIO {
-
+namespace FileIO
+{
 GMSHAdaptiveMeshDensity::GMSHAdaptiveMeshDensity(double pnt_density, double station_density,
 				size_t max_pnts_per_leaf) :
 	_pnt_density(pnt_density), _station_density(station_density),
@@ -35,32 +38,26 @@ GMSHAdaptiveMeshDensity::~GMSHAdaptiveMeshDensity()
 void GMSHAdaptiveMeshDensity::init(std::vector<GeoLib::Point const*> const& pnts)
 {
 	// *** QuadTree - determining bounding box
-#ifndef NDEBUG
-	std::cout << "[GMSHAdaptiveMeshDensity::init]" << std::endl;
-	std::cout << "\tcomputing axis aligned bounding box (2D) for quadtree ... " << std::flush;
-#endif
+	DBUG("GMSHAdaptiveMeshDensity::init(): computing axis aligned bounding box (2D) for quadtree.");
+
 	GeoLib::Point min(pnts[0]->getCoords()), max(pnts[0]->getCoords());
 	size_t n_pnts(pnts.size());
 	for (size_t k(1); k<n_pnts; k++) {
-		for (size_t j(0); j<2; j++)
-			if ((*(pnts[k]))[j] < min[j]) min[j] = (*(pnts[k]))[j];
-		for (size_t j(0); j<2; j++)
-			if ((*(pnts[k]))[j] > max[j]) max[j] = (*(pnts[k]))[j];
+		for (size_t j(0); j < 2; j++)
+			if ((*(pnts[k]))[j] < min[j])
+				min[j] = (*(pnts[k]))[j];
+		for (size_t j(0); j < 2; j++)
+			if ((*(pnts[k]))[j] > max[j])
+				max[j] = (*(pnts[k]))[j];
 	}
 	min[2] = 0.0;
 	max[2] = 0.0;
-#ifndef NDEBUG
-	std::cout << "ok" << std::endl;
-#endif
+	DBUG("GMSHAdaptiveMeshDensity::init(): \tok");
 
 	// *** QuadTree - create object
-#ifndef NDEBUG
-	std::cout << "\tcreating quadtree ... " << std::flush;
-#endif
+	DBUG("GMSHAdaptiveMeshDensity::init(): Creating quadtree.");
 	_quad_tree = new GeoLib::QuadTree<GeoLib::Point> (min, max, _max_pnts_per_leaf);
-#ifndef NDEBUG
-	std::cout << "ok" << std::endl;
-#endif
+	DBUG("GMSHAdaptiveMeshDensity::init(): \tok.");
 
 	// *** QuadTree - insert points
 	addPoints(pnts);
@@ -70,33 +67,29 @@ void GMSHAdaptiveMeshDensity::addPoints(std::vector<GeoLib::Point const*> const&
 {
 	// *** QuadTree - insert points
 	const size_t n_pnts(pnts.size());
-#ifndef NDEBUG
-	std::cout << "\tinserting " << n_pnts << " points into quadtree ... " <<
-	std::flush;
-#endif
+	DBUG("GMSHAdaptiveMeshDensity::addPoints(): Inserting %d points into quadtree.", n_pnts);
 	for (size_t k(0); k < n_pnts; k++)
 		_quad_tree->addPoint(pnts[k]);
-#ifndef NDEBUG
-	std::cout << "ok" << std::endl;
-#endif
+	DBUG("GMSHAdaptiveMeshDensity::addPoints(): \tok.");
 	_quad_tree->balance();
 }
 
-double GMSHAdaptiveMeshDensity::getMeshDensityAtPoint(GeoLib::Point const*const pnt) const
+double GMSHAdaptiveMeshDensity::getMeshDensityAtPoint(GeoLib::Point const* const pnt) const
 {
 	GeoLib::Point ll, ur;
 	_quad_tree->getLeaf(*pnt, ll, ur);
 	return _pnt_density * (ur[0] - ll[0]);
 }
 
-double GMSHAdaptiveMeshDensity::getMeshDensityAtStation(GeoLib::Point const*const pnt) const
+double GMSHAdaptiveMeshDensity::getMeshDensityAtStation(GeoLib::Point const* const pnt) const
 {
 	GeoLib::Point ll, ur;
 	_quad_tree->getLeaf(*pnt, ll, ur);
 	return (_station_density * (ur[0] - ll[0]));
 }
 
-void GMSHAdaptiveMeshDensity::getSteinerPoints (std::vector<GeoLib::Point*> & pnts, size_t additional_levels) const
+void GMSHAdaptiveMeshDensity::getSteinerPoints (std::vector<GeoLib::Point*> & pnts,
+                                                size_t additional_levels) const
 {
 	// get Steiner points
 	size_t max_depth(0);
@@ -144,11 +137,11 @@ void GMSHAdaptiveMeshDensity::getQuadTreeGeometry(std::vector<GeoLib::Point*> &p
 		pnts.push_back(ur);
 		pnts.push_back(new GeoLib::Point((*ll)[0], (*ur)[1], 0.0));
 		plys.push_back(new GeoLib::Polyline(pnts));
-		plys[plys.size()-1]->addPoint(pnt_offset);
-		plys[plys.size()-1]->addPoint(pnt_offset+1);
-		plys[plys.size()-1]->addPoint(pnt_offset+2);
-		plys[plys.size()-1]->addPoint(pnt_offset+3);
-		plys[plys.size()-1]->addPoint(pnt_offset);
+		plys[plys.size() - 1]->addPoint(pnt_offset);
+		plys[plys.size() - 1]->addPoint(pnt_offset + 1);
+		plys[plys.size() - 1]->addPoint(pnt_offset + 2);
+		plys[plys.size() - 1]->addPoint(pnt_offset + 3);
+		plys[plys.size() - 1]->addPoint(pnt_offset);
 	}
 }
 #endif
diff --git a/FileIO/MeshIO/GMSHInterface.cpp b/FileIO/MeshIO/GMSHInterface.cpp
index 66053345474aca57fe1201f6af905354c99d8426..da88f0049d922c307de11c80c69a3f8da01e8d96 100644
--- a/FileIO/MeshIO/GMSHInterface.cpp
+++ b/FileIO/MeshIO/GMSHInterface.cpp
@@ -1,51 +1,60 @@
 /**
- * Copyright (c) 2012, OpenGeoSys Community (http://www.opengeosys.net)
+ * @copyright
+ * Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.net)
  *            Distributed under a Modified BSD License.
  *              See accompanying file LICENSE.txt or
  *              http://www.opengeosys.org/LICENSE.txt
  *
- * \file GMSHInterface.cpp
- *
- *  Created on 2010-04-29 by Thomas Fischer
+ * @file GMSHInterface.cpp
+ * @date 2010-04-29
+ * @author Thomas Fischer
  */
 
 #include <fstream>
 #include <vector>
 
+// ThirdParty/logog
+#include "logog/include/logog.hpp"
+
 // BaseLib
-#include "StringTools.h"
-#include "Configure.h"
 #include "BuildInfo.h"
+#include "Configure.h"
 #include "FileTools.h"
+#include "StringTools.h"
 
 // FileIO
+#include "GMSHAdaptiveMeshDensity.h"
+#include "GMSHFixedMeshDensity.h"
 #include "GMSHInterface.h"
 #include "GMSHNoMeshDensity.h"
-#include "GMSHFixedMeshDensity.h"
-#include "GMSHAdaptiveMeshDensity.h"
 
 // GeoLib
 #include "Point.h"
 #include "Polygon.h"
 #include "Polyline.h"
-#include "QuadTree.h"
 #include "PolylineWithSegmentMarker.h"
+#include "QuadTree.h"
 
 // MSH
-#include "Mesh.h"
-#include "Node.h"
 #include "Elements/Edge.h"
-#include "Elements/Tri.h"
-#include "Elements/Quad.h"
-#include "Elements/Tet.h"
 #include "Elements/Hex.h"
-#include "Elements/Pyramid.h"
 #include "Elements/Prism.h"
+#include "Elements/Pyramid.h"
+#include "Elements/Quad.h"
+#include "Elements/Tet.h"
+#include "Elements/Tri.h"
+#include "Mesh.h"
+#include "Node.h"
 
-namespace FileIO {
-GMSHInterface::GMSHInterface(GeoLib::GEOObjects & geo_objs, bool include_stations_as_constraints,
-				GMSH::MeshDensityAlgorithm mesh_density_algorithm, double param1, double param2,
-				size_t param3, std::vector<std::string>& selected_geometries) :
+namespace FileIO
+{
+GMSHInterface::GMSHInterface(GeoLib::GEOObjects & geo_objs,
+                             bool include_stations_as_constraints,
+                             GMSH::MeshDensityAlgorithm mesh_density_algorithm,
+                             double param1,
+                             double param2,
+                             size_t param3,
+                             std::vector<std::string>& selected_geometries) :
 	_n_lines(0), _n_plane_sfc(0), _geo_objs(geo_objs), _selected_geometries(selected_geometries),
 	_include_stations_as_constraints(include_stations_as_constraints)
 {
@@ -67,7 +76,7 @@ bool GMSHInterface::isGMSHMeshFile(const std::string& fname)
 	std::ifstream input(fname.c_str());
 
 	if (!input) {
-		std::cerr << "GMSHInterface::isGMSHMeshFile could not open file " << fname << std::endl;
+		ERR("GMSHInterface::isGMSHMeshFile(): Could not open file %s.", fname.c_str());
 		return false;
 	}
 
@@ -78,7 +87,8 @@ bool GMSHInterface::isGMSHMeshFile(const std::string& fname)
 		std::string version;
 		getline(input, version);
 		getline(input, version);
-		std::cerr << "found GMSH mesh file version: " << version << std::endl;
+		INFO("GMSHInterface::isGMSHMeshFile(): Found GMSH mesh file version: %s.",
+		     version.c_str());
 		input.close();
 		return true;
 	}
@@ -135,74 +145,74 @@ MeshLib::Mesh* GMSHInterface::readGMSHMesh(std::string const& fname)
 
 				switch (type)
 				{
-					case 1: {
-						readNodeIDs(in, 2, node_ids, id_map);
-						// edge_nodes array will be deleted from Edge object
-						MeshLib::Node** edge_nodes = new MeshLib::Node*[2];
-						edge_nodes[0] = nodes[node_ids[0]];
-						edge_nodes[1] = nodes[node_ids[1]];
-						elem = new MeshLib::Edge(edge_nodes, 0);
-						break;
-					}
-					case 2: {
-						readNodeIDs(in, 3, node_ids, id_map);
-						MeshLib::Node** tri_nodes = new MeshLib::Node*[3];
-						tri_nodes[0] = nodes[node_ids[2]];
-						tri_nodes[1] = nodes[node_ids[1]];
-						tri_nodes[2] = nodes[node_ids[0]];
-						elem = new MeshLib::Tri(tri_nodes, mat_id);
-						break;
-					}
-					case 3: {
-						readNodeIDs(in, 4, node_ids, id_map);
-						MeshLib::Node** quad_nodes = new MeshLib::Node*[4];
-						for (unsigned k(0); k<4; k++)
-							quad_nodes[k] = nodes[node_ids[k]];
-						elem = new MeshLib::Quad(quad_nodes, mat_id);
-						break;
-					}
-					case 4: {
-						readNodeIDs(in, 4, node_ids, id_map);
-						MeshLib::Node** tet_nodes = new MeshLib::Node*[5];
-						for (unsigned k(0); k<4; k++)
-							tet_nodes[k] = nodes[node_ids[k]];
-						elem = new MeshLib::Tet(tet_nodes, mat_id);
-						break;
-					}
-					case 5: {
-						readNodeIDs(in, 8, node_ids, id_map);
-						MeshLib::Node** hex_nodes = new MeshLib::Node*[8];
-						for (unsigned k(0); k<8; k++)
-							hex_nodes[k] = nodes[node_ids[k]];
-						elem = new MeshLib::Hex(hex_nodes, mat_id);
-						break;
-					}
-					case 6: {
-						readNodeIDs(in, 6, node_ids, id_map);
-						MeshLib::Node** prism_nodes = new MeshLib::Node*[6];
-						for (unsigned k(0); k<6; k++)
-							prism_nodes[k] = nodes[node_ids[k]];
-						elem = new MeshLib::Prism(prism_nodes, mat_id);
-						break;
-					}
-					case 7: {
-						readNodeIDs(in, 5, node_ids, id_map);
-						MeshLib::Node** pyramid_nodes = new MeshLib::Node*[5];
-						for (unsigned k(0); k<5; k++)
-							pyramid_nodes[k] = nodes[node_ids[k]];
-						elem = new MeshLib::Pyramid(pyramid_nodes, mat_id);
-						break;
-					}
-					case 15:
-						in >> dummy; // skip rest of line
-						continue;
-						break;
-					default:
-						std::cout << "Error in GMSHInterface::readGMSHMesh() - Unknown element type " << type << "." << std::endl;
+				case 1: {
+					readNodeIDs(in, 2, node_ids, id_map);
+					// edge_nodes array will be deleted from Edge object
+					MeshLib::Node** edge_nodes = new MeshLib::Node*[2];
+					edge_nodes[0] = nodes[node_ids[0]];
+					edge_nodes[1] = nodes[node_ids[1]];
+					elem = new MeshLib::Edge(edge_nodes, 0);
+					break;
+				}
+				case 2: {
+					readNodeIDs(in, 3, node_ids, id_map);
+					MeshLib::Node** tri_nodes = new MeshLib::Node*[3];
+					tri_nodes[0] = nodes[node_ids[2]];
+					tri_nodes[1] = nodes[node_ids[1]];
+					tri_nodes[2] = nodes[node_ids[0]];
+					elem = new MeshLib::Tri(tri_nodes, mat_id);
+					break;
+				}
+				case 3: {
+					readNodeIDs(in, 4, node_ids, id_map);
+					MeshLib::Node** quad_nodes = new MeshLib::Node*[4];
+					for (unsigned k(0); k < 4; k++)
+						quad_nodes[k] = nodes[node_ids[k]];
+					elem = new MeshLib::Quad(quad_nodes, mat_id);
+					break;
+				}
+				case 4: {
+					readNodeIDs(in, 4, node_ids, id_map);
+					MeshLib::Node** tet_nodes = new MeshLib::Node*[5];
+					for (unsigned k(0); k < 4; k++)
+						tet_nodes[k] = nodes[node_ids[k]];
+					elem = new MeshLib::Tet(tet_nodes, mat_id);
+					break;
+				}
+				case 5: {
+					readNodeIDs(in, 8, node_ids, id_map);
+					MeshLib::Node** hex_nodes = new MeshLib::Node*[8];
+					for (unsigned k(0); k < 8; k++)
+						hex_nodes[k] = nodes[node_ids[k]];
+					elem = new MeshLib::Hex(hex_nodes, mat_id);
+					break;
+				}
+				case 6: {
+					readNodeIDs(in, 6, node_ids, id_map);
+					MeshLib::Node** prism_nodes = new MeshLib::Node*[6];
+					for (unsigned k(0); k < 6; k++)
+						prism_nodes[k] = nodes[node_ids[k]];
+					elem = new MeshLib::Prism(prism_nodes, mat_id);
+					break;
+				}
+				case 7: {
+					readNodeIDs(in, 5, node_ids, id_map);
+					MeshLib::Node** pyramid_nodes = new MeshLib::Node*[5];
+					for (unsigned k(0); k < 5; k++)
+						pyramid_nodes[k] = nodes[node_ids[k]];
+					elem = new MeshLib::Pyramid(pyramid_nodes, mat_id);
+					break;
+				}
+				case 15:
+					in >> dummy; // skip rest of line
+					continue;
+					break;
+				default:
+						WARN("GMSHInterface::readGMSHMesh(): Unknown element type %d.", type);
 				}
 				in >> std::ws;
 
-				if (type>0 && type<8)
+				if (type > 0 && type < 8)
 					elements.push_back(elem);
 			}
 
@@ -213,10 +223,13 @@ MeshLib::Mesh* GMSHInterface::readGMSHMesh(std::string const& fname)
 	return new MeshLib::Mesh(BaseLib::extractBaseNameWithoutExtension(fname), nodes, elements);
 }
 
-void GMSHInterface::readNodeIDs(std::ifstream &in, unsigned n_nodes, std::vector<unsigned> &node_ids, std::map<unsigned, unsigned> &id_map)
+void GMSHInterface::readNodeIDs(std::ifstream &in,
+                                unsigned n_nodes,
+                                std::vector<unsigned> &node_ids,
+                                std::map<unsigned, unsigned> &id_map)
 {
 	unsigned idx;
-	for (unsigned i=0; i<n_nodes; i++)
+	for (unsigned i = 0; i < n_nodes; i++)
 	{
 		in >> idx;
 		node_ids.push_back(id_map[idx]);
@@ -225,9 +238,9 @@ void GMSHInterface::readNodeIDs(std::ifstream &in, unsigned n_nodes, std::vector
 
 int GMSHInterface::write(std::ostream& out)
 {
-	out << "// GMSH input file created by OpenGeoSys " << OGS_VERSION_AND_PERSONS << " built on ";
+	out << "// GMSH input file created by OpenGeoSys " << OGS_VERSION_AND_PERSONS;
 #ifdef BUILD_TIMESTAMP
-	out << BUILD_TIMESTAMP;
+	out << " built on " << BUILD_TIMESTAMP;
 #endif
 	out << std::endl << std::endl;
 
@@ -237,15 +250,14 @@ int GMSHInterface::write(std::ostream& out)
 
 void GMSHInterface::writeGMSHInputFile(std::ostream& out)
 {
-#ifndef NDEBUG
-	std::cerr << "[GMSHInterface::writeGMSHInputFile] get data from GEOObjects ... " << std::flush;
-#endif
+	DBUG("GMSHInterface::writeGMSHInputFile(): get data from GEOObjects.");
+
 	// *** get and merge data from _geo_objs
 	_gmsh_geo_name = "GMSHGeometry";
 	_geo_objs.mergeGeometries(_selected_geometries, _gmsh_geo_name);
 	std::vector<GeoLib::Point*> * merged_pnts(const_cast<std::vector<GeoLib::Point*> *>(_geo_objs.getPointVec(_gmsh_geo_name)));
 	if (! merged_pnts) {
-		std::cerr << "[GMSHInterface::writeGMSHInputFile] did not found any points" << std::endl;
+		ERR("GMSHInterface::writeGMSHInputFile(): Did not found any points.");
 		return;
 	} else {
 		const size_t n_pnts(merged_pnts->size());
@@ -254,9 +266,7 @@ void GMSHInterface::writeGMSHInputFile(std::ostream& out)
 		}
 	}
 	std::vector<GeoLib::Polyline*> const* merged_plys(_geo_objs.getPolylineVec(_gmsh_geo_name));
-#ifndef NDEBUG
-	std::cerr << "ok" << std::endl;
-#endif
+	DBUG("GMSHInterface::writeGMSHInputFile(): \t ok.");
 
 	// *** compute topological hierarchy of polygons
 	if (merged_plys) {
@@ -266,11 +276,9 @@ void GMSHInterface::writeGMSHInputFile(std::ostream& out)
 				_polygon_tree_list.push_back(new GMSHPolygonTree(new GeoLib::Polygon(*(*it), true), NULL, _geo_objs, _gmsh_geo_name, _mesh_density_strategy));
 			}
 		}
-		std::cout << "[GMSHInterface::writeGMSHInputFile] compute topological hierarchy - detected "
-						<< _polygon_tree_list.size() << " polygons" << std::endl;
+		DBUG("GMSHInterface::writeGMSHInputFile(): Compute topological hierarchy - detected %d polygons.", _polygon_tree_list.size());
 		GeoLib::createPolygonTrees<FileIO::GMSHPolygonTree>(_polygon_tree_list);
-		std::cout << "[GMSHInterface::writeGMSHInputFile] compute topological hierarchy - calculated "
-								<< _polygon_tree_list.size() << " polygon trees" << std::endl;
+		DBUG("GMSHInterface::writeGMSHInputFile(): Compute topological hierarchy - calculated %d polygon trees.", _polygon_tree_list.size());
 	} else {
 		return;
 	}