diff --git a/BaseLib/BuildInfo.h.in b/BaseLib/BuildInfo.h.in
index 4a6b12c45b87f0b607cf801fefc8e08d8fae3e1f..0547210518db050b64c57de34387577c74e98663 100644
--- a/BaseLib/BuildInfo.h.in
+++ b/BaseLib/BuildInfo.h.in
@@ -36,4 +36,8 @@
 #define CMAKE_CXX_FLAGS ""
 #endif // CMAKE_CXX_FLAGS
 
+#ifndef CMAKE_CXX_FLAGS_DEBUG
+#define CMAKE_CXX_FLAGS_DEBUG ""
+#endif // CMAKE_CXX_FLAGS_DEBUG
+
 #endif // BUILDINFO_H
diff --git a/BaseLib/Configure.h.in b/BaseLib/Configure.h.in
index f042ad5e435700c2529c94fd60ce548dd580c879..854b9fc4223e84df8add91ea8738602a9a84a142 100644
--- a/BaseLib/Configure.h.in
+++ b/BaseLib/Configure.h.in
@@ -36,6 +36,7 @@
 #define SOURCEPATH "${CMAKE_SOURCE_DIR}"
 
 #cmakedefine OGS_VERSION "${OGS_VERSION}"
+#cmakedefine OGS_VERSION_AND_PERSONS "${OGS_VERSION_AND_PERSONS}"
 #cmakedefine OGS_DATE "${OGS_DATE}"
 
 // for tests
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 30b3b673bd1d7e08efce9051b1d7765b63d664d3..83d15a730449d4c53bfae09fa2bfa5b19f4e3a01 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -3,10 +3,10 @@
 #####################
 
 # Specify minimum CMake version
-cmake_minimum_required(VERSION 2.8.3)
+CMAKE_MINIMUM_REQUIRED(VERSION 2.8.3)
 
 # Project name
-project( OGS-6 )
+PROJECT( OGS-6 )
 
 ### CMake includes ###
 INCLUDE(scripts/cmake/CheckTypeSizes.cmake)
@@ -19,6 +19,10 @@ IF(NOT OGS_NO_EXTERNAL_LIBS)
 ENDIF() # NOT OGS_NO_EXTERNAL_LIBS
 INCLUDE(scripts/cmake/ProjectSetup.cmake)
 INCLUDE(scripts/cmake/DocumentationSetup.cmake)
+INCLUDE(scripts/cmake/Test.cmake)
+IF(OGS_COVERAGE)
+	INCLUDE(scripts/cmake/Coverage.cmake)
+ENDIF()
 
 ###########################################################################
 ### OGS version information. Adjust these if you release a new version. ###
@@ -27,7 +31,8 @@ SET (OGS_VERSION_MAJOR 0)
 SET (OGS_VERSION_MINOR 0)
 SET (OGS_VERSION_PATCH 1)
 SET (OGS_RELEASE_PERSONS "LB/TF/KR")
-SET (OGS_VERSION "${OGS_VERSION_MAJOR}.${OGS_VERSION_MINOR}.${OGS_VERSION_PATCH}(${OGS_RELEASE_PERSONS})")
+SET (OGS_VERSION "${OGS_VERSION_MAJOR}.${OGS_VERSION_MINOR}.${OGS_VERSION_PATCH}")
+SET (OGS_VERSION_AND_PERSONS "${OGS_VERSION} (${OGS_RELEASE_PERSONS})")
 SET (OGS_DATE "2012-08-20")
 
 ###############
@@ -66,6 +71,9 @@ IF (OGS_BUILD_INFO)
 	ADD_DEFINITIONS (-DOGS_BUILD_INFO)
 ENDIF (OGS_BUILD_INFO)
 
+# Code coverage
+OPTION(OGS_COVERAGE "Enables code coverage measurements with gcov/lcov." OFF)
+
 # Packaging
 OPTION(OGS_PACKAGING "Creating installers / packages" OFF)
 OPTION_REQUIRES(OGS_PACKAGING_ZIP "Do you want to package as zip?" OGS_PACKAGING)
diff --git a/FileIO/GMSInterface.cpp b/FileIO/GMSInterface.cpp
index 63c0ab4bb7ab306923a4974385220af2f40e1eaa..63e03d1fb3889ce6caae17f25398b0d39cfb6e78 100644
--- a/FileIO/GMSInterface.cpp
+++ b/FileIO/GMSInterface.cpp
@@ -280,7 +280,6 @@ MeshLib::Mesh* GMSInterface::readGMS3DMMesh(std::string filename)
 	unsigned node_idx[6], mat_id;
 	while ( getline(in, line) )
 	{
-		MeshLib::Element* elem (NULL);
 		std::string element_id(line.substr(0,3));
 		std::stringstream str(line);
 
@@ -288,25 +287,29 @@ MeshLib::Mesh* GMSInterface::readGMS3DMMesh(std::string filename)
 		{
 			str >> dummy >> id >> node_idx[0] >> node_idx[1] >> node_idx[2] >> node_idx[3]
 			    >> node_idx[4] >> node_idx[5] >> mat_id;
-			elem = new MeshLib::Prism(nodes[id_map.find(node_idx[0])->second], nodes[id_map.find(node_idx[1])->second],
-									  nodes[id_map.find(node_idx[2])->second], nodes[id_map.find(node_idx[3])->second],
-									  nodes[id_map.find(node_idx[4])->second], nodes[id_map.find(node_idx[5])->second], mat_id);
-			elements.push_back(elem);
+			MeshLib::Node** prism_nodes(new MeshLib::Node*[6]);
+			for (unsigned k(0); k<6; k++) {
+				prism_nodes[k] = nodes[id_map.find(node_idx[k])->second];
+			}
+			elements.push_back(new MeshLib::Prism(prism_nodes, mat_id));
 		}
 		else if (element_id.compare("E4T") == 0) // Tet
 		{
 			str >> dummy >> id >> node_idx[0] >> node_idx[1] >> node_idx[2] >> node_idx[3] >> mat_id;
-			elem = new MeshLib::Tet(nodes[id_map.find(node_idx[0])->second], nodes[id_map.find(node_idx[1])->second],
-				                    nodes[id_map.find(node_idx[2])->second], nodes[id_map.find(node_idx[3])->second], mat_id);
-			elements.push_back(elem);
+			MeshLib::Node** tet_nodes(new MeshLib::Node*[4]);
+			for (unsigned k(0); k<4; k++) {
+				tet_nodes[k] = nodes[id_map.find(node_idx[k])->second];
+			}
+			elements.push_back(new MeshLib::Tet(tet_nodes, mat_id));
 		}
 		else if ((element_id.compare("E4P") == 0) || (element_id.compare("E5P") == 0)) // Pyramid (both do exist for some reason)
 		{
 			str >> dummy >> id >> node_idx[0] >> node_idx[1] >> node_idx[2] >> node_idx[3] >> node_idx[4] >> mat_id;
-			elem = new MeshLib::Pyramid(nodes[id_map.find(node_idx[0])->second], nodes[id_map.find(node_idx[1])->second],
-				                        nodes[id_map.find(node_idx[2])->second], nodes[id_map.find(node_idx[3])->second],
-										nodes[id_map.find(node_idx[4])->second], mat_id);
-			elements.push_back(elem);
+			MeshLib::Node** pyramid_nodes(new MeshLib::Node*[5]);
+			for (unsigned k(0); k<5; k++) {
+				pyramid_nodes[k] = nodes[id_map.find(node_idx[k])->second];
+			}
+			elements.push_back(new MeshLib::Pyramid(pyramid_nodes, mat_id));
 		}
 		else if (element_id.compare("ND ") == 0) // Node
 		{
diff --git a/FileIO/Legacy/MeshIO.cpp b/FileIO/Legacy/MeshIO.cpp
index 31068bc462903879985cec25a96b073c8fc725c4..ecbaed45653b672df3735322b049c1d45074b124 100644
--- a/FileIO/Legacy/MeshIO.cpp
+++ b/FileIO/Legacy/MeshIO.cpp
@@ -133,41 +133,71 @@ MeshLib::Element* MeshIO::readElement(const std::string& line, const std::vector
 
 	switch(elem_type)
 	{
-	case MshElemType::EDGE:
+	case MshElemType::EDGE: {
 		for (int i = 0; i < 2; i++)
 			ss >> idx[i];
-		elem = new MeshLib::Edge(nodes[idx[1]], nodes[idx[0]], patch_index);
+		// edge_nodes array will be deleted from Edge object
+		MeshLib::Node** edge_nodes(new MeshLib::Node*[2]);
+		edge_nodes[0] = nodes[idx[1]];
+		edge_nodes[1] = nodes[idx[0]];
+		elem = new MeshLib::Edge(edge_nodes, patch_index);
 		break;
-	case MshElemType::TRIANGLE:
+	}
+	case MshElemType::TRIANGLE: {
 		for (int i = 0; i < 3; i++)
 			ss >> idx[i];
-		elem = new MeshLib::Tri(nodes[idx[2]], nodes[idx[1]], nodes[idx[0]], patch_index);
+		MeshLib::Node** tri_nodes(new MeshLib::Node*[3]);
+		tri_nodes[0] = nodes[idx[2]];
+		tri_nodes[1] = nodes[idx[1]];
+		tri_nodes[2] = nodes[idx[0]];
+		elem = new MeshLib::Tri(tri_nodes, patch_index);
 		break;
-	case MshElemType::QUAD:
+	}
+	case MshElemType::QUAD: {
 		for (int i = 0; i < 4; i++)
 			ss >> idx[i];
-		elem = new MeshLib::Quad(nodes[idx[3]], nodes[idx[2]], nodes[idx[1]], nodes[idx[0]], patch_index);
+		MeshLib::Node** quad_nodes(new MeshLib::Node*[4]);
+		for (unsigned k(0); k<4; k++)
+			quad_nodes[k] = nodes[idx[4-(k+1)]];
+		elem = new MeshLib::Quad(quad_nodes, patch_index);
 		break;
-	case MshElemType::TETRAHEDRON:
+	}
+	case MshElemType::TETRAHEDRON: {
 		for (int i = 0; i < 4; i++)
 			ss >> idx[i];
-		elem = new MeshLib::Tet(nodes[idx[3]], nodes[idx[2]], nodes[idx[1]], nodes[idx[0]], patch_index);
+		MeshLib::Node** tet_nodes(new MeshLib::Node*[4]);
+		for (unsigned k(0); k<4; k++)
+			tet_nodes[k] = nodes[idx[4-(k+1)]];
+		elem = new MeshLib::Tet(tet_nodes, patch_index);
 		break;
-	case MshElemType::HEXAHEDRON:
+	}
+	case MshElemType::HEXAHEDRON: {
 		for (int i = 0; i < 8; i++)
 			ss >> idx[i];
-		elem = new MeshLib::Hex(nodes[idx[7]], nodes[idx[6]], nodes[idx[5]], nodes[idx[4]], nodes[idx[3]], nodes[idx[2]], nodes[idx[1]], nodes[idx[0]], patch_index);
+		MeshLib::Node** hex_nodes(new MeshLib::Node*[8]);
+		for (unsigned k(0); k<8; k++)
+			hex_nodes[k] = nodes[idx[8-(k+1)]];
+		elem = new MeshLib::Hex(hex_nodes, patch_index);
 		break;
-	case MshElemType::PYRAMID:
+	}
+	case MshElemType::PYRAMID: {
 		for (int i = 0; i < 5; i++)
 			ss >> idx[i];
-		elem = new MeshLib::Pyramid(nodes[idx[4]], nodes[idx[3]], nodes[idx[2]], nodes[idx[1]], nodes[idx[0]], patch_index);
+		MeshLib::Node** pyramid_nodes(new MeshLib::Node*[5]);
+		for (unsigned k(0); k<5; k++)
+			pyramid_nodes[k] = nodes[idx[5-(k+1)]];
+		elem = new MeshLib::Pyramid(pyramid_nodes, patch_index);
 		break;
-	case MshElemType::PRISM:
+	}
+	case MshElemType::PRISM: {
 		for (int i = 0; i < 6; i++)
 			ss >> idx[i];
-		elem = new MeshLib::Prism(nodes[idx[5]], nodes[idx[4]], nodes[idx[3]], nodes[idx[2]], nodes[idx[1]], nodes[idx[0]], patch_index);
+		MeshLib::Node** prism_nodes(new MeshLib::Node*[6]);
+		for (unsigned k(0); k<6; k++)
+			prism_nodes[k] = nodes[idx[6-(k+1)]];
+		elem = new MeshLib::Prism(prism_nodes, patch_index);
 		break;
+	}
 	default:
 		elem = NULL;
 	}
diff --git a/FileIO/MeshIO/GMSHInterface.cpp b/FileIO/MeshIO/GMSHInterface.cpp
index 32ff2a5e8a3f7082893b2f73b95cefd8013be100..143004cec140cedec18643c9d278217f524e7e54 100644
--- a/FileIO/MeshIO/GMSHInterface.cpp
+++ b/FileIO/MeshIO/GMSHInterface.cpp
@@ -136,37 +136,64 @@ MeshLib::Mesh* GMSHInterface::readGMSHMesh(std::string const& fname)
 
 				switch (type)
 				{
-					case 1:
+					case 1: {
 						readNodeIDs(in, 2, node_ids, id_map);
-						elem = new MeshLib::Edge(nodes[node_ids[0]], nodes[node_ids[1]], 0);
+						// 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:
+					}
+					case 2: {
 						readNodeIDs(in, 3, node_ids, id_map);
-						elem = new MeshLib::Tri(nodes[node_ids[2]], nodes[node_ids[1]], nodes[node_ids[0]], mat_id);
+						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:
+					}
+					case 3: {
 						readNodeIDs(in, 4, node_ids, id_map);
-						elem = new MeshLib::Quad(nodes[node_ids[0]], nodes[node_ids[1]], nodes[node_ids[2]], nodes[node_ids[3]], mat_id);
+						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:
+					}
+					case 4: {
 						readNodeIDs(in, 4, node_ids, id_map);
-						elem = new MeshLib::Tet(nodes[node_ids[0]], nodes[node_ids[1]], nodes[node_ids[2]], nodes[node_ids[3]], mat_id);
+						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:
+					}
+					case 5: {
 						readNodeIDs(in, 8, node_ids, id_map);
-						elem = new MeshLib::Hex(nodes[node_ids[0]], nodes[node_ids[1]], nodes[node_ids[2]], nodes[node_ids[3]],
-							                    nodes[node_ids[4]], nodes[node_ids[5]], nodes[node_ids[6]], nodes[node_ids[7]], mat_id);
+						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:
+					}
+					case 6: {
 						readNodeIDs(in, 6, node_ids, id_map);
-						elem = new MeshLib::Prism(nodes[node_ids[0]], nodes[node_ids[1]], nodes[node_ids[2]],
-							                      nodes[node_ids[3]], nodes[node_ids[4]], nodes[node_ids[5]], mat_id);
+						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:
+					}
+					case 7: {
 						readNodeIDs(in, 5, node_ids, id_map);
-						elem = new MeshLib::Pyramid(nodes[node_ids[0]], nodes[node_ids[1]], nodes[node_ids[2]],
-							                        nodes[node_ids[3]], nodes[node_ids[4]], mat_id);
+						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;
@@ -199,7 +226,7 @@ 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 << " built on ";
+	out << "// GMSH input file created by OpenGeoSys " << OGS_VERSION_AND_PERSONS << " built on ";
 #ifdef BUILD_TIMESTAMP
 	out << BUILD_TIMESTAMP;
 #endif
diff --git a/FileIO/MeshIO/TetGenInterface.cpp b/FileIO/MeshIO/TetGenInterface.cpp
index 144f98e693bf5961544fb166581f0d241ecab79e..dd2b9b07bf77e5c7b408f495ee595b62c034a87c 100644
--- a/FileIO/MeshIO/TetGenInterface.cpp
+++ b/FileIO/MeshIO/TetGenInterface.cpp
@@ -326,7 +326,11 @@ bool TetGenInterface::parseElements(std::ifstream& ins, size_t n_tets, size_t n_
 					}
 				}
 				// insert new element into vector
-				_elements.push_back (new MeshLib::Tet(_nodes[ids[0]], _nodes[ids[1]], _nodes[ids[2]], _nodes[ids[3]], region));
+				MeshLib::Node** tet_nodes(new MeshLib::Node*[4]);
+				for (unsigned k(0); k<4; k++) {
+					tet_nodes[k] = _nodes[ids[k]];
+				}
+				_elements.push_back (new MeshLib::Tet(tet_nodes, region));
 
 			}
 		}
diff --git a/FileIO/XmlIO/VTKInterface.cpp b/FileIO/XmlIO/VTKInterface.cpp
index 5a02297e5b76fefa9b2205d4ec5c1c9ae53758b2..35720b60545012714f0a6a3a3961990b2fbfb1d5 100644
--- a/FileIO/XmlIO/VTKInterface.cpp
+++ b/FileIO/XmlIO/VTKInterface.cpp
@@ -32,7 +32,7 @@
 namespace FileIO {
 
 using namespace rapidxml;
-	
+
 VTKInterface::VTKInterface()
 : _export_name(""), _mesh(NULL), _doc(new xml_document<>), _use_compressor(false)
 {
@@ -96,7 +96,7 @@ MeshLib::Mesh* VTKInterface::readVTUFile(const std::string &file_name)
 			std::vector<unsigned> cell_types(nElems);
 
 			const rapidxml::xml_node<>* mat_id_node (piece_node->first_node("CellData")->first_node("DataArray"));
-			if (mat_id_node && 
+			if (mat_id_node &&
 				((std::string(mat_id_node->first_attribute("Name")->value()).compare("MaterialIDs") == 0) ||
 				 (std::string(mat_id_node->first_attribute("Name")->value()).compare("MatGroup") == 0)))
 			{
@@ -112,9 +112,9 @@ MeshLib::Mesh* VTKInterface::readVTUFile(const std::string &file_name)
 
 
 			const rapidxml::xml_node<>* points_node (piece_node->first_node("Points")->first_node("DataArray"));
-			// This _may_ have an attribute "Name" with the value "Points" but you cannot count on it. 
+			// This _may_ have an attribute "Name" with the value "Points" but you cannot count on it.
 			// However, there shouldn't be any other DataArray nodes so most likely not checking the name isn't a problem.
-			if (points_node) 
+			if (points_node)
 			{
 				if (std::string(points_node->first_attribute("format")->value()).compare("ascii") == 0)
 				{
@@ -185,46 +185,87 @@ MeshLib::Element* VTKInterface::readElement(std::stringstream &iss, const std::v
 	unsigned node_ids[8];
 	switch (type)
 	{
-	case 3: //line
+	case 3: { //line
 		for (unsigned i(0); i<2; i++) iss >> node_ids[i];
-		return new MeshLib::Edge(nodes[node_ids[0]], nodes[node_ids[1]], material);
+		MeshLib::Node** edge_nodes(new MeshLib::Node*[2]);
+		edge_nodes[0] = nodes[node_ids[0]];
+		edge_nodes[1] = nodes[node_ids[1]];
+		return new MeshLib::Edge(edge_nodes, material);
 		break;
-	case 5: //triangle
+	}
+	case 5: { //triangle
 		for (unsigned i(0); i<3; i++) iss >> node_ids[i];
-		return new MeshLib::Tri(nodes[node_ids[0]], nodes[node_ids[1]], nodes[node_ids[2]], material);
+		MeshLib::Node** tri_nodes(new MeshLib::Node*[3]);
+		tri_nodes[0] = nodes[node_ids[0]];
+		tri_nodes[1] = nodes[node_ids[1]];
+		tri_nodes[2] = nodes[node_ids[2]];
+		return new MeshLib::Tri(tri_nodes, material);
 		break;
-	case 9: //quad
+	}
+	case 9: { //quad
 		for (unsigned i(0); i<4; i++) iss >> node_ids[i];
-		return new MeshLib::Quad(nodes[node_ids[0]], nodes[node_ids[1]], nodes[node_ids[2]], nodes[node_ids[3]], material);
+		MeshLib::Node** quad_nodes(new MeshLib::Node*[4]);
+		for (unsigned k(0); k<4; k++)
+			quad_nodes[k] = nodes[node_ids[k]];
+		return new MeshLib::Quad(quad_nodes, material);
 		break;
-	case 8: //pixel
+	}
+	case 8: { //pixel
 		for (unsigned i(0); i<4; i++) iss >> node_ids[i];
-		return new MeshLib::Quad(nodes[node_ids[0]], nodes[node_ids[1]], nodes[node_ids[3]], nodes[node_ids[2]], material);
+		MeshLib::Node** quad_nodes(new MeshLib::Node*[4]);
+		quad_nodes[0] = nodes[node_ids[0]];
+		quad_nodes[1] = nodes[node_ids[1]];
+		quad_nodes[2] = nodes[node_ids[3]];
+		quad_nodes[3] = nodes[node_ids[2]];
+		return new MeshLib::Quad(quad_nodes, material);
 		break;
-	case 10:
+	}
+	case 10: {
 		for (unsigned i(0); i<4; i++) iss >> node_ids[i];
-		return new MeshLib::Tet(nodes[node_ids[0]], nodes[node_ids[1]], nodes[node_ids[2]], nodes[node_ids[3]], material);
+		MeshLib::Node** tet_nodes(new MeshLib::Node*[4]);
+		for (unsigned k(0); k<4; k++)
+			tet_nodes[k] = nodes[node_ids[k]];
+		return new MeshLib::Tet(tet_nodes, material);
 		break;
-	case 12: //hexahedron
+	}
+	case 12: { //hexahedron
 		for (unsigned i(0); i<8; i++) iss >> node_ids[i];
-		return new MeshLib::Hex(nodes[node_ids[0]], nodes[node_ids[1]], nodes[node_ids[2]], nodes[node_ids[3]],
-				                nodes[node_ids[4]], nodes[node_ids[5]], nodes[node_ids[6]], nodes[node_ids[7]], material);
+		MeshLib::Node** hex_nodes(new MeshLib::Node*[8]);
+		for (unsigned k(0); k<8; k++)
+			hex_nodes[k] = nodes[node_ids[k]];
+		return new MeshLib::Hex(hex_nodes, material);
 		break;
-	case 11: //voxel
+	}
+	case 11: { //voxel
 		for (unsigned i(0); i<8; i++) iss >> node_ids[i];
-		return new MeshLib::Hex(nodes[node_ids[0]], nodes[node_ids[1]], nodes[node_ids[3]], nodes[node_ids[2]],
-				                nodes[node_ids[4]], nodes[node_ids[5]], nodes[node_ids[7]], nodes[node_ids[6]], material);
+		MeshLib::Node** voxel_nodes(new MeshLib::Node*[8]);
+		voxel_nodes[0] = nodes[node_ids[0]];
+		voxel_nodes[1] = nodes[node_ids[1]];
+		voxel_nodes[2] = nodes[node_ids[3]];
+		voxel_nodes[3] = nodes[node_ids[2]];
+		voxel_nodes[4] = nodes[node_ids[4]];
+		voxel_nodes[5] = nodes[node_ids[5]];
+		voxel_nodes[6] = nodes[node_ids[7]];
+		voxel_nodes[7] = nodes[node_ids[6]];
+		return new MeshLib::Hex(voxel_nodes, material);
 		break;
-	case 14: //pyramid
+	}
+	case 14: { //pyramid
 		for (unsigned i(0); i<5; i++) iss >> node_ids[i];
-		return new MeshLib::Pyramid(nodes[node_ids[0]], nodes[node_ids[1]], nodes[node_ids[2]],
-				                    nodes[node_ids[3]], nodes[node_ids[4]], material);
+		MeshLib::Node** pyramid_nodes(new MeshLib::Node*[5]);
+		for (unsigned k(0); k<5; k++)
+			pyramid_nodes[k] = nodes[node_ids[k]];
+		return new MeshLib::Pyramid(pyramid_nodes, material);
 		break;
-	case 13: //wedge
+	}
+	case 13: { //wedge
 		for (unsigned i(0); i<6; i++) iss >> node_ids[i];
-		return new MeshLib::Prism(nodes[node_ids[0]], nodes[node_ids[1]], nodes[node_ids[2]],
-				                    nodes[node_ids[3]], nodes[node_ids[4]], nodes[node_ids[5]], material);
+		MeshLib::Node** prism_nodes(new MeshLib::Node*[6]);
+		for (unsigned k(0); k<6; k++)
+			prism_nodes[k] = nodes[node_ids[k]];
+		return new MeshLib::Prism(prism_nodes, material);
 		break;
+	}
 	default:
 		std::cout << "Error in VTKInterface::readElement() - Unknown mesh element type \"" << type << "\" ..." << std::endl;
 		return NULL;
diff --git a/Gui/DataView/CondFromRasterDialog.cpp b/Gui/DataView/CondFromRasterDialog.cpp
index 984fc350e8471542449dbab600c2dc5254c04f9a..7c8736b6aaed11dc7b4e6c26b7e6228c8836d877 100644
--- a/Gui/DataView/CondFromRasterDialog.cpp
+++ b/Gui/DataView/CondFromRasterDialog.cpp
@@ -42,7 +42,7 @@ CondFromRasterDialog::~CondFromRasterDialog()
 
 void CondFromRasterDialog::on_selectButton_pressed()
 {
-	QSettings settings("UFZ", "OpenGeoSys-5");
+	QSettings settings;
 #ifdef libgeotiff_FOUND
 	QString geotiffExtension(" *.tif");
 #else
diff --git a/Gui/DataView/ConditionWriterDialog.cpp b/Gui/DataView/ConditionWriterDialog.cpp
index 323188840b76b67264b6a43e6055c693f8f13283..5821a79c89e1b80ef7b7df8b6b17b8ac8204a2d2 100644
--- a/Gui/DataView/ConditionWriterDialog.cpp
+++ b/Gui/DataView/ConditionWriterDialog.cpp
@@ -44,7 +44,7 @@ void ConditionWriterDialog::on_fileNameButton_pressed()
 	else if ((geo_idx != 0) && (cnd_idx == 2)) filetypes = "OpenGeoSys FEM Condition file (*.cnd);;GeoSys Initial Condition (*.ic)";
 	else if ((geo_idx != 0) && (cnd_idx == 3)) filetypes = "OpenGeoSys FEM Condition file (*.cnd);;GeoSys Source Condition (*.st)";
 
-	QSettings settings("UFZ", "OpenGeoSys-5");
+	QSettings settings;
 	QString fileName = QFileDialog::getSaveFileName(this, "Select path",
 					settings.value("lastOpenedConditionsFileDirectory").toString(), filetypes);
 
diff --git a/Gui/DataView/DataView.cpp b/Gui/DataView/DataView.cpp
index f0a9b1a6a2564336be88e39e778bb781a090a592..60342a3a436731a5a4ac40536d893b9cb576751d 100644
--- a/Gui/DataView/DataView.cpp
+++ b/Gui/DataView/DataView.cpp
@@ -51,7 +51,7 @@ void DataView::updateView()
 
 void DataView::addMeshAction()
 {
-	QSettings settings("UFZ", "OpenGeoSys-5");
+	QSettings settings;
 	QString fileName =
 	        QFileDialog::getOpenFileName(this, "Select mesh file", settings.value(
 	                                             "lastOpenedFileDirectory").toString(),
@@ -128,7 +128,7 @@ int DataView::writeMeshToFile() const
 
 	if (mesh)
 	{
-		QSettings settings("UFZ", "OpenGeoSys-5");
+		QSettings settings;
 		QString mshName = QString::fromStdString(
 		        static_cast<MshModel*>(this->model())->getMesh(index)->getName());
 		QString fileName = QFileDialog::getSaveFileName(NULL, "Save mesh as",
diff --git a/Gui/DataView/DiagramView/DetailWindow.cpp b/Gui/DataView/DiagramView/DetailWindow.cpp
index ea5cdcfe85ff98fa9bfa2386f79f63b8bc3ee08d..f361d474155a2f3519201b9ac67f95e72bd047c1 100644
--- a/Gui/DataView/DiagramView/DetailWindow.cpp
+++ b/Gui/DataView/DiagramView/DetailWindow.cpp
@@ -136,7 +136,7 @@ void DetailWindow::addList(DiagramList* list, QColor c)
 
 void DetailWindow::on_addDataButton_clicked()
 {
-	QSettings settings("UFZ", "OpenGeoSys-5");
+	QSettings settings;
 	QString fileName = QFileDialog::getOpenFileName( this, "Select data file to open",
 	                                                 settings.value(
 	                                                         "lastOpenedFileDirectory").
diff --git a/Gui/DataView/MshEditDialog.cpp b/Gui/DataView/MshEditDialog.cpp
index e5019f191ff8e60a9ab7861e2a1b2866b7b86d5f..e9836806603d3d4a69b613bb5fbe04736d2a7a26 100644
--- a/Gui/DataView/MshEditDialog.cpp
+++ b/Gui/DataView/MshEditDialog.cpp
@@ -26,8 +26,8 @@
 #include <QVBoxLayout>
 
 MshEditDialog::MshEditDialog(const MeshLib::Mesh* mesh, QDialog* parent)
-	: QDialog(parent), _msh(mesh), _noDataDeleteBox(NULL), 
-	  _nLayerLabel (new QLabel("Please specify the number of layers to add:")),  
+	: QDialog(parent), _msh(mesh), _noDataDeleteBox(NULL),
+	  _nLayerLabel (new QLabel("Please specify the number of layers to add:")),
 	  _nLayerExplanation (new QLabel("(select \"0\" for surface mapping)")),
 	  _layerEdit (new QLineEdit("0")),
 	  _nextButton (new QPushButton("Next")),
@@ -58,7 +58,7 @@ MshEditDialog::~MshEditDialog()
 	}
 	for (int i = 0; i < _buttons.size(); ++i)
 		delete _buttons[i];
-	
+
 	delete _nLayerLabel;
 	delete _nLayerExplanation;
 	delete _layerEdit;
@@ -94,13 +94,13 @@ void MshEditDialog::nextButtonPressed()
 		_radioButtonBox->setLayout(_radiobuttonLayout);
 		gridLayoutLayerMapping->addWidget(_radioButtonBox, 2, 0, 1, 3);
 		// add an empty line to better arrange the following information
-		gridLayoutLayerMapping->addWidget(_nLayerExplanation, 3, 0); 
+		gridLayoutLayerMapping->addWidget(_nLayerExplanation, 3, 0);
 		connect(_selectButton1, SIGNAL(pressed()), this, SLOT(createWithRasters()));
 		connect(_selectButton2, SIGNAL(pressed()), this, SLOT(createStatic()));
 	}
 	else
 		this->createWithRasters();
-			
+
 
 }
 
@@ -133,7 +133,7 @@ void MshEditDialog::createWithRasters()
 		this->_layerSelectionLayout->addWidget(_buttons[i], i, 2);
 
 		// don't add bottom layer if mesh contains only surface
-		if (this->_n_layers==0) break; 
+		if (this->_n_layers==0) break;
 	}
 	this->_layerBox->setLayout(this->_layerSelectionLayout);
 	this->_noDataDeleteBox = new QCheckBox("Remove mesh nodes at NoData values");
@@ -224,7 +224,7 @@ void MshEditDialog::accept()
 
 			if (new_mesh)
 				emit mshEditFinished(new_mesh);
-			
+
 			if (!new_mesh || result==0)
 				OGSError::box("Error creating mesh");
 
@@ -245,7 +245,7 @@ void MshEditDialog::reject()
 void MshEditDialog::getFileName()
 {
 	QPushButton* button = dynamic_cast<QPushButton*>(this->sender());
-	QSettings settings("UFZ", "OpenGeoSys-5");
+	QSettings settings;
 	QString filename = QFileDialog::getOpenFileName(this, "Select raster file to open",
 	                                                settings.value("lastOpenedFileDirectory").toString(),
 	                                                "ASCII raster files (*.asc);;All files (* *.*)");
diff --git a/Gui/Vrpn/QVrpnArtTrackingClient.cpp b/Gui/Vrpn/QVrpnArtTrackingClient.cpp
index 7755261a401a4a5c2447ee3487c67c09247df269..062c78cabfaeadf8b3332177b30c82885127b25e 100644
--- a/Gui/Vrpn/QVrpnArtTrackingClient.cpp
+++ b/Gui/Vrpn/QVrpnArtTrackingClient.cpp
@@ -32,7 +32,7 @@ QVrpnArtTrackingClient::~QVrpnArtTrackingClient()
 {
 	QStringList list = _deviceName.split("@");
 
-	QSettings settings("UFZ", "OpenGeoSys-5");
+	QSettings settings;
 	settings.beginGroup("Tracking");
 	settings.setValue("artDeviceName", list.at(0));
 	settings.setValue("artDeviceNameAt", list.at(1));
diff --git a/Gui/VtkVis/VisualizationWidget.cpp b/Gui/VtkVis/VisualizationWidget.cpp
index fa4add6c7fad61ca7969aff093944ca37e13dab9..939a3a1229dc2c690684b2d2a51bd6e2928ddcd0 100644
--- a/Gui/VtkVis/VisualizationWidget.cpp
+++ b/Gui/VtkVis/VisualizationWidget.cpp
@@ -72,7 +72,7 @@ VisualizationWidget::VisualizationWidget( QWidget* parent /*= 0*/ )
 	_interactorStyle->SetDefaultRenderer(_vtkRender);
 #endif // OGS_VRED_PLUGIN
 
-	QSettings settings("UFZ", "OpenGeoSys-5");
+	QSettings settings;
 
 #ifdef OGS_USE_VRPN
 	VtkTrackedCamera* cam = new VtkTrackedCamera(this);
@@ -145,7 +145,7 @@ VisualizationWidget::VisualizationWidget( QWidget* parent /*= 0*/ )
 VisualizationWidget::~VisualizationWidget()
 {
 	// Write settings
-	QSettings settings("UFZ", "OpenGeoSys-5");
+	QSettings settings;
 	settings.setValue("stereoEnabled", stereoToolButton->isChecked());
 	settings.setValue("stereoEyeAngle", _vtkRender->GetActiveCamera()->GetEyeAngle());
 
@@ -249,7 +249,7 @@ void VisualizationWidget::on_orthogonalProjectionToolButton_toggled( bool checke
 
 void VisualizationWidget::on_screenshotPushButton_pressed()
 {
-	QSettings settings("UFZ", "OpenGeoSys-5");
+	QSettings settings;
 	QString filename = QFileDialog::getSaveFileName(this, tr("Save screenshot"),
 	                                                settings.value(
 	                                                        "lastScreenshotDir").toString(),
diff --git a/Gui/VtkVis/VtkCompositeColormapToImageFilter.cpp b/Gui/VtkVis/VtkCompositeColormapToImageFilter.cpp
index b0daf3b630ffcb269fdc8000a3f7251c71956906..5716dd702138e631ce950060789a1dd8b07d673e 100644
--- a/Gui/VtkVis/VtkCompositeColormapToImageFilter.cpp
+++ b/Gui/VtkVis/VtkCompositeColormapToImageFilter.cpp
@@ -42,7 +42,7 @@ void VtkCompositeColormapToImageFilter::init()
 	vtkSmartPointer<VtkColorLookupTable> colormap;
 
 	QWidget* parent = 0;
-	QSettings settings("UFZ", "OpenGeoSys-5");
+	QSettings settings;
 	QString fileName = QFileDialog::getOpenFileName(parent,
 	                                                "Select color lookup table",
 	                                                settings.value("lastOpenedLookupTableFileDirectory").
diff --git a/Gui/VtkVis/VtkCompositeTextureOnSurfaceFilter.cpp b/Gui/VtkVis/VtkCompositeTextureOnSurfaceFilter.cpp
index 427b85017f84d5b50a97a17dbaf6ffba30eab1a3..dd34ba0140ddb9bd0bf6cacd99b0928aa8c442eb 100644
--- a/Gui/VtkVis/VtkCompositeTextureOnSurfaceFilter.cpp
+++ b/Gui/VtkVis/VtkCompositeTextureOnSurfaceFilter.cpp
@@ -55,7 +55,7 @@ void VtkCompositeTextureOnSurfaceFilter::init()
 		surface->SetInputConnection(_inputAlgorithm->GetOutputPort());
 
 	QWidget* parent = 0;
-	QSettings settings("UFZ", "OpenGeoSys-5");
+	QSettings settings;
 	QString fileName = QFileDialog::getOpenFileName(parent,
 	                                                "Select raster file to apply as texture",
 	                                                settings.value("lastOpenedTextureFileDirectory").
diff --git a/Gui/VtkVis/VtkMeshConverter.cpp b/Gui/VtkVis/VtkMeshConverter.cpp
index 169c0c4353e2d3360064878fc6bb90c5a7e92090..db3b59119f6d27c3d4293cdc67d13df4262d0371 100644
--- a/Gui/VtkVis/VtkMeshConverter.cpp
+++ b/Gui/VtkVis/VtkMeshConverter.cpp
@@ -207,18 +207,27 @@ MeshLib::Mesh* VtkMeshConverter::constructMesh(const double* pixVal,
 				const int mat = (intensity_type != UseIntensityAs::MATERIAL) ? 0 : static_cast<int>(pixVal[index+incHeight]);
 				if (elem_type == MshElemType::TRIANGLE)
 				{
-					MeshLib::Tri* tri1 (new MeshLib::Tri(nodes[node_idx_map[index]], nodes[node_idx_map[index + 1]],
-						                                 nodes[node_idx_map[index + incHeight]], mat));	// upper left triangle
-					MeshLib::Tri* tri2 (new MeshLib::Tri(nodes[node_idx_map[index + 1]], nodes[node_idx_map[index + incHeight + 1]],
-						                                 nodes[node_idx_map[index + incHeight]], mat));	// lower right triangle
-					elements.push_back(tri1);
-					elements.push_back(tri2);
+					MeshLib::Node** tri1_nodes (new MeshLib::Node*[3]);
+					tri1_nodes[0] = nodes[node_idx_map[index]];
+					tri1_nodes[1] = nodes[node_idx_map[index+1]];
+					tri1_nodes[2] = nodes[node_idx_map[index+incHeight]];
+
+					MeshLib::Node** tri2_nodes (new MeshLib::Node*[3]);
+					tri2_nodes[0] = nodes[node_idx_map[index+1]];
+					tri2_nodes[1] = nodes[node_idx_map[index+incHeight+1]];
+					tri2_nodes[2] = nodes[node_idx_map[index+incHeight]];
+
+					elements.push_back(new MeshLib::Tri(tri1_nodes, mat)); // upper left triangle
+					elements.push_back(new MeshLib::Tri(tri2_nodes, mat)); // lower right triangle
 				}
 				if (elem_type == MshElemType::QUAD)
 				{
-					MeshLib::Quad* quad (new MeshLib::Quad(nodes[node_idx_map[index]], nodes[node_idx_map[index + 1]],
-													 nodes[node_idx_map[index + incHeight + 1]], nodes[node_idx_map[index + incHeight]], mat));
-					elements.push_back(quad);
+					MeshLib::Node** quad_nodes (new MeshLib::Node*[4]);
+					quad_nodes[0] = nodes[node_idx_map[index]];
+					quad_nodes[1] = nodes[node_idx_map[index + 1]];
+					quad_nodes[2] = nodes[node_idx_map[index + incHeight + 1]];
+					quad_nodes[3] = nodes[node_idx_map[index + incHeight]];
+					elements.push_back(new MeshLib::Quad(quad_nodes, mat));
 				}
 			}
 		}
@@ -257,34 +266,70 @@ MeshLib::Mesh* VtkMeshConverter::convertUnstructuredGrid(vtkUnstructuredGrid* gr
 		int cell_type = grid->GetCellType(i);
 		switch (cell_type)
 		{
-		case VTK_TRIANGLE:
-			elem = new MeshLib::Tri(nodes[node_ids[0]], nodes[node_ids[1]], nodes[node_ids[2]], material);
+		case VTK_TRIANGLE: {
+			MeshLib::Node** tri_nodes(new MeshLib::Node*[3]);
+			for (unsigned k(0); k<3; k++)
+				tri_nodes[k] = nodes[node_ids[k]];
+			elem = new MeshLib::Tri(tri_nodes, material);
 			break;
-		case VTK_QUAD:
-			elem = new MeshLib::Quad(nodes[node_ids[0]], nodes[node_ids[1]], nodes[node_ids[2]], nodes[node_ids[3]], material);
+		}
+		case VTK_QUAD: {
+			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, material);
 			break;
-		case VTK_PIXEL:
-			elem = new MeshLib::Quad(nodes[node_ids[0]], nodes[node_ids[1]], nodes[node_ids[3]], nodes[node_ids[2]], material);
+		}
+		case VTK_PIXEL: {
+			MeshLib::Node** quad_nodes(new MeshLib::Node*[4]);
+			quad_nodes[0] = nodes[node_ids[0]];
+			quad_nodes[1] = nodes[node_ids[1]];
+			quad_nodes[2] = nodes[node_ids[3]];
+			quad_nodes[3] = nodes[node_ids[2]];
+			elem = new MeshLib::Quad(quad_nodes, material);
 			break;
-		case VTK_TETRA:
-			elem = new MeshLib::Tet(nodes[node_ids[0]], nodes[node_ids[1]], nodes[node_ids[2]], nodes[node_ids[3]], material);
+		}
+		case VTK_TETRA: {
+			MeshLib::Node** tet_nodes(new MeshLib::Node*[4]);
+			for (unsigned k(0); k<4; k++)
+				tet_nodes[k] = nodes[node_ids[k]];
+			elem = new MeshLib::Tet(tet_nodes, material);
 			break;
-		case VTK_HEXAHEDRON:
-			elem = new MeshLib::Hex(nodes[node_ids[0]], nodes[node_ids[1]], nodes[node_ids[2]], nodes[node_ids[3]],
-				                    nodes[node_ids[4]], nodes[node_ids[5]], nodes[node_ids[6]], nodes[node_ids[7]], material);
+		}
+		case VTK_HEXAHEDRON: {
+			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, material);
 			break;
-		case VTK_VOXEL:
-			elem = new MeshLib::Hex(nodes[node_ids[0]], nodes[node_ids[1]], nodes[node_ids[3]], nodes[node_ids[2]],
-				                    nodes[node_ids[4]], nodes[node_ids[5]], nodes[node_ids[7]], nodes[node_ids[6]], material);
+		}
+		case VTK_VOXEL: {
+			MeshLib::Node** voxel_nodes(new MeshLib::Node*[8]);
+			voxel_nodes[0] = nodes[node_ids[0]];
+			voxel_nodes[1] = nodes[node_ids[1]];
+			voxel_nodes[2] = nodes[node_ids[3]];
+			voxel_nodes[3] = nodes[node_ids[2]];
+			voxel_nodes[4] = nodes[node_ids[4]];
+			voxel_nodes[5] = nodes[node_ids[5]];
+			voxel_nodes[6] = nodes[node_ids[7]];
+			voxel_nodes[7] = nodes[node_ids[6]];
+			elem = new MeshLib::Hex(voxel_nodes, material);
 			break;
-		case VTK_PYRAMID:
-			elem = new MeshLib::Pyramid(nodes[node_ids[0]], nodes[node_ids[1]], nodes[node_ids[2]],
-				                        nodes[node_ids[3]], nodes[node_ids[4]], material);
+		}
+		case VTK_PYRAMID: {
+			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, material);
 			break;
-		case VTK_WEDGE:
-			elem = new MeshLib::Prism(nodes[node_ids[0]], nodes[node_ids[1]], nodes[node_ids[2]],
-				                      nodes[node_ids[3]], nodes[node_ids[4]], nodes[node_ids[5]], material);
+		}
+		case VTK_WEDGE: {
+			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, material);
 			break;
+		}
 		default:
 			std::cout << "Error in GridAdapter::convertUnstructuredGrid() - Unknown mesh element type \"" << cell_type << "\" ..." << std::endl;
 			return NULL;
diff --git a/Gui/VtkVis/VtkTrackedCamera.cpp b/Gui/VtkVis/VtkTrackedCamera.cpp
index 04c518ba36964b20c9fb0e6bcf5052b2016c1a88..1fbb0df2d7ef7e441719a362eb91fa91d818c1a5 100644
--- a/Gui/VtkVis/VtkTrackedCamera.cpp
+++ b/Gui/VtkVis/VtkTrackedCamera.cpp
@@ -20,7 +20,7 @@
 VtkTrackedCamera::VtkTrackedCamera(QObject* parent)
 	: QObject(parent), vtkOpenGLCamera()
 {
-	QSettings settings("UFZ", "OpenGeoSys-5");
+	QSettings settings;
 	settings.beginGroup("Tracking");
 	// Read from settings
 	if (settings.contains("artOffsetX"))
@@ -58,7 +58,7 @@ VtkTrackedCamera::VtkTrackedCamera(QObject* parent)
 
 VtkTrackedCamera::~VtkTrackedCamera()
 {
-	QSettings settings("UFZ", "OpenGeoSys-5");
+	QSettings settings;
 	settings.beginGroup("Tracking");
 	settings.setValue("artOffsetX", _trackedPositionOffset[0]);
 	settings.setValue("artOffsetY", _trackedPositionOffset[1]);
diff --git a/Gui/VtkVis/VtkVisPipeline.cpp b/Gui/VtkVis/VtkVisPipeline.cpp
index 5efc349f765991900e268c90fae449ffdf410adf..c3b0819353f61fc162d12b8138cb1a2a5acac3e9 100644
--- a/Gui/VtkVis/VtkVisPipeline.cpp
+++ b/Gui/VtkVis/VtkVisPipeline.cpp
@@ -71,7 +71,7 @@ VtkVisPipeline::VtkVisPipeline( vtkRenderer* renderer, QObject* parent /*= 0*/ )
 	delete _rootItem;
 	_rootItem = new TreeItem(rootData, NULL);
 
-	QSettings settings("UFZ", "OpenGeoSys-5");
+	QSettings settings;
 	QVariant backgroundColorVariant = settings.value("VtkBackgroundColor");
 	if (backgroundColorVariant != QVariant())
 		this->setBGColor(backgroundColorVariant.value<QColor>());
@@ -141,7 +141,7 @@ const QColor VtkVisPipeline::getBGColor() const
 
 void VtkVisPipeline::setBGColor(const QColor &color)
 {
-	QSettings settings("UFZ", "OpenGeoSys-5");
+	QSettings settings;
 	settings.setValue("VtkBackgroundColor", color);
 	_renderer->SetBackground(color.redF(), color.greenF(), color.blueF());
 }
diff --git a/Gui/VtkVis/VtkVisPipelineView.cpp b/Gui/VtkVis/VtkVisPipelineView.cpp
index aa7011617bcff4a3ce0c4dbbf727e70b7965669c..f8eafc0b6af9a9dd98ff02bedfbeb26165d3fead 100644
--- a/Gui/VtkVis/VtkVisPipelineView.cpp
+++ b/Gui/VtkVis/VtkVisPipelineView.cpp
@@ -124,7 +124,7 @@ void VtkVisPipelineView::contextMenuEvent( QContextMenuEvent* event )
 
 void VtkVisPipelineView::exportSelectedPipelineItemAsVtk()
 {
-	QSettings settings("UFZ", "OpenGeoSys-5");
+	QSettings settings;
 	QModelIndex idx = this->selectionModel()->currentIndex();
 	QString filename = QFileDialog::getSaveFileName(this, "Export object to vtk-file",
 	                                settings.value("lastExportedFileDirectory").toString(),
@@ -140,7 +140,7 @@ void VtkVisPipelineView::exportSelectedPipelineItemAsVtk()
 
 void VtkVisPipelineView::exportSelectedPipelineItemAsOsg()
 {
-	QSettings settings("UFZ", "OpenGeoSys-5");
+	QSettings settings;
 	QModelIndex idx = this->selectionModel()->currentIndex();
 	QString filename = QFileDialog::getSaveFileName(this, "Export object to OpenSG file",
 	                                                settings.value("lastExportedFileDirectory").
@@ -263,7 +263,7 @@ void VtkVisPipelineView::addColorTable()
 	                                                 this->selectionModel()->currentIndex())) );
 	const QString array_name = item->GetActiveAttribute();
 
-	QSettings settings("UFZ", "OpenGeoSys-5");
+	QSettings settings;
 	QString filename = QFileDialog::getOpenFileName(this, "Select color table",
 	                                                settings.value("lastOpenedLutFileDirectory"). toString(),
 													"Color table files (*.xml);;");
diff --git a/Gui/main.cpp b/Gui/main.cpp
index d7934a0e3482b09dc763a17f1046c0cff535446b..c810a3755a175242fcd3b21b4c659f4f4fb49124 100644
--- a/Gui/main.cpp
+++ b/Gui/main.cpp
@@ -17,9 +17,13 @@ int main(int argc, char* argv[])
 	BaseLib::LogogSimpleFormatter* formatter = new BaseLib::LogogSimpleFormatter;
 	logogCout->SetFormatter(*formatter);
 	QApplication a(argc, argv);
+	QApplication::setApplicationName("OpenGeoSys - Data Explorer");
+	QApplication::setApplicationVersion(QString(OGS_VERSION));
+	QApplication::setOrganizationName("OpenGeoSys Community");
+	QApplication::setOrganizationDomain("opengeosys.org");
 	setlocale(LC_NUMERIC,"C");
 	MainWindow* w = new MainWindow();
-	w->setWindowTitle( w->windowTitle() + " - " + QString(OGS_VERSION) + " - FirstFloor");
+	w->setWindowTitle( w->windowTitle() + " - " + QString(OGS_VERSION_AND_PERSONS) + " - FirstFloor");
 	w->show();
 	int returncode = a.exec();
 	delete w;
diff --git a/Gui/mainwindow.cpp b/Gui/mainwindow.cpp
index 438449cedb8bf729314530ae4123f3a812ac1928..4386667c0f6c51fb977ae13e4226c320e4489da7 100644
--- a/Gui/mainwindow.cpp
+++ b/Gui/mainwindow.cpp
@@ -395,7 +395,7 @@ void MainWindow::showVisDockWidget(bool show)
 
 void MainWindow::open(int file_type)
 {
-	QSettings settings("UFZ", "OpenGeoSys-5");
+	QSettings settings;
 	ImportFileType::type t = static_cast<ImportFileType::type>(file_type);
 	QString type_str = QString::fromStdString((ImportFileType::convertImportFileTypeToString(t)));
 	QString fileName = QFileDialog::getOpenFileName(this,
@@ -642,7 +642,7 @@ void MainWindow::loadFile(ImportFileType::type t, const QString &fileName)
 #endif
 	else if (t == ImportFileType::TETGEN)
 	{
-		QSettings settings("UFZ", "OpenGeoSys-5");
+		QSettings settings;
 		QString element_fname = QFileDialog::getOpenFileName(this, "Select TetGen element file",
 						                                     settings.value("lastOpenedTetgenFileDirectory").toString(),
 						                                     "TetGen element files (*.ele);;");
@@ -682,7 +682,7 @@ void MainWindow::updateDataViews()
 
 void MainWindow::readSettings()
 {
-	QSettings settings("UFZ", "OpenGeoSys-5");
+	QSettings settings;
 
 	restoreGeometry(settings.value("windowGeometry").toByteArray());
 	restoreState(settings.value("windowState").toByteArray());
@@ -690,7 +690,7 @@ void MainWindow::readSettings()
 
 void MainWindow::writeSettings()
 {
-	QSettings settings("UFZ", "OpenGeoSys-5");
+	QSettings settings;
 
 	settings.setValue("windowGeometry", saveGeometry());
 	settings.setValue("windowState", saveState());
@@ -698,7 +698,7 @@ void MainWindow::writeSettings()
 
 void MainWindow::about()
 {
-	QString ogsVersion = QString(OGS_VERSION);
+	QString ogsVersion = QString(OGS_VERSION_AND_PERSONS);
 
 	QString about = tr("Built on %1\nOGS Version: %2\n\n").arg(
 		QDate::currentDate().toString(Qt::ISODate)).arg(ogsVersion);
@@ -765,7 +765,7 @@ QMenu* MainWindow::createImportFilesMenu()
 
 void MainWindow::loadPetrelFiles()
 {
-	QSettings settings("UFZ", "OpenGeoSys-5");
+	QSettings settings;
 	QStringList sfc_file_names = QFileDialog::getOpenFileNames(
 	        this, "Select surface data file(s) to import", "", "Petrel files (*)");
 	QStringList well_path_file_names = QFileDialog::getOpenFileNames(
@@ -817,7 +817,7 @@ void MainWindow::showAddPipelineFilterItemDialog(QModelIndex parentIndex)
 
 void MainWindow::loadFEMConditions(std::string geoName)
 {
-	QSettings settings("UFZ", "OpenGeoSys-5");
+	QSettings settings;
 	QString fileName = QFileDialog::getOpenFileName( this, "Select data file to open",
 														settings.value("lastOpenedFileDirectory").toString(),
 														"Geosys FEM condition files (*.cnd *.bc *.ic *.st);;All files (* *.*)");
@@ -1031,7 +1031,7 @@ void MainWindow::showDiagramPrefsDialog(QModelIndex &index)
 
 void MainWindow::showDiagramPrefsDialog()
 {
-	QSettings settings("UFZ", "OpenGeoSys-5");
+	QSettings settings;
 	QString fileName = QFileDialog::getOpenFileName( this, "Select data file to open",
 	                                                 settings.value(
 	                                                         "lastOpenedFileDirectory").
@@ -1192,7 +1192,7 @@ void MainWindow::HideWindow()
 void MainWindow::on_actionExportVTK_triggered(bool checked /*= false*/)
 {
 	Q_UNUSED(checked)
-	QSettings settings("UFZ", "OpenGeoSys-5");
+	QSettings settings;
 	int count = 0;
 	QString filename = QFileDialog::getSaveFileName(this,
 	                                                "Export object to vtk-files",
@@ -1222,7 +1222,7 @@ void MainWindow::on_actionExportVTK_triggered(bool checked /*= false*/)
 void MainWindow::on_actionExportVRML2_triggered(bool checked /*= false*/)
 {
 	Q_UNUSED(checked)
-	QSettings settings("UFZ", "OpenGeoSys-5");
+	QSettings settings;
 	QString fileName = QFileDialog::getSaveFileName(this,
 	                                                "Save scene to VRML file", settings.value(
 	                                                        "lastExportedFileDirectory").
@@ -1245,7 +1245,7 @@ void MainWindow::on_actionExportVRML2_triggered(bool checked /*= false*/)
 void MainWindow::on_actionExportObj_triggered(bool checked /*= false*/)
 {
 	Q_UNUSED(checked)
-	QSettings settings("UFZ", "OpenGeoSys-5");
+	QSettings settings;
 	QString fileName = QFileDialog::getSaveFileName(this,
 	                                                "Save scene to Wavefront OBJ files",
 	                                                settings.value(
@@ -1270,7 +1270,7 @@ void MainWindow::on_actionExportOpenSG_triggered(bool checked /*= false*/)
 {
 	Q_UNUSED(checked)
 #ifdef OGS_USE_OPENSG
-	QSettings settings("UFZ", "OpenGeoSys-5");
+	QSettings settings;
 	QString filename = QFileDialog::getSaveFileName(
 	        this, "Export scene to OpenSG binary file", settings.value(
 	                "lastExportedFileDirectory").toString(), "OpenSG files (*.osb);;");
@@ -1382,7 +1382,7 @@ void MainWindow::quitPresentationMode()
 
 QString MainWindow::getLastUsedDir()
 {
-	QSettings settings("UFZ", "OpenGeoSys-5");
+	QSettings settings;
 	QString fileName("");
 	QStringList files = settings.value("recentFileList").toStringList();
 	if (files.size() != 0)
diff --git a/MeshLib/Elements/Edge.cpp b/MeshLib/Elements/Edge.cpp
deleted file mode 100644
index 2722b9a13def6939e8fcc2d48f97c7df5e5531d3..0000000000000000000000000000000000000000
--- a/MeshLib/Elements/Edge.cpp
+++ /dev/null
@@ -1,77 +0,0 @@
-/**
- * Copyright (c) 2012, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- *
- * \file Edge.cpp
- *
- * Created on 2012-05-02 by Karsten Rink
- */
-
-#include "Edge.h"
-#include "Node.h"
-
-#include "MathTools.h"
-
-namespace MeshLib {
-
-Edge::Edge(Node* nodes[2], unsigned value)
-	: Element(value)
-{
-	_nodes = nodes;
-	this->_length = this->computeVolume();
-}
-
-Edge::Edge(Node* n0, Node* n1, unsigned value)
-	: Element(value)
-{
-	_nodes = new Node*[2];
-	_nodes[0] = n0;
-	_nodes[1] = n1;
-
-	this->_length = this->computeVolume();
-}
-
-Edge::Edge(const Edge &edge)
-	: Element(edge.getValue())
-{
-	_nodes = new Node*[2];
-	_nodes[0] = edge._nodes[0];
-	_nodes[1] = edge._nodes[1];
-	_length = edge.getLength();
-}
-
-Edge::~Edge()
-{
-}
-
-double Edge::computeVolume()
-{
-	return sqrt(MathLib::sqrDist(_nodes[0]->getCoords(), _nodes[1]->getCoords()));
-}
-
-bool Edge::isEdge(unsigned idx1, unsigned idx2) const
-{
-	if (0==idx1 && 1==idx2) return true;
-	if (1==idx1 && 0==idx2) return true;
-	return false;
-}
-
-Element* Edge::clone() const
-{
-	return new Edge(*this);
-}
-
-Element* Edge::reviseElement() const
-{
-	if (_nodes[0] == _nodes[1]) {
-		return NULL;
-	}
-
-	return NULL;
-}
-
-}
-
diff --git a/MeshLib/Elements/Edge.h b/MeshLib/Elements/Edge.h
index 4d78168ff2b5756a88c6473b70abdb7e50e2aa89..598b5dae57ecae4b12648e16eb6aa1c250f49b5b 100644
--- a/MeshLib/Elements/Edge.h
+++ b/MeshLib/Elements/Edge.h
@@ -1,113 +1,24 @@
 /**
- * Copyright (c) 2012, OpenGeoSys Community (http://www.opengeosys.org)
+ * Copyright (c) 2012, OpenGeoSys Community (http://www.opengeosys.net)
  *            Distributed under a Modified BSD License.
  *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
+ *              http://www.opengeosys.org/LICENSE.txt
  *
  * \file Edge.h
  *
- * Created on 2012-05-02 by Karsten Rink
+ *  Created on  Sep 27, 2012 by Thomas Fischer
  */
 
 #ifndef EDGE_H_
 #define EDGE_H_
 
-#include <limits>
-
-#include "Element.h"
+#include "TemplateEdge.h"
 
 namespace MeshLib {
 
-class Node;
-
-/**
- * A 1d Edge or Line Element.
- * @code
- *  0--------1
- * @endcode
- */
-class Edge : public Element
-{
-public:
-	/// Constructor with an array of mesh nodes.
-	Edge(Node* nodes[2], unsigned value = 0);
-
-	/// Constructor using single nodes
-	Edge(Node* n1, Node* n2, unsigned value = 0);
-
-	/// Copy constructor
-	Edge(const Edge &edge);
-
-	/// Destructor
-	virtual ~Edge();
-
-	/// Returns the length, area or volume of a 1D, 2D or 3D element
-	double getContent() const { return _length; };
-
-	/// Returns the edge i of the element.
-	const Element* getEdge(unsigned i) const { (void)i; return NULL; };
-
-	/// Returns the face i of the element.
-	const Element* getFace(unsigned i) const { (void)i; return NULL; };
-
-	/// Compute the minimum and maximum squared edge length for this element
-	void computeSqrEdgeLengthRange(double &min, double &max) const { min = _length; max = _length; };
-
-	/// 1D elements have no edges
-	unsigned getNEdges() const { return 0; };
-
-	/// Get the number of nodes for face i.
-	unsigned getNFaceNodes(unsigned i) const { (void)i; return 0; };
+typedef TemplateEdge<1,2> Edge;
 
-	/// Get the number of faces for this element.
-	unsigned getNFaces() const { return 0; };
+}
 
-	/// Get the length of this 1d element.
-	double getLength() const { return _length; };
-
-	/// Get dimension of the mesh element.
-	unsigned getDimension() const { return 1; };
-
-	/// Get the number of neighbors for this element.
-	unsigned getNNeighbors() const { return 0; };
-
-	/// Get the number of nodes for this element.
-	virtual unsigned getNNodes() const { return 2; };
-
-	/**
-	 * Method returns the type of the element. In this case EDGE will be returned.
-	 * @return MshElemType::EDGE
-	 */
-	virtual MshElemType::type getType() const { return MshElemType::EDGE; }
-
-	/// Returns true if these two indeces form an edge and false otherwise
-	bool isEdge(unsigned i, unsigned j) const;
-
-	virtual Element* clone() const;
-
-	/**
-	 * If for instance two nodes of the element are collapsed the Edge element disappears.
-	 * @return NULL
-	 */
-	virtual Element* reviseElement() const;
-
-protected:
-	double computeVolume();
-	/// 1D elements have no edges.
-	Node* getEdgeNode(unsigned edge_id, unsigned node_id) const { (void)edge_id; (void)node_id; return NULL; };
-
-	/// 1D elements have no faces.
-	Node* getFaceNode(unsigned face_id, unsigned node_id) const { (void)face_id; (void)node_id; return NULL; };
-
-	/// Returns the ID of a face given an array of nodes (but is not applicable for edges!).
-	unsigned identifyFace(Node* [3]/*nodes[3]*/) const { return std::numeric_limits<unsigned>::max(); };
-
-	double _length;
-
-}; /* class */
-
-} /* namespace */
 
 #endif /* EDGE_H_ */
-
diff --git a/MeshLib/Elements/Element.h b/MeshLib/Elements/Element.h
index 2b896aa6a1dc3dcbef40c23a0113575a13bff4d7..3f43902829a2f99c17632e03cf6e047f4ce309c4 100644
--- a/MeshLib/Elements/Element.h
+++ b/MeshLib/Elements/Element.h
@@ -88,8 +88,12 @@ public:
 	/// Get the number of neighbors for this element.
 	virtual unsigned getNNeighbors() const = 0;
 
-	/// Get the number of nodes for this element.
-	virtual unsigned getNNodes() const = 0;
+	/**
+	 * Get the number of nodes for this element with respect to the order.
+	 * @param order (default = 1)
+	 * @return the number of nodes with respect to the order.
+	 */
+	virtual unsigned getNNodes(unsigned order = 1) const = 0;
 
 	/// Returns the position of the given node in the node array of this element.
 	virtual unsigned getNodeIDinElement(const MeshLib::Node* node) const;
diff --git a/MeshLib/Elements/Hex.h b/MeshLib/Elements/Hex.h
index 14943ec47a3b1f5f95f75c57dfa68a3c330417cd..e21fbf2c31f7bc0a9e9ab4cd63196207ea9ad810 100644
--- a/MeshLib/Elements/Hex.h
+++ b/MeshLib/Elements/Hex.h
@@ -1,117 +1,21 @@
 /**
- * Copyright (c) 2012, OpenGeoSys Community (http://www.opengeosys.org)
+ * Copyright (c) 2012, OpenGeoSys Community (http://www.opengeosys.net)
  *            Distributed under a Modified BSD License.
  *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
+ *              http://www.opengeosys.org/LICENSE.txt
  *
  * \file Hex.h
  *
- * Created on 2012-05-02 by Karsten Rink
+ *  Created on  Sep 27, 2012 by Thomas Fischer
  */
 
 #ifndef HEX_H_
 #define HEX_H_
 
-#include "Cell.h"
+#include "TemplateHex.h"
 
 namespace MeshLib {
-
-/**
- * A 3d Hexahedron Element.
- * @code
- *
- *  Hex:
- *                10
- *          7-----------6
- *         /:          /|
- *        / :         / |
- *     11/  :        /9 |
- *      /  7:       /   | 6
- *     /    : 8    /    |
- *    4-----------5     |
- *    |     :     | 2   |
- *    |     3.....|.....2
- *    |    .      |    /
- *  4 |   .       |5  /
- *    | 3.        |  / 1
- *    | .         | /
- *    |.          |/
- *    0-----------1
- *          0
- *
- * @endcode
- */
-class Hex : public Cell
-{
-public:
-	/// Constructor with an array of mesh nodes.
-	Hex(Node* nodes[8], unsigned value = 0);
-
-	/// Constructor using single mesh nodes.
-	Hex(Node* n0, Node* n1, Node* n2, Node* n3, Node* n4, Node* n5, Node* n6, Node* n7, unsigned value);
-
-	/// Copy constructor
-	Hex(const Hex &hex);
-
-	/// Destructor
-	virtual ~Hex();
-
-	/// Returns the face i of the element.
-	const Element* getFace(unsigned i) const;
-
-	/// Get the number of edges for this element.
-	unsigned getNEdges() const { return 12; };
-
-	/// Get the number of nodes for face i.
-	unsigned getNFaceNodes(unsigned i) const { (void)i; return 4; };
-
-	/// Get the number of faces for this element.
-	unsigned getNFaces() const { return 6; };
-
-	/// Get the number of neighbors for this element.
-	unsigned getNNeighbors() const { return 6; };
-
-	/// Get the number of nodes for this element.
-	virtual unsigned getNNodes() const { return 8; };
-
-	/**
-	 * Method returns the type of the element. In this case HEXAHEDRON will be returned.
-	 * @return MshElemType::HEXAHEDRON
-	 */
-	virtual MshElemType::type getType() const { return MshElemType::HEXAHEDRON; }
-
-	/// Returns true if these two indeces form an edge and false otherwise
-	bool isEdge(unsigned i, unsigned j) const;
-
-	/**
-	 * Method clone is inherited from class Element. It makes a deep copy of the Hex instance.
-	 * @return an exact copy of the object
-	 */
-	virtual Element* clone() const;
-
-	/**
-	 * Change the element type from hexahedron to a prism if two appropriate edges of the hexahedron are collapsed.
-	 * @return a prism element with nice properties or NULL
-	 */
-	virtual Element* reviseElement() const;
-
-protected:
-	/// Calculates the volume of a convex hexahedron by partitioning it into six tetrahedra.
-	double computeVolume();
-
-	/// Return a specific edge node.
-	inline Node* getEdgeNode(unsigned edge_id, unsigned node_id) const { return _nodes[_edge_nodes[edge_id][node_id]]; };
-
-	/// Returns the ID of a face given an array of nodes.
-	unsigned identifyFace(Node* nodes[3]) const;
-
-	static const unsigned _face_nodes[6][4];
-	static const unsigned _edge_nodes[12][2];
-
-}; /* class */
-
-} /* namespace */
+typedef TemplateHex<1,8> Hex;
+}
 
 #endif /* HEX_H_ */
-
diff --git a/MeshLib/Elements/Prism.h b/MeshLib/Elements/Prism.h
index 785283715fb1db3c70157c1ced2d6133ef7f4f3a..0a442bc296570f5cd65333faea62fab7319de8c7 100644
--- a/MeshLib/Elements/Prism.h
+++ b/MeshLib/Elements/Prism.h
@@ -1,122 +1,23 @@
 /**
- * Copyright (c) 2012, OpenGeoSys Community (http://www.opengeosys.org)
+ * Copyright (c) 2012, OpenGeoSys Community (http://www.opengeosys.net)
  *            Distributed under a Modified BSD License.
  *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
+ *              http://www.opengeosys.org/LICENSE.txt
  *
  * \file Prism.h
  *
- * Created on 2012-05-02 by Karsten Rink
+ *  Created on  Sep 27, 2012 by Thomas Fischer
  */
 
 #ifndef PRISM_H_
 #define PRISM_H_
 
-#include "Cell.h"
+#include "TemplatePrism.h"
 
 namespace MeshLib {
 
-/**
- * This class represents a 3d prism element. The following sketch shows the node and edge numbering.
- * @anchor PrismNodeAndEdgeNumbering
- * @code
- *            5
- *           / \
- *          / : \
- *        8/  :  \7
- *        /   :5  \
- *       /    :  6 \
- *      3-----------4
- *      |     :     |
- *      |     2     |
- *      |    . .    |
- *     3|   .   .   |4
- *      | 2.     .1 |
- *      | .       . |
- *      |.         .|
- *      0-----------1
- *            0
- *
- * @endcode
- */
-class Prism : public Cell
-{
-public:
-	/// Constructor with an array of mesh nodes.
-	Prism(Node* nodes[6], unsigned value = 0);
-
-	/// Constructor using single mesh nodes.
-	Prism(Node* n0, Node* n1, Node* n2, Node* n3, Node* n4, Node* n5, unsigned value = 0);
-
-	/// Copy constructor
-	Prism(const Prism &prism);
-
-	/// Destructor
-	virtual ~Prism();
-
-	/// Returns the face i of the element.
-	const Element* getFace(unsigned i) const;
-
-	/// Get the number of edges for this element.
-	unsigned getNEdges() const { return 9; };
-
-	/// Get the number of nodes for face i.
-	unsigned getNFaceNodes(unsigned i) const;
-
-	/// Get the number of faces for this element.
-	unsigned getNFaces() const { return 5; };
-
-	/// Get the number of neighbors for this element.
-	unsigned getNNeighbors() const { return 5; };
+typedef TemplatePrism<1,6> Prism;
 
-	/// Get the number of nodes for this element.
-	virtual unsigned getNNodes() const { return 6; };
-
-	/**
-	 * Method returns the type of the element. In this case PRISM will be returned.
-	 * @return MshElemType::PRISM
-	 */
-	virtual MshElemType::type getType() const { return MshElemType::PRISM; }
-
-	/// Returns true if these two indeces form an edge and false otherwise
-	bool isEdge(unsigned i, unsigned j) const;
-
-	/**
-	 * Method clone is inherited from class Element. It makes a deep copy of the
-	 * Hex instance employing the copy constructor of class Prism.
-	 * @return an exact copy of the object
-	 */
-	virtual Element* clone() const;
-
-	/**
-	 * This method should be called after at least two nodes of the prism
-	 * element are collapsed. As a consequence of the node collapsing an edge
-	 * of the prism will be collapsed. If one of the edges 3, 4 or 5 (see
-	 * sketch @ref PrismNodeAndEdgeNumbering) is collapsed we obtain a
-	 * pyramid. In this case the method will create the appropriate
-	 * object of class Pyramid.
-	 * @return a pyramid object or NULL
-	 */
-	virtual Element* reviseElement() const;
-
-protected:
-	/// Calculates the volume of a prism by subdividing it into three tetrahedra.
-	double computeVolume();
-
-	/// Return a specific edge node.
-	inline Node* getEdgeNode(unsigned edge_id, unsigned node_id) const { return _nodes[_edge_nodes[edge_id][node_id]]; };
-
-	/// Returns the ID of a face given an array of nodes.
-	unsigned identifyFace(Node* nodes[3]) const;
-
-	static const unsigned _face_nodes[5][4];
-	static const unsigned _edge_nodes[9][2];
-	static const unsigned _n_face_nodes[5];
-
-}; /* class */
-
-} /* namespace */
+}
 
 #endif /* PRISM_H_ */
-
diff --git a/MeshLib/Elements/Pyramid.h b/MeshLib/Elements/Pyramid.h
index 97d818035bb80325571784703fdc28838067654e..2dcaac5af6ac800e406a7d707127501e27e564f4 100644
--- a/MeshLib/Elements/Pyramid.h
+++ b/MeshLib/Elements/Pyramid.h
@@ -38,20 +38,18 @@ namespace MeshLib {
  *        0
  * @endcode
  */
-class Pyramid : public Cell
+template <unsigned ORDER, unsigned NNODES>
+class TemplatePyramid : public Cell
 {
 public:
 	/// Constructor with an array of mesh nodes.
-	Pyramid(Node* nodes[5], unsigned value = 0);
-
-	/// Constructor using single mesh nodes.
-	Pyramid(Node* n0, Node* n1, Node* n2, Node* n3, Node* n4, unsigned value = 0);
+	TemplatePyramid(Node* nodes[NNODES], unsigned value = 0);
 
 	/// Copy constructor
-	Pyramid(const Pyramid &pyramid);
+	TemplatePyramid(const TemplatePyramid<ORDER,NNODES> &pyramid);
 
 	/// Destructor
-	virtual ~Pyramid();
+	virtual ~TemplatePyramid();
 
 	/// Returns the face i of the element.
 	const Element* getFace(unsigned i) const;
@@ -69,7 +67,10 @@ public:
 	unsigned getNNeighbors() const { return 5; };
 
 	/// Get the number of nodes for this element.
-	virtual unsigned getNNodes() const { return 5; };
+	virtual unsigned getNNodes(unsigned order = 1) const
+	{
+		return order == ORDER ? NNODES : 5;
+	}
 
 	/**
 	 * Method returns the type of the element. In this case PYRAMID will be returned.
@@ -82,7 +83,7 @@ public:
 
 	/**
 	 * Method clone is inherited from class Element. It makes a deep copy of the
-	 * Pyramid instance employing the copy constructor of class Pyramid.
+	 * TemplatePyramid instance employing the copy constructor of class TemplatePyramid.
 	 * @return an exact copy of the object
 	 */
 	virtual Element* clone() const;
@@ -114,7 +115,11 @@ protected:
 
 }; /* class */
 
+typedef TemplatePyramid<1,5> Pyramid;
+
 } /* namespace */
 
+#include "Pyramid.hpp"
+
 #endif /* PYRAMID_H_ */
 
diff --git a/MeshLib/Elements/Pyramid.cpp b/MeshLib/Elements/Pyramid.hpp
similarity index 58%
rename from MeshLib/Elements/Pyramid.cpp
rename to MeshLib/Elements/Pyramid.hpp
index 89bbe6753bb36738d6ba520f0e832353e41bba39..2ab54693881b216218dfd9e7672bbb127bf43144 100644
--- a/MeshLib/Elements/Pyramid.cpp
+++ b/MeshLib/Elements/Pyramid.hpp
@@ -5,12 +5,11 @@
  *              http://www.opengeosys.org/project/license
  *
  *
- * \file Pyramid.cpp
+ * \file Pyramid.hpp
  *
  * Created on 2012-05-02 by Karsten Rink
  */
 
-#include "Pyramid.h"
 #include "Node.h"
 #include "Tri.h"
 #include "Tet.h"
@@ -20,7 +19,8 @@
 
 namespace MeshLib {
 
-const unsigned Pyramid::_face_nodes[5][4] =
+template <unsigned ORDER, unsigned NNODES>
+const unsigned TemplatePyramid<ORDER,NNODES>::_face_nodes[5][4] =
 {
 	{0, 1, 4, 99}, // Face 0
 	{1, 2, 4, 99}, // Face 1
@@ -29,7 +29,8 @@ const unsigned Pyramid::_face_nodes[5][4] =
 	{0, 3, 2,  1}  // Face 4
 };
 
-const unsigned Pyramid::_edge_nodes[8][2] =
+template <unsigned ORDER, unsigned NNODES>
+const unsigned TemplatePyramid<ORDER,NNODES>::_edge_nodes[8][2] =
 {
 	{0, 1}, // Edge 0
 	{1, 2}, // Edge 1
@@ -41,28 +42,15 @@ const unsigned Pyramid::_edge_nodes[8][2] =
 	{3, 4}  // Edge 7
 };
 
-const unsigned Pyramid::_n_face_nodes[5] = { 3, 3, 3, 3, 4 };
+template <unsigned ORDER, unsigned NNODES>
+const unsigned TemplatePyramid<ORDER,NNODES>::_n_face_nodes[5] = { 3, 3, 3, 3, 4 };
 
-
-Pyramid::Pyramid(Node* nodes[5], unsigned value)
+template <unsigned ORDER, unsigned NNODES>
+TemplatePyramid<ORDER,NNODES>::TemplatePyramid(Node* nodes[NNODES], unsigned value)
 	: Cell(value)
 {
 	_nodes = nodes;
-	_neighbors = new Element*[5];
-	for (unsigned i=0; i<5; i++)
-		_neighbors[i] = NULL;
-	this->_volume = this->computeVolume();
-}
 
-Pyramid::Pyramid(Node* n0, Node* n1, Node* n2, Node* n3, Node* n4, unsigned value)
-	: Cell(value)
-{
-	_nodes = new Node*[5];
-	_nodes[0] = n0;
-	_nodes[1] = n1;
-	_nodes[2] = n2;
-	_nodes[3] = n3;
-	_nodes[4] = n4;
 	_neighbors = new Element*[5];
 	for (unsigned i=0; i<5; i++)
 		_neighbors[i] = NULL;
@@ -70,30 +58,37 @@ Pyramid::Pyramid(Node* n0, Node* n1, Node* n2, Node* n3, Node* n4, unsigned valu
 	this->_volume = this->computeVolume();
 }
 
-Pyramid::Pyramid(const Pyramid &pyramid)
+template <unsigned ORDER, unsigned NNODES>
+TemplatePyramid<ORDER,NNODES>::TemplatePyramid(const TemplatePyramid<ORDER,NNODES> &pyramid)
 	: Cell(pyramid.getValue())
 {
-	_nodes = new Node*[5];
-	_neighbors = new Element*[5];
-	for (unsigned i=0; i<5; i++)
-	{
+	_nodes = new Node*[NNODES];
+	for (unsigned i=0; i<NNODES; i++) {
 		_nodes[i] = pyramid._nodes[i];
+	}
+
+	_neighbors = new Element*[5];
+	for (unsigned i=0; i<5; i++) {
 		_neighbors[i] = pyramid._neighbors[i];
 	}
+
 	_volume = pyramid.getVolume();
 }
 
-Pyramid::~Pyramid()
+template <unsigned ORDER, unsigned NNODES>
+TemplatePyramid<ORDER,NNODES>::~TemplatePyramid()
 {
 }
 
-double Pyramid::computeVolume()
+template <unsigned ORDER, unsigned NNODES>
+double TemplatePyramid<ORDER,NNODES>::computeVolume()
 {
 	return MathLib::calcTetrahedronVolume(_nodes[0]->getCoords(), _nodes[1]->getCoords(), _nodes[2]->getCoords(), _nodes[4]->getCoords())
 		 + MathLib::calcTetrahedronVolume(_nodes[2]->getCoords(), _nodes[3]->getCoords(), _nodes[0]->getCoords(), _nodes[4]->getCoords());
 }
 
-const Element* Pyramid::getFace(unsigned i) const
+template <unsigned ORDER, unsigned NNODES>
+const Element* TemplatePyramid<ORDER,NNODES>::getFace(unsigned i) const
 {
 	if (i<this->getNFaces())
 	{
@@ -111,7 +106,8 @@ const Element* Pyramid::getFace(unsigned i) const
 	return NULL;
 }
 
-unsigned Pyramid::getNFaceNodes(unsigned i) const
+template <unsigned ORDER, unsigned NNODES>
+unsigned TemplatePyramid<ORDER,NNODES>::getNFaceNodes(unsigned i) const
 {
 	if (i<5)
 		return _n_face_nodes[i];
@@ -119,7 +115,8 @@ unsigned Pyramid::getNFaceNodes(unsigned i) const
 	return 0;
 }
 
-bool Pyramid::isEdge(unsigned idx1, unsigned idx2) const
+template <unsigned ORDER, unsigned NNODES>
+bool TemplatePyramid<ORDER,NNODES>::isEdge(unsigned idx1, unsigned idx2) const
 {
 	for (unsigned i(0); i<8; i++)
 	{
@@ -129,12 +126,14 @@ bool Pyramid::isEdge(unsigned idx1, unsigned idx2) const
 	return false;
 }
 
-Element* Pyramid::clone() const
+template <unsigned ORDER, unsigned NNODES>
+Element* TemplatePyramid<ORDER,NNODES>::clone() const
 {
-	return new Pyramid(*this);
+	return new TemplatePyramid<ORDER,NNODES>(*this);
 }
 
-unsigned Pyramid::identifyFace(Node* nodes[3]) const
+template <unsigned ORDER, unsigned NNODES>
+unsigned TemplatePyramid<ORDER,NNODES>::identifyFace(Node* nodes[3]) const
 {
 	for (unsigned i=0; i<5; i++)
 	{
@@ -149,17 +148,23 @@ unsigned Pyramid::identifyFace(Node* nodes[3]) const
 	return std::numeric_limits<unsigned>::max();
 }
 
-Element* Pyramid::reviseElement() const
+template <unsigned ORDER, unsigned NNODES>
+Element* TemplatePyramid<ORDER,NNODES>::reviseElement() const
 {
 	// try to create tetrahedron
 	if (_nodes[_edge_nodes[0][0]] == _nodes[_edge_nodes[0][1]]
 		|| _nodes[_edge_nodes[1][0]] == _nodes[_edge_nodes[1][1]]) {
-		return new Tet(_nodes[0], _nodes[2], _nodes[3], _nodes[4], _value);
+		Node** tet_nodes(new Node*[4]);
+		for (unsigned k(0); k<4; k++) tet_nodes[k] = _nodes[k];
+		return new Tet(tet_nodes, _value);
 	}
 
 	if (_nodes[_edge_nodes[2][0]] == _nodes[_edge_nodes[2][1]]
 	|| _nodes[_edge_nodes[3][0]] == _nodes[_edge_nodes[3][1]]) {
-		return new Tet(_nodes[0], _nodes[1], _nodes[2], _nodes[4], _value);
+		Node** tet_nodes(new Node*[4]);
+		for (unsigned k(0); k<3; k++) tet_nodes[k] = _nodes[k];
+		tet_nodes[3] = _nodes[4];
+		return new Tet(tet_nodes, _value);
 	}
 
 	return NULL;
diff --git a/MeshLib/Elements/Quad.h b/MeshLib/Elements/Quad.h
index e5a8d798b4cfa3793775805a4a9ddaea4199fbb5..0cc34b19bc9244384da5985538d827c44a7e261b 100644
--- a/MeshLib/Elements/Quad.h
+++ b/MeshLib/Elements/Quad.h
@@ -1,103 +1,23 @@
 /**
- * Copyright (c) 2012, OpenGeoSys Community (http://www.opengeosys.org)
+ * Copyright (c) 2012, OpenGeoSys Community (http://www.opengeosys.net)
  *            Distributed under a Modified BSD License.
  *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
+ *              http://www.opengeosys.org/LICENSE.txt
  *
  * \file Quad.h
  *
- * Created on 2012-05-02 by Karsten Rink
+ *  Created on  Sep 27, 2012 by Thomas Fischer
  */
 
 #ifndef QUAD_H_
 #define QUAD_H_
 
-#include "Face.h"
+#include "TemplateQuad.h"
 
 namespace MeshLib {
 
-/**
- * This class represents a 2d quadliteral element. The following sketch shows the node and edge numbering.
- * @anchor QuadNodeAndEdgeNumbering
- * @code
- *              2
- *        3-----------2
- *        |           |
- *        |           |
- *       3|           |1
- *        |           |
- *        |           |
- *        0-----------1
- *              0
- * @endcode
- */
-class Quad : public Face
-{
-public:
-	/// Constructor with an array of mesh nodes.
-	Quad(Node* nodes[4], unsigned value = 0);
-
-	/// Constructor using single mesh nodes.
-	Quad(Node* n0, Node* n1, Node* n2, Node* n3, unsigned value = 0);
-
-	/// Copy constructor
-	Quad(const Quad &quad);
-
-	/// Destructor
-	virtual ~Quad();
-
-	/// Get the number of edges for this element.
-	unsigned getNEdges() const { return 4; };
-
-	/// Get the number of neighbors for this element.
-	unsigned getNNeighbors() const { return 4; };
-
-	/// Get the number of nodes for this element.
-	virtual unsigned getNNodes() const { return 4; };
-
-	/**
-	 * Method returns the type of the element. In this case QUAD will be returned.
-	 * @return MshElemType::QUAD
-	 */
-	virtual MshElemType::type getType() const { return MshElemType::QUAD; }
+typedef TemplateQuad<1,4> Quad;
 
-	/// Returns true if these two indeces form an edge and false otherwise
-	bool isEdge(unsigned i, unsigned j) const;
-
-	/**
-	 * Method clone is inherited from class Element. It makes a deep copy of the Quad instance.
-	 * @return an exact copy of the object
-	 */
-	virtual Element* clone() const;
-
-	/**
-	 * This method should be called after at least two nodes of the quad
-	 * element are collapsed. As a consequence of the node collapsing an edge
-	 * of the quad will be collapsed. If one of the edges (see
-	 * sketch @ref PyramidNodeAndEdgeNumbering) is collapsed we obtain a
-	 * triangle. In this case the method will create the appropriate
-	 * object of class Tri.
-	 * @return a Tri object or NULL
-	 */
-	virtual Element* reviseElement() const;
-
-protected:
-	/// Calculates the area of a convex quadliteral by dividing it into two triangles.
-	double computeVolume();
-
-protected:
-	/// Return a specific edge node.
-	inline Node* getEdgeNode(unsigned edge_id, unsigned node_id) const { return _nodes[_edge_nodes[edge_id][node_id]]; };
-
-	/// Returns the ID of a face given an array of nodes.
-	unsigned identifyFace(Node* nodes[3]) const;
-
-	static const unsigned _edge_nodes[4][2];
-
-}; /* class */
-
-} /* namespace */
+}
 
 #endif /* QUAD_H_ */
-
diff --git a/MeshLib/Elements/TemplateEdge.h b/MeshLib/Elements/TemplateEdge.h
new file mode 100644
index 0000000000000000000000000000000000000000..149a94d739bdbc81d39dbbbac69a2263cc4955df
--- /dev/null
+++ b/MeshLib/Elements/TemplateEdge.h
@@ -0,0 +1,137 @@
+/**
+ * Copyright (c) 2012, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ *
+ * \file Edge.h
+ *
+ * Created on 2012-05-02 by Karsten Rink
+ */
+
+#ifndef TEMPLATEEDGE_H_
+#define TEMPLATEEDGE_H_
+
+#include <limits>
+
+#include "Element.h"
+#include "Node.h"
+
+#include "MathTools.h"
+
+
+namespace MeshLib {
+
+/**
+ * A 1d Edge or Line Element.
+ * @code
+ *  0--------1
+ * @endcode
+ */
+template <unsigned ORDER, unsigned NNODES>
+class TemplateEdge : public Element
+{
+public:
+	/// Constructor with an array of mesh nodes.
+	TemplateEdge(Node* nodes[NNODES], unsigned value = 0);
+
+	/// Copy constructor
+	TemplateEdge(const TemplateEdge<ORDER, NNODES> &edge);
+
+	/// Destructor
+	virtual ~TemplateEdge();
+
+	/// Returns the length, area or volume of a 1D, 2D or 3D element
+	double getContent() const { return _length; };
+
+	/// Returns the edge i of the element.
+	const Element* getEdge(unsigned i) const { (void)i; return NULL; };
+
+	/// Returns the face i of the element.
+	const Element* getFace(unsigned i) const { (void)i; return NULL; };
+
+	/// Compute the minimum and maximum squared edge length for this element
+	void computeSqrEdgeLengthRange(double &min, double &max) const { min = _length; max = _length; };
+
+	/// 1D elements have no edges
+	unsigned getNEdges() const { return 0; };
+
+	/// Get the number of nodes for face i.
+	unsigned getNFaceNodes(unsigned i) const { (void)i; return 0; };
+
+	/// Get the number of faces for this element.
+	unsigned getNFaces() const { return 0; };
+
+	/// Get the length of this 1d element.
+	double getLength() const { return _length; };
+
+	/// Get dimension of the mesh element.
+	unsigned getDimension() const { return 1; };
+
+	/// Get the number of neighbors for this element.
+	unsigned getNNeighbors() const { return 0; };
+
+	/// Get the number of nodes for this element.
+	virtual unsigned getNNodes(unsigned order) const
+	{
+		return order == ORDER ? NNODES : 2;
+	}
+
+	/**
+	 * Method returns the type of the element. In this case EDGE will be returned.
+	 * @return MshElemType::EDGE
+	 */
+	virtual MshElemType::type getType() const { return MshElemType::EDGE; }
+
+	/// Returns true if these two indices form an edge and false otherwise
+	bool isEdge(unsigned idx1, unsigned idx2) const
+	{
+		if (0==idx1 && 1==idx2) return true;
+		if (1==idx1 && 0==idx2) return true;
+		return false;
+	}
+
+	virtual Element* clone() const
+	{
+		return new TemplateEdge<ORDER,NNODES>(*this);
+	}
+
+	/**
+	 * If for instance two nodes of the element are collapsed the Edge element disappears.
+	 * @return NULL
+	 */
+	virtual Element* reviseElement() const
+	{
+		if (_nodes[0] == _nodes[1]) {
+			return NULL;
+		}
+
+		return NULL;
+	}
+
+protected:
+	double computeVolume()
+	{
+		return sqrt(MathLib::sqrDist(_nodes[0]->getCoords(), _nodes[1]->getCoords()));
+	}
+
+	/// 1D elements have no edges.
+	Node* getEdgeNode(unsigned edge_id, unsigned node_id) const { (void)edge_id; (void)node_id; return NULL; };
+
+	/// 1D elements have no faces.
+	Node* getFaceNode(unsigned face_id, unsigned node_id) const { (void)face_id; (void)node_id; return NULL; };
+
+	/// Returns the ID of a face given an array of nodes (but is not applicable for edges!).
+	unsigned identifyFace(Node* [3]/*nodes[3]*/) const { return std::numeric_limits<unsigned>::max(); };
+
+	double _length;
+
+}; /* class */
+
+} /* namespace */
+
+#include "TemplateEdge.hpp"
+
+#endif /* TEMPLATEEDGE_H_ */
+
diff --git a/MeshLib/Elements/TemplateEdge.hpp b/MeshLib/Elements/TemplateEdge.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..265b1010bd1fcb862f6b0081bbc1e6dcee483ddb
--- /dev/null
+++ b/MeshLib/Elements/TemplateEdge.hpp
@@ -0,0 +1,41 @@
+/**
+ * Copyright (c) 2012, 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 TemplateEdge.hpp
+ *
+ *  Created on  Sep 27, 2012 by Thomas Fischer
+ */
+
+#ifndef TEMPLATEEDGE_HPP_
+#define TEMPLATEEDGE_HPP_
+
+namespace MeshLib {
+
+template <unsigned ORDER, unsigned NNODES>
+TemplateEdge<ORDER,NNODES>::TemplateEdge(Node* nodes[NNODES], unsigned value) :
+	Element(value)
+{
+	_nodes = nodes;
+	this->_length = this->computeVolume();
+}
+
+template <unsigned ORDER, unsigned NNODES>
+TemplateEdge<ORDER,NNODES>::TemplateEdge(const TemplateEdge<ORDER, NNODES> &edge) :
+	Element(edge.getValue())
+{
+	_nodes = new Node*[NNODES];
+	for (unsigned k(0); k<NNODES; k++)
+		_nodes[k] = edge._nodes[k];
+	_length = edge.getLength();
+}
+
+template <unsigned ORDER, unsigned NNODES>
+TemplateEdge<ORDER,NNODES>::~TemplateEdge()
+{}
+
+} // namespace MeshLib
+
+#endif /* TEMPLATEEDGE_HPP_ */
diff --git a/MeshLib/Elements/TemplateHex.h b/MeshLib/Elements/TemplateHex.h
new file mode 100644
index 0000000000000000000000000000000000000000..8b40b39226d0784a7bf93a4b0a367f2897bdb665
--- /dev/null
+++ b/MeshLib/Elements/TemplateHex.h
@@ -0,0 +1,120 @@
+/**
+ * Copyright (c) 2012, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ *
+ * \file TemplateHex.h
+ *
+ * Created on 2012-05-02 by Karsten Rink
+ */
+
+#ifndef TEMPLATEHEX_H_
+#define TEMPLATEHEX_H_
+
+#include "Cell.h"
+
+namespace MeshLib {
+
+/**
+ * A 3d Hexahedron Element.
+ * @code
+ *
+ *  Hex:
+ *                10
+ *          7-----------6
+ *         /:          /|
+ *        / :         / |
+ *     11/  :        /9 |
+ *      /  7:       /   | 6
+ *     /    : 8    /    |
+ *    4-----------5     |
+ *    |     :     | 2   |
+ *    |     3.....|.....2
+ *    |    .      |    /
+ *  4 |   .       |5  /
+ *    | 3.        |  / 1
+ *    | .         | /
+ *    |.          |/
+ *    0-----------1
+ *          0
+ *
+ * @endcode
+ */
+template <unsigned ORDER, unsigned NNODES>
+class TemplateHex : public Cell
+{
+public:
+	/// Constructor with an array of mesh nodes.
+	TemplateHex(Node* nodes[8], unsigned value = 0);
+
+	/// Copy constructor
+	TemplateHex(const TemplateHex<ORDER,NNODES> &hex);
+
+	/// Destructor
+	virtual ~TemplateHex();
+
+	/// Returns the face i of the element.
+	const Element* getFace(unsigned i) const;
+
+	/// Get the number of edges for this element.
+	unsigned getNEdges() const { return 12; };
+
+	/// Get the number of nodes for face i.
+	unsigned getNFaceNodes(unsigned i) const { (void)i; return 4; };
+
+	/// Get the number of faces for this element.
+	unsigned getNFaces() const { return 6; };
+
+	/// Get the number of neighbors for this element.
+	unsigned getNNeighbors() const { return 6; };
+
+	/// Get the number of nodes for this element.
+	virtual unsigned getNNodes(unsigned order = 1) const
+	{
+		return order == ORDER ? NNODES : 8;
+	}
+
+	/**
+	 * Method returns the type of the element. In this case HEXAHEDRON will be returned.
+	 * @return MshElemType::HEXAHEDRON
+	 */
+	virtual MshElemType::type getType() const { return MshElemType::HEXAHEDRON; }
+
+	/// Returns true if these two indices form an edge and false otherwise
+	bool isEdge(unsigned i, unsigned j) const;
+
+	/**
+	 * Method clone is inherited from class Element. It makes a deep copy of the Hex instance.
+	 * @return an exact copy of the object
+	 */
+	virtual Element* clone() const;
+
+	/**
+	 * Change the element type from hexahedron to a prism if two appropriate edges of the hexahedron are collapsed.
+	 * @return a prism element with nice properties or NULL
+	 */
+	virtual Element* reviseElement() const;
+
+protected:
+	/// Calculates the volume of a convex hexahedron by partitioning it into six tetrahedra.
+	double computeVolume();
+
+	/// Return a specific edge node.
+	inline Node* getEdgeNode(unsigned edge_id, unsigned node_id) const { return _nodes[_edge_nodes[edge_id][node_id]]; };
+
+	/// Returns the ID of a face given an array of nodes.
+	unsigned identifyFace(Node* nodes[3]) const;
+
+	static const unsigned _face_nodes[6][4];
+	static const unsigned _edge_nodes[12][2];
+
+}; /* class */
+
+} /* namespace */
+
+#include "TemplateHex.hpp"
+
+#endif /* TEMPLATEHEX_H_ */
+
diff --git a/MeshLib/Elements/Hex.cpp b/MeshLib/Elements/TemplateHex.hpp
similarity index 50%
rename from MeshLib/Elements/Hex.cpp
rename to MeshLib/Elements/TemplateHex.hpp
index 2ec5bf0a108bb6eb511677c2db295787f5ad3c5e..92e1a821868c5321c7b8a0cbba5241c38dd9e999 100644
--- a/MeshLib/Elements/Hex.cpp
+++ b/MeshLib/Elements/TemplateHex.hpp
@@ -5,22 +5,21 @@
  *              http://www.opengeosys.org/project/license
  *
  *
- * \file Hex.cpp
+ * \file TemplateHex.hpp
  *
  * Created on 2012-05-02 by Karsten Rink
  */
 
-#include "Hex.h"
 #include "Node.h"
 #include "Quad.h"
 #include "Prism.h"
 
 #include "MathTools.h"
 
-
 namespace MeshLib {
 
-const unsigned Hex::_face_nodes[6][4] =
+template <unsigned ORDER, unsigned NNODES>
+const unsigned TemplateHex<ORDER,NNODES>::_face_nodes[6][4] =
 {
 	{0, 3, 2, 1}, // Face 0
 	{0, 1, 5, 4}, // Face 1
@@ -30,7 +29,8 @@ const unsigned Hex::_face_nodes[6][4] =
 	{4, 5, 6, 7}  // Face 5
 };
 
-const unsigned Hex::_edge_nodes[12][2] =
+template <unsigned ORDER, unsigned NNODES>
+const unsigned TemplateHex<ORDER,NNODES>::_edge_nodes[12][2] =
 {
 	{0, 1}, // Edge 0
 	{1, 2}, // Edge 1
@@ -46,52 +46,41 @@ const unsigned Hex::_edge_nodes[12][2] =
 	{4, 7}  // Edge 11
 };
 
-
-Hex::Hex(Node* nodes[8], unsigned value)
+template <unsigned ORDER, unsigned NNODES>
+TemplateHex<ORDER,NNODES>::TemplateHex(Node* nodes[8], unsigned value)
 	: Cell(value)
 {
 	_nodes = nodes;
-	_neighbors = new Element*[6];
-	for (unsigned i=0; i<6; i++)
-		_neighbors[i] = NULL;
-	this->_volume = this->computeVolume();
-}
 
-Hex::Hex(Node* n0, Node* n1, Node* n2, Node* n3, Node* n4, Node* n5, Node* n6, Node* n7, unsigned value)
-	: Cell(value)
-{
-	_nodes = new Node*[8];
-	_nodes[0] = n0;
-	_nodes[1] = n1;
-	_nodes[2] = n2;
-	_nodes[3] = n3;
-	_nodes[4] = n4;
-	_nodes[5] = n5;
-	_nodes[6] = n6;
-	_nodes[7] = n7;
 	_neighbors = new Element*[6];
 	for (unsigned i=0; i<6; i++)
 		_neighbors[i] = NULL;
+
 	this->_volume = this->computeVolume();
 }
 
-Hex::Hex(const Hex &hex)
+template <unsigned ORDER, unsigned NNODES>
+TemplateHex<ORDER,NNODES>::TemplateHex(const TemplateHex<ORDER,NNODES> &hex)
 	: Cell(hex.getValue())
 {
-	_nodes = new Node*[8];
-	for (unsigned i=0; i<8; i++)
+	_nodes = new Node*[NNODES];
+	for (unsigned i=0; i<NNODES; i++)
 		_nodes[i] = hex._nodes[i];
+
 	_neighbors = new Element*[6];
 	for (unsigned i=0; i<6; i++)
 		_neighbors[i] = hex._neighbors[i];
+
 	_volume = hex.getVolume();
 }
 
-Hex::~Hex()
+template <unsigned ORDER, unsigned NNODES>
+TemplateHex<ORDER,NNODES>::~TemplateHex()
 {
 }
 
-double Hex::computeVolume()
+template <unsigned ORDER, unsigned NNODES>
+double TemplateHex<ORDER,NNODES>::computeVolume()
 {
 	return MathLib::calcTetrahedronVolume(_nodes[4]->getCoords(), _nodes[7]->getCoords(), _nodes[5]->getCoords(), _nodes[0]->getCoords())
 		 + MathLib::calcTetrahedronVolume(_nodes[5]->getCoords(), _nodes[3]->getCoords(), _nodes[1]->getCoords(), _nodes[0]->getCoords())
@@ -101,7 +90,8 @@ double Hex::computeVolume()
 		 + MathLib::calcTetrahedronVolume(_nodes[3]->getCoords(), _nodes[7]->getCoords(), _nodes[5]->getCoords(), _nodes[2]->getCoords());
 }
 
-const Element* Hex::getFace(unsigned i) const
+template <unsigned ORDER, unsigned NNODES>
+const Element* TemplateHex<ORDER,NNODES>::getFace(unsigned i) const
 {
 	if (i<this->getNFaces())
 	{
@@ -115,7 +105,8 @@ const Element* Hex::getFace(unsigned i) const
 	return NULL;
 }
 
-bool Hex::isEdge(unsigned idx1, unsigned idx2) const
+template <unsigned ORDER, unsigned NNODES>
+bool TemplateHex<ORDER,NNODES>::isEdge(unsigned idx1, unsigned idx2) const
 {
 	for (unsigned i(0); i<12; i++)
 	{
@@ -125,12 +116,14 @@ bool Hex::isEdge(unsigned idx1, unsigned idx2) const
 	return false;
 }
 
-Element* Hex::clone() const
+template <unsigned ORDER, unsigned NNODES>
+Element* TemplateHex<ORDER,NNODES>::clone() const
 {
-	return new Hex(*this);
+	return new TemplateHex<ORDER,NNODES>(*this);
 }
 
-unsigned Hex::identifyFace(Node* nodes[3]) const
+template <unsigned ORDER, unsigned NNODES>
+unsigned TemplateHex<ORDER,NNODES>::identifyFace(Node* nodes[3]) const
 {
 	for (unsigned i=0; i<6; i++)
 	{
@@ -145,7 +138,8 @@ unsigned Hex::identifyFace(Node* nodes[3]) const
 	return std::numeric_limits<unsigned>::max();
 }
 
-Element* Hex::reviseElement() const
+template <unsigned ORDER, unsigned NNODES>
+Element* TemplateHex<ORDER,NNODES>::reviseElement() const
 {
 	std::vector<size_t> collapsed_edges;
 	for (size_t edge(0); edge<getNEdges(); edge++) {
@@ -155,35 +149,91 @@ Element* Hex::reviseElement() const
 	}
 
 	if (collapsed_edges.size() == 1) {
-		std::cerr << "[Hex::reviseElement()] collapsing of one edge in hexahedron not handled" << std::endl;
+		std::cerr << "[TemplateHex<ORDER,NNODES>::reviseElement()] collapsing of one edge in hexahedron not handled" << std::endl;
 		return NULL;
 	}
 
 	if (collapsed_edges.size() == 2) {
 		// try to create a prism out of the hex
 		if (collapsed_edges[0] == 0 && collapsed_edges[1] == 2) {
-			return new Prism(_nodes[0],_nodes[4], _nodes[5], _nodes[3], _nodes[7], _nodes[6], _value);
+			Node** prism_nodes(new Node*[6]);
+			prism_nodes[0] = _nodes[0];
+			prism_nodes[1] = _nodes[4];
+			prism_nodes[2] = _nodes[5];
+			prism_nodes[3] = _nodes[3];
+			prism_nodes[4] = _nodes[7];
+			prism_nodes[5] = _nodes[6];
+			return new Prism(prism_nodes, _value);
 		}
 		if (collapsed_edges[0] == 1 && collapsed_edges[1] == 3) {
-			return new Prism(_nodes[0],_nodes[4], _nodes[7], _nodes[1], _nodes[5], _nodes[6], _value);
+			Node** prism_nodes(new Node*[6]);
+			prism_nodes[0] = _nodes[0];
+			prism_nodes[1] = _nodes[4];
+			prism_nodes[2] = _nodes[7];
+			prism_nodes[3] = _nodes[1];
+			prism_nodes[4] = _nodes[5];
+			prism_nodes[5] = _nodes[6];
+			return new Prism(prism_nodes, _value);
 		}
 		if (collapsed_edges[0] == 4 && collapsed_edges[1] == 5) {
-			return new Prism(_nodes[0],_nodes[7], _nodes[3], _nodes[1], _nodes[6], _nodes[2], _value);
+			Node** prism_nodes(new Node*[6]);
+			prism_nodes[0] = _nodes[0];
+			prism_nodes[1] = _nodes[7];
+			prism_nodes[2] = _nodes[3];
+			prism_nodes[3] = _nodes[1];
+			prism_nodes[4] = _nodes[6];
+			prism_nodes[5] = _nodes[2];
+			return new Prism(prism_nodes, _value);
 		}
 		if (collapsed_edges[0] == 5 && collapsed_edges[1] == 6) {
-			return new Prism(_nodes[0],_nodes[1], _nodes[4], _nodes[3], _nodes[2], _nodes[7], _value);
+			Node** prism_nodes(new Node*[6]);
+			prism_nodes[0] = _nodes[0];
+			prism_nodes[1] = _nodes[1];
+			prism_nodes[2] = _nodes[4];
+			prism_nodes[3] = _nodes[3];
+			prism_nodes[4] = _nodes[2];
+			prism_nodes[5] = _nodes[7];
+			return new Prism(prism_nodes, _value);
 		}
 		if (collapsed_edges[0] == 6 && collapsed_edges[1] == 7) {
-			return new Prism(_nodes[0],_nodes[3], _nodes[4], _nodes[1], _nodes[2], _nodes[5], _value);
+			Node** prism_nodes(new Node*[6]);
+			prism_nodes[0] = _nodes[0];
+			prism_nodes[1] = _nodes[3];
+			prism_nodes[2] = _nodes[4];
+			prism_nodes[3] = _nodes[1];
+			prism_nodes[4] = _nodes[2];
+			prism_nodes[5] = _nodes[5];
+			return new Prism(prism_nodes, _value);
 		}
 		if (collapsed_edges[0] == 7 && collapsed_edges[1] == 4) {
-			return new Prism(_nodes[0],_nodes[1], _nodes[5], _nodes[3], _nodes[2], _nodes[6], _value);
+			Node** prism_nodes(new Node*[6]);
+			prism_nodes[0] = _nodes[0];
+			prism_nodes[1] = _nodes[1];
+			prism_nodes[2] = _nodes[5];
+			prism_nodes[3] = _nodes[3];
+			prism_nodes[4] = _nodes[2];
+			prism_nodes[5] = _nodes[6];
+			return new Prism(prism_nodes, _value);
 		}
 		if (collapsed_edges[0] == 8 && collapsed_edges[1] == 10) {
-			return new Prism(_nodes[0],_nodes[1], _nodes[4], _nodes[3], _nodes[2], _nodes[7], _value);
+			Node** prism_nodes(new Node*[6]);
+			prism_nodes[0] = _nodes[0];
+			prism_nodes[1] = _nodes[1];
+			prism_nodes[2] = _nodes[4];
+			prism_nodes[3] = _nodes[3];
+			prism_nodes[4] = _nodes[2];
+			prism_nodes[5] = _nodes[7];
+			return new Prism(prism_nodes, _value);
 		}
 		if (collapsed_edges[0] == 9 && collapsed_edges[1] == 11) {
-			return new Prism(_nodes[0],_nodes[3], _nodes[4], _nodes[1], _nodes[2], _nodes[5], _value);
+			Node** prism_nodes(new Node*[6]);
+			prism_nodes[0] = _nodes[0];
+			prism_nodes[1] = _nodes[3];
+			prism_nodes[2] = _nodes[4];
+			prism_nodes[3] = _nodes[1];
+			prism_nodes[4] = _nodes[2];
+			prism_nodes[5] = _nodes[5];
+			return new Prism(prism_nodes, _value);
 		}
 		return NULL;
 	}
diff --git a/MeshLib/Elements/TemplatePrism.h b/MeshLib/Elements/TemplatePrism.h
new file mode 100644
index 0000000000000000000000000000000000000000..84b0e9752b4a2a4a3585683339cd18f4ba81a955
--- /dev/null
+++ b/MeshLib/Elements/TemplatePrism.h
@@ -0,0 +1,125 @@
+/**
+ * Copyright (c) 2012, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ *
+ * \file TemplatePrism.h
+ *
+ * Created on 2012-05-02 by Karsten Rink
+ */
+
+#ifndef TEMPLATEPRISM_H_
+#define TEMPLATEPRISM_H_
+
+#include "Cell.h"
+
+namespace MeshLib {
+
+/**
+ * This class represents a 3d prism element. The following sketch shows the node and edge numbering.
+ * @anchor PrismNodeAndEdgeNumbering
+ * @code
+ *            5
+ *           / \
+ *          / : \
+ *        8/  :  \7
+ *        /   :5  \
+ *       /    :  6 \
+ *      3-----------4
+ *      |     :     |
+ *      |     2     |
+ *      |    . .    |
+ *     3|   .   .   |4
+ *      | 2.     .1 |
+ *      | .       . |
+ *      |.         .|
+ *      0-----------1
+ *            0
+ *
+ * @endcode
+ */
+template <unsigned ORDER, unsigned NNODES>
+class TemplatePrism : public Cell
+{
+public:
+	/// Constructor with an array of mesh nodes.
+	TemplatePrism(Node* nodes[6], unsigned value = 0);
+
+	/// Copy constructor
+	TemplatePrism(const TemplatePrism<ORDER,NNODES> &prism);
+
+	/// Destructor
+	virtual ~TemplatePrism();
+
+	/// Returns the face i of the element.
+	const Element* getFace(unsigned i) const;
+
+	/// Get the number of edges for this element.
+	unsigned getNEdges() const { return 9; };
+
+	/// Get the number of nodes for face i.
+	unsigned getNFaceNodes(unsigned i) const;
+
+	/// Get the number of faces for this element.
+	unsigned getNFaces() const { return 5; };
+
+	/// Get the number of neighbors for this element.
+	unsigned getNNeighbors() const { return 5; };
+
+	/// Get the number of nodes for this element.
+	virtual unsigned getNNodes(unsigned order) const
+	{
+		return order == ORDER ? NNODES : 6;
+	}
+
+	/**
+	 * Method returns the type of the element. In this case PRISM will be returned.
+	 * @return MshElemType::PRISM
+	 */
+	virtual MshElemType::type getType() const { return MshElemType::PRISM; }
+
+	/// Returns true if these two indeces form an edge and false otherwise
+	bool isEdge(unsigned i, unsigned j) const;
+
+	/**
+	 * Method clone is inherited from class Element. It makes a deep copy of the
+	 * Hex instance employing the copy constructor of class TemplatePrism.
+	 * @return an exact copy of the object
+	 */
+	virtual Element* clone() const;
+
+	/**
+	 * This method should be called after at least two nodes of the prism
+	 * element are collapsed. As a consequence of the node collapsing an edge
+	 * of the prism will be collapsed. If one of the edges 3, 4 or 5 (see
+	 * sketch @ref PrismNodeAndEdgeNumbering) is collapsed we obtain a
+	 * pyramid. In this case the method will create the appropriate
+	 * object of class Pyramid.
+	 * @return a pyramid object or NULL
+	 */
+	virtual Element* reviseElement() const;
+
+protected:
+	/// Calculates the volume of a prism by subdividing it into three tetrahedra.
+	double computeVolume();
+
+	/// Return a specific edge node.
+	inline Node* getEdgeNode(unsigned edge_id, unsigned node_id) const { return _nodes[_edge_nodes[edge_id][node_id]]; };
+
+	/// Returns the ID of a face given an array of nodes.
+	unsigned identifyFace(Node* nodes[3]) const;
+
+	static const unsigned _face_nodes[5][4];
+	static const unsigned _edge_nodes[9][2];
+	static const unsigned _n_face_nodes[5];
+
+}; /* class */
+
+} /* namespace */
+
+#include "TemplatePrism.hpp"
+
+#endif /* TEMPLATEPRISM_H_ */
+
diff --git a/MeshLib/Elements/Prism.cpp b/MeshLib/Elements/TemplatePrism.hpp
similarity index 55%
rename from MeshLib/Elements/Prism.cpp
rename to MeshLib/Elements/TemplatePrism.hpp
index 76198ffdfb127d58fc124388014c101cf2b4f729..b81db0ff8a9c5485608f434bb8b53176f5580c22 100644
--- a/MeshLib/Elements/Prism.cpp
+++ b/MeshLib/Elements/TemplatePrism.hpp
@@ -5,12 +5,11 @@
  *              http://www.opengeosys.org/project/license
  *
  *
- * \file Prism.cpp
+ * \file TemplatePrism.hpp
  *
  * Created on 2012-05-02 by Karsten Rink
  */
 
-#include "Prism.h"
 #include "Node.h"
 #include "Tri.h"
 #include "Pyramid.h"
@@ -20,7 +19,8 @@
 
 namespace MeshLib {
 
-const unsigned Prism::_face_nodes[5][4] =
+template <unsigned ORDER,unsigned NNODES>
+const unsigned TemplatePrism<ORDER,NNODES>::_face_nodes[5][4] =
 {
 	{0, 2, 1, 99}, // Face 0
 	{0, 1, 4,  3}, // Face 1
@@ -29,7 +29,8 @@ const unsigned Prism::_face_nodes[5][4] =
 	{3, 4, 5, 99}  // Face 4
 };
 
-const unsigned Prism::_edge_nodes[9][2] =
+template <unsigned ORDER,unsigned NNODES>
+const unsigned TemplatePrism<ORDER,NNODES>::_edge_nodes[9][2] =
 {
 	{0, 1}, // Edge 0
 	{1, 2}, // Edge 1
@@ -42,10 +43,11 @@ const unsigned Prism::_edge_nodes[9][2] =
 	{3, 5}  // Edge 8
 };
 
-const unsigned Prism::_n_face_nodes[5] = { 3, 4, 4, 4, 3 };
+template <unsigned ORDER,unsigned NNODES>
+const unsigned TemplatePrism<ORDER,NNODES>::_n_face_nodes[5] = { 3, 4, 4, 4, 3 };
 
-
-Prism::Prism(Node* nodes[6], unsigned value)
+template <unsigned ORDER,unsigned NNODES>
+TemplatePrism<ORDER,NNODES>::TemplatePrism(Node* nodes[6], unsigned value)
 	: Cell(value)
 {
 	_nodes = nodes;
@@ -55,46 +57,36 @@ Prism::Prism(Node* nodes[6], unsigned value)
 	this->_volume = this->computeVolume();
 }
 
-Prism::Prism(Node* n0, Node* n1, Node* n2, Node* n3, Node* n4, Node* n5, unsigned value)
-	: Cell(value)
-{
-	_nodes = new Node*[6];
-	_nodes[0] = n0;
-	_nodes[1] = n1;
-	_nodes[2] = n2;
-	_nodes[3] = n3;
-	_nodes[4] = n4;
-	_nodes[5] = n5;
-	_neighbors = new Element*[5];
-	for (unsigned i=0; i<5; i++)
-		_neighbors[i] = NULL;
-	this->_volume = this->computeVolume();
-}
-
-Prism::Prism(const Prism &prism)
+template <unsigned ORDER,unsigned NNODES>
+TemplatePrism<ORDER,NNODES>::TemplatePrism(const TemplatePrism<ORDER,NNODES> &prism)
 	: Cell(prism.getValue())
 {
-	_nodes = new Node*[6];
-	for (unsigned i=0; i<6; i++)
+	_nodes = new Node*[NNODES];
+	for (unsigned i=0; i<NNODES; i++)
 		_nodes[i] = prism._nodes[i];
+
 	_neighbors = new Element*[5];
 	for (unsigned i=0; i<5; i++)
 		_neighbors[i] = prism._neighbors[i];
+
 	_volume = prism.getVolume();
 }
 
-Prism::~Prism()
+template <unsigned ORDER,unsigned NNODES>
+TemplatePrism<ORDER,NNODES>::~TemplatePrism()
 {
 }
 
-double Prism::computeVolume()
+template <unsigned ORDER,unsigned NNODES>
+double TemplatePrism<ORDER,NNODES>::computeVolume()
 {
 	return MathLib::calcTetrahedronVolume(_nodes[0]->getCoords(), _nodes[1]->getCoords(), _nodes[2]->getCoords(), _nodes[3]->getCoords())
 		 + MathLib::calcTetrahedronVolume(_nodes[1]->getCoords(), _nodes[4]->getCoords(), _nodes[2]->getCoords(), _nodes[3]->getCoords())
 		 + MathLib::calcTetrahedronVolume(_nodes[2]->getCoords(), _nodes[4]->getCoords(), _nodes[5]->getCoords(), _nodes[3]->getCoords());
 }
 
-const Element* Prism::getFace(unsigned i) const
+template <unsigned ORDER,unsigned NNODES>
+const Element* TemplatePrism<ORDER,NNODES>::getFace(unsigned i) const
 {
 	if (i<this->getNFaces())
 	{
@@ -112,7 +104,8 @@ const Element* Prism::getFace(unsigned i) const
 	return NULL;
 }
 
-unsigned Prism::getNFaceNodes(unsigned i) const
+template <unsigned ORDER,unsigned NNODES>
+unsigned TemplatePrism<ORDER,NNODES>::getNFaceNodes(unsigned i) const
 {
 	if (i<5)
 		return _n_face_nodes[i];
@@ -120,7 +113,8 @@ unsigned Prism::getNFaceNodes(unsigned i) const
 	return 0;
 }
 
-bool Prism::isEdge(unsigned idx1, unsigned idx2) const
+template <unsigned ORDER,unsigned NNODES>
+bool TemplatePrism<ORDER,NNODES>::isEdge(unsigned idx1, unsigned idx2) const
 {
 	for (unsigned i(0); i<9; i++)
 	{
@@ -130,12 +124,14 @@ bool Prism::isEdge(unsigned idx1, unsigned idx2) const
 	return false;
 }
 
-Element* Prism::clone() const
+template <unsigned ORDER,unsigned NNODES>
+Element* TemplatePrism<ORDER,NNODES>::clone() const
 {
-	return new Prism(*this);
+	return new TemplatePrism<ORDER,NNODES>(*this);
 }
 
-unsigned Prism::identifyFace(Node* nodes[3]) const
+template <unsigned ORDER,unsigned NNODES>
+unsigned TemplatePrism<ORDER,NNODES>::identifyFace(Node* nodes[3]) const
 {
 	for (unsigned i=0; i<5; i++)
 	{
@@ -150,19 +146,38 @@ unsigned Prism::identifyFace(Node* nodes[3]) const
 	return std::numeric_limits<unsigned>::max();
 }
 
-Element* Prism::reviseElement() const
+template <unsigned ORDER,unsigned NNODES>
+Element* TemplatePrism<ORDER,NNODES>::reviseElement() const
 {
 	// try to create Pyramid
 	if (_nodes[_edge_nodes[3][0]] == _nodes[_edge_nodes[3][1]]) {
-		return new Pyramid(_nodes[1], _nodes[4], _nodes[5], _nodes[2], _nodes[0], _value);
+		Node** pyramid_nodes(new Node*[5]);
+		pyramid_nodes[0] = _nodes[1];
+		pyramid_nodes[1] = _nodes[4];
+		pyramid_nodes[2] = _nodes[5];
+		pyramid_nodes[3] = _nodes[2];
+		pyramid_nodes[4] = _nodes[0];
+		return new Pyramid(pyramid_nodes, _value);
 	}
 
 	if (_nodes[_edge_nodes[4][0]] == _nodes[_edge_nodes[4][1]]) {
-		return new Pyramid(_nodes[0], _nodes[2], _nodes[5], _nodes[3], _nodes[1], _value);
+		Node** pyramid_nodes(new Node*[5]);
+		pyramid_nodes[0] = _nodes[0];
+		pyramid_nodes[1] = _nodes[2];
+		pyramid_nodes[2] = _nodes[5];
+		pyramid_nodes[3] = _nodes[3];
+		pyramid_nodes[4] = _nodes[1];
+		return new Pyramid(pyramid_nodes, _value);
 	}
 
 	if (_nodes[_edge_nodes[5][0]] == _nodes[_edge_nodes[5][1]]) {
-		return new Pyramid(_nodes[0], _nodes[1], _nodes[4], _nodes[3], _nodes[2], _value);
+		Node** pyramid_nodes(new Node*[5]);
+		pyramid_nodes[0] = _nodes[0];
+		pyramid_nodes[1] = _nodes[1];
+		pyramid_nodes[2] = _nodes[4];
+		pyramid_nodes[3] = _nodes[3];
+		pyramid_nodes[4] = _nodes[2];
+		return new Pyramid(pyramid_nodes, _value);
 	}
 
 	return NULL;
diff --git a/MeshLib/Elements/TemplateQuad.h b/MeshLib/Elements/TemplateQuad.h
new file mode 100644
index 0000000000000000000000000000000000000000..2dcae2b4612b85d100ceceef113a70af020265d4
--- /dev/null
+++ b/MeshLib/Elements/TemplateQuad.h
@@ -0,0 +1,115 @@
+/**
+ * Copyright (c) 2012, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ *
+ * \file Quad.h
+ *
+ * Created on 2012-05-02 by Karsten Rink
+ */
+
+#ifndef TEMPLATEQUAD_H_
+#define TEMPLATEQUAD_H_
+
+#include "Face.h"
+
+namespace MeshLib {
+
+/**
+ * This class represents a 2d quadliteral element. The following sketch shows the node and edge numbering.
+ * @anchor QuadNodeAndEdgeNumbering
+ * @code
+ *              2
+ *        3-----------2
+ *        |           |
+ *        |           |
+ *       3|           |1
+ *        |           |
+ *        |           |
+ *        0-----------1
+ *              0
+ * @endcode
+ */
+template <unsigned ORDER, unsigned NNODES>
+class TemplateQuad : public Face
+{
+public:
+	/// Constructor with an array of mesh nodes.
+	TemplateQuad(Node* nodes[NNODES], unsigned value = 0);
+
+	/// Copy constructor
+	TemplateQuad(const TemplateQuad &quad);
+
+	/// Destructor
+	virtual ~TemplateQuad();
+
+	/// Get the number of edges for this element.
+	unsigned getNEdges() const { return 4; };
+
+	/// Get the number of neighbors for this element.
+	unsigned getNNeighbors() const { return 4; };
+
+	/// Get the number of nodes for this element.
+	virtual unsigned getNNodes(unsigned order = 1) const
+	{
+		return order == ORDER ? NNODES : 4;
+	}
+
+	/**
+	 * Method returns the type of the element. In this case QUAD will be returned.
+	 * @return MshElemType::QUAD
+	 */
+	virtual MshElemType::type getType() const { return MshElemType::QUAD; }
+
+	/// Returns true if these two indeces form an edge and false otherwise
+	bool isEdge(unsigned i, unsigned j) const;
+
+	/**
+	 * Method clone is inherited from class Element. It makes a deep copy of the TemplateQuad instance.
+	 * @return an exact copy of the object
+	 */
+	virtual Element* clone() const;
+
+	/**
+	 * This method should be called after at least two nodes of the quad
+	 * element are collapsed. As a consequence of the node collapsing an edge
+	 * of the quad will be collapsed. If one of the edges (see
+	 * sketch @ref PyramidNodeAndEdgeNumbering) is collapsed we obtain a
+	 * triangle. In this case the method will create the appropriate
+	 * object of class Tri.
+	 * @return a Tri object or NULL
+	 */
+	virtual Element* reviseElement() const;
+
+protected:
+	/// Calculates the area of a convex quadliteral by dividing it into two triangles.
+	double computeVolume();
+
+protected:
+	/// Return a specific edge node.
+	inline Node* getEdgeNode(unsigned edge_id, unsigned node_id) const { return _nodes[_edge_nodes[edge_id][node_id]]; };
+
+	/// Returns the ID of a face given an array of nodes.
+	unsigned identifyFace(Node* nodes[3]) const;
+
+	static const unsigned _edge_nodes[4][2];
+
+}; /* class */
+
+template <unsigned ORDER, unsigned NNODES>
+const unsigned TemplateQuad<ORDER,NNODES>::_edge_nodes[4][2] =
+{
+	{0, 1}, // Edge 0
+	{1, 2}, // Edge 1
+	{2, 3}, // Edge 2
+	{0, 3}  // Edge 3
+};
+
+} /* namespace */
+
+#include "TemplateQuad.hpp"
+
+#endif /* TEMPLATEQUAD_H_ */
+
diff --git a/MeshLib/Elements/Quad.cpp b/MeshLib/Elements/TemplateQuad.hpp
similarity index 53%
rename from MeshLib/Elements/Quad.cpp
rename to MeshLib/Elements/TemplateQuad.hpp
index cc2515c3529883b85acd68c16fee80c816461922..f8523d97fa296ce4496a04b56e00a6ec2c464c88 100644
--- a/MeshLib/Elements/Quad.cpp
+++ b/MeshLib/Elements/TemplateQuad.hpp
@@ -5,12 +5,12 @@
  *              http://www.opengeosys.org/project/license
  *
  *
- * \file Quad.cpp
+ * \file Quad.hpp
  *
  * Created on 2012-05-02 by Karsten Rink
  */
 
-#include "Quad.h"
+//#include "Quad.h"
 #include "Node.h"
 #include "Tri.h"
 
@@ -18,64 +18,50 @@
 
 namespace MeshLib {
 
-
-const unsigned Quad::_edge_nodes[4][2] =
-{
-	{0, 1}, // Edge 0
-	{1, 2}, // Edge 1
-	{2, 3}, // Edge 2
-	{0, 3}  // Edge 3
-};
-
-
-Quad::Quad(Node* nodes[4], unsigned value)
+template <unsigned ORDER, unsigned NNODES>
+TemplateQuad<ORDER,NNODES>::TemplateQuad(Node* nodes[NNODES], unsigned value)
 	: Face(value)
 {
 	_nodes = nodes;
-	_neighbors = new Element*[4];
-	for (unsigned i=0; i<4; i++)
-		_neighbors[i] = NULL;
-	this->_area = this->computeVolume();
-}
 
-Quad::Quad(Node* n0, Node* n1, Node* n2, Node* n3, unsigned value)
-	: Face(value)
-{
-	_nodes = new Node*[4];
-	_nodes[0] = n0;
-	_nodes[1] = n1;
-	_nodes[2] = n2;
-	_nodes[3] = n3;
 	_neighbors = new Element*[4];
 	for (unsigned i=0; i<4; i++)
 		_neighbors[i] = NULL;
+
 	this->_area = this->computeVolume();
 }
 
-Quad::Quad(const Quad &quad)
+template <unsigned ORDER, unsigned NNODES>
+TemplateQuad<ORDER,NNODES>::TemplateQuad(const TemplateQuad<ORDER,NNODES> &quad)
 	: Face(quad.getValue())
 {
-	_nodes = new Node*[4];
-	_neighbors = new Element*[4];
-	for (unsigned i=0; i<4; i++)
-	{
+	_nodes = new Node*[NNODES];
+	for (unsigned i=0; i<NNODES; i++) {
 		_nodes[i] = quad._nodes[i];
+	}
+
+	_neighbors = new Element*[4];
+	for (unsigned i=0; i<4; i++) {
 		_neighbors[i] = quad._neighbors[i];
 	}
+
 	_area = quad.getArea();
 }
 
-Quad::~Quad()
+template <unsigned ORDER, unsigned NNODES>
+TemplateQuad<ORDER,NNODES>::~TemplateQuad()
 {
 }
 
-double Quad::computeVolume()
+template <unsigned ORDER, unsigned NNODES>
+double TemplateQuad<ORDER,NNODES>::computeVolume()
 {
 	return MathLib::calcTriangleArea(_nodes[0]->getCoords(), _nodes[1]->getCoords(), _nodes[2]->getCoords())
          + MathLib::calcTriangleArea(_nodes[2]->getCoords(), _nodes[3]->getCoords(), _nodes[0]->getCoords());
 }
 
-bool Quad::isEdge(unsigned idx1, unsigned idx2) const
+template <unsigned ORDER, unsigned NNODES>
+bool TemplateQuad<ORDER,NNODES>::isEdge(unsigned idx1, unsigned idx2) const
 {
 	for (unsigned i(0); i<4; i++)
 	{
@@ -85,12 +71,14 @@ bool Quad::isEdge(unsigned idx1, unsigned idx2) const
 	return false;
 }
 
-Element* Quad::clone() const
+template <unsigned ORDER, unsigned NNODES>
+Element* TemplateQuad<ORDER,NNODES>::clone() const
 {
-	return new Quad(*this);
+	return new TemplateQuad(*this);
 }
 
-unsigned Quad::identifyFace(Node* nodes[3]) const
+template <unsigned ORDER, unsigned NNODES>
+unsigned TemplateQuad<ORDER,NNODES>::identifyFace(Node* nodes[3]) const
 {
 	for (unsigned i=0; i<4; i++)
 	{
@@ -104,14 +92,24 @@ unsigned Quad::identifyFace(Node* nodes[3]) const
 	}
 	return std::numeric_limits<unsigned>::max();
 }
-Element* Quad::reviseElement() const
+
+template <unsigned ORDER, unsigned NNODES>
+Element* TemplateQuad<ORDER,NNODES>::reviseElement() const
 {
 	if (_nodes[0] == _nodes[1] || _nodes[1] == _nodes[2]) {
-		return new Tri(_nodes[0], _nodes[2], _nodes[3], _value);
+		MeshLib::Node** tri_nodes(new MeshLib::Node*[3]);
+		tri_nodes[0] = _nodes[0];
+		tri_nodes[1] = _nodes[2];
+		tri_nodes[2] = _nodes[3];
+		return new Tri(tri_nodes, _value);
 	}
 
 	if (_nodes[2] == _nodes[3] || _nodes[3] == _nodes[0]) {
-		return new Tri(_nodes[0], _nodes[1], _nodes[2], _value);
+		MeshLib::Node** tri_nodes(new MeshLib::Node*[3]);
+		tri_nodes[0] = _nodes[0];
+		tri_nodes[1] = _nodes[1];
+		tri_nodes[2] = _nodes[2];
+		return new Tri(tri_nodes, _value);
 	}
 
 	// this should not happen
diff --git a/MeshLib/Elements/TemplateTet.h b/MeshLib/Elements/TemplateTet.h
new file mode 100644
index 0000000000000000000000000000000000000000..2c2868cd512f07cf428b30a265184ae0b21ae11c
--- /dev/null
+++ b/MeshLib/Elements/TemplateTet.h
@@ -0,0 +1,126 @@
+/**
+ * Copyright (c) 2012, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ *
+ * \file TemplateTet.h
+ *
+ * Created on 2012-05-02 by Karsten Rink
+ */
+
+#ifndef TEMPLATETET_H_
+#define TEMPLATETET_H_
+
+#include "Cell.h"
+
+namespace MeshLib {
+
+/**
+ * This class represents a 3d tetrahedron element. The following sketch shows the node and edge numbering.
+ * @anchor TetrahedronNodeAndEdgeNumbering
+ * @code
+ *          3
+ *         /|\
+ *        / | \
+ *      3/  |  \5
+ *      /   |4  \
+ *     /    |    \
+ *    0.....|.....2
+ *     \    |  2 /
+ *      \   |   /
+ *      0\  |  /1
+ *        \ | /
+ *         \|/
+ *          1
+ *
+ * @endcode
+ */
+template <unsigned ORDER, unsigned NNODES>
+class TemplateTet : public Cell
+{
+public:
+	/// Constructor with an array of mesh nodes.
+	TemplateTet(Node* nodes[NNODES], unsigned value = 0);
+
+	/// Copy constructor
+	TemplateTet(const TemplateTet<ORDER,NNODES> &tet);
+
+	/// Destructor
+	virtual ~TemplateTet();
+
+	/// Returns the face i of the element.
+	const Element* getFace(unsigned i) const;
+
+	/// Get the number of edges for this element.
+	unsigned getNEdges() const { return 6; };
+
+	/// Get the number of nodes for face i.
+	unsigned getNFaceNodes(unsigned i) const { (void)i; return 3; };
+
+	/// Get the number of faces for this element.
+	unsigned getNFaces() const { return 4; };
+
+	/// Get the number of neighbors for this element.
+	unsigned getNNeighbors() const { return 4; };
+
+	/// Get the number of nodes for this element.
+	virtual unsigned getNNodes(unsigned order = 1) const
+	{
+		return order == ORDER ? NNODES : 4;
+	};
+
+	/**
+	 * Method returns the type of the element. In this case TETRAHEDRON will be returned.
+	 * @return MshElemType::TETRAHEDRON
+	 */
+	virtual MshElemType::type getType() const { return MshElemType::TETRAHEDRON; }
+
+	/// Returns true if these two indeces form an edge and false otherwise
+	bool isEdge(unsigned i, unsigned j) const;
+
+	/**
+	 * Method clone is inherited from class Element. It makes a deep copy of the TemplateTet instance.
+	 * @return an exact copy of the object
+	 */
+	virtual Element* clone() const;
+
+	/**
+	 * This method should be called after at least two nodes of the tetrahedron
+	 * element are collapsed. As a consequence of the node collapsing an edge
+	 * of the tetrahedron will be collapsed. If one edge is collapsed we obtain
+	 * a triangle.
+	 * @return a Triangle object or NULL
+	 */
+	virtual Element* reviseElement() const;
+
+protected:
+	/// Calculates the volume of a tetrahedron via the determinant of the matrix given by its four points.
+	double computeVolume();
+
+	/**
+	 * Return a specific edge node, see @ref TetrahedronNodeAndEdgeNumbering for numbering
+	 * @param edge_id the id/number of the edge, have to be an integer value in the interval [0,6]
+	 * @param node_id the id of the node within the edge (either 0 or 1)
+	 * @return a pointer to the internal Node
+	 */
+	inline Node* getEdgeNode(unsigned edge_id, unsigned node_id) const {
+		return _nodes[_edge_nodes[edge_id][node_id]];
+	}
+
+	/// Returns the ID of a face given an array of nodes.
+	unsigned identifyFace(Node* nodes[3]) const;
+
+
+	static const unsigned _face_nodes[4][3];
+	static const unsigned _edge_nodes[6][2];
+
+}; /* class */
+
+} /* namespace */
+
+#include "TemplateTet.hpp"
+
+#endif /* TEMPLATETET_H_ */
+
diff --git a/MeshLib/Elements/Tet.cpp b/MeshLib/Elements/TemplateTet.hpp
similarity index 54%
rename from MeshLib/Elements/Tet.cpp
rename to MeshLib/Elements/TemplateTet.hpp
index c34f03c58ed7a94942a744499ec81300d40e9c50..17e37fa42d571cae1be23c52f9f29990735d6429 100644
--- a/MeshLib/Elements/Tet.cpp
+++ b/MeshLib/Elements/TemplateTet.hpp
@@ -5,12 +5,11 @@
  *              http://www.opengeosys.org/project/license
  *
  *
- * \file Tet.cpp
+ * \file TemplateTet.hpp
  *
  * Created on 2012-05-02 by Karsten Rink
  */
 
-#include "Tet.h"
 #include "Node.h"
 #include "Tri.h"
 
@@ -18,8 +17,8 @@
 
 namespace MeshLib {
 
-
-const unsigned Tet::_face_nodes[4][3] =
+template <unsigned ORDER, unsigned NNODES>
+const unsigned TemplateTet<ORDER,NNODES>::_face_nodes[4][3] =
 {
 	{0, 2, 1}, // Face 0
 	{0, 1, 3}, // Face 1
@@ -27,7 +26,8 @@ const unsigned Tet::_face_nodes[4][3] =
 	{2, 0, 3}  // Face 3
 };
 
-const unsigned Tet::_edge_nodes[6][2] =
+template <unsigned ORDER, unsigned NNODES>
+const unsigned TemplateTet<ORDER, NNODES>::_edge_nodes[6][2] =
 {
 	{0, 1}, // Edge 0
 	{1, 2}, // Edge 1
@@ -37,61 +37,50 @@ const unsigned Tet::_edge_nodes[6][2] =
 	{2, 3}  // Edge 5
 };
 
-Tet::Tet(Node* nodes[4], unsigned value)
+template <unsigned ORDER, unsigned NNODES>
+TemplateTet<ORDER, NNODES>::TemplateTet(Node* nodes[4], unsigned value)
 	: Cell(value)
 {
 	_nodes = nodes;
-	_neighbors = new Element*[4];
-	for (unsigned i=0; i<4; i++)
-		_neighbors[i] = NULL;
-	this->_volume = this->computeVolume();
-}
 
-Tet::Tet(Node* n0, Node* n1, Node* n2, Node* n3, unsigned value)
-	: Cell(value)
-{
-	_nodes = new Node*[4];
-	_nodes[0] = n0;
-	_nodes[1] = n1;
-	_nodes[2] = n2;
-	_nodes[3] = n3;
 	_neighbors = new Element*[4];
 	for (unsigned i=0; i<4; i++)
 		_neighbors[i] = NULL;
-	this->_volume = this->computeVolume();
-}
 
-Tet::Tet(unsigned value)
-	: Cell(value)
-{
-	_neighbors = new Element*[4];
-	for (unsigned i=0; i<4; i++)
-		_neighbors[i] = NULL;
+	this->_volume = this->computeVolume();
 }
 
-Tet::Tet(const Tet &tet)
+template <unsigned ORDER, unsigned NNODES>
+TemplateTet<ORDER, NNODES>::TemplateTet(const TemplateTet<ORDER, NNODES> &tet)
 	: Cell(tet.getValue())
 {
-	_nodes = new Node*[4];
+	_nodes = new Node*[NNODES];
+	for (unsigned i=0; i<NNODES; i++) {
+		_nodes[i] = tet._nodes[i];
+	}
+
 	_neighbors = new Element*[4];
 	for (unsigned i=0; i<4; i++)
 	{
-		_nodes[i] = tet._nodes[i];
 		_neighbors[i] = tet._neighbors[i];
 	}
+
 	_volume = tet.getVolume();
 }
 
-Tet::~Tet()
+template <unsigned ORDER, unsigned NNODES>
+TemplateTet<ORDER, NNODES>::~TemplateTet()
 {
 }
 
-double Tet::computeVolume()
+template <unsigned ORDER, unsigned NNODES>
+double TemplateTet<ORDER, NNODES>::computeVolume()
 {
 	return MathLib::calcTetrahedronVolume(_nodes[0]->getCoords(), _nodes[1]->getCoords(), _nodes[2]->getCoords(), _nodes[3]->getCoords());
 }
 
-const Element* Tet::getFace(unsigned i) const
+template <unsigned ORDER, unsigned NNODES>
+const Element* TemplateTet<ORDER, NNODES>::getFace(unsigned i) const
 {
 	if (i<this->getNFaces())
 	{
@@ -105,7 +94,8 @@ const Element* Tet::getFace(unsigned i) const
 	return NULL;
 }
 
-bool Tet::isEdge(unsigned idx1, unsigned idx2) const
+template <unsigned ORDER, unsigned NNODES>
+bool TemplateTet<ORDER, NNODES>::isEdge(unsigned idx1, unsigned idx2) const
 {
 	for (unsigned i(0); i<6; i++)
 	{
@@ -115,12 +105,14 @@ bool Tet::isEdge(unsigned idx1, unsigned idx2) const
 	return false;
 }
 
-Element* Tet::clone() const
+template <unsigned ORDER, unsigned NNODES>
+Element* TemplateTet<ORDER, NNODES>::clone() const
 {
-	return new Tet(*this);
+	return new TemplateTet<ORDER,NNODES>(*this);
 }
 
-unsigned Tet::identifyFace(Node* nodes[3]) const
+template <unsigned ORDER, unsigned NNODES>
+unsigned TemplateTet<ORDER, NNODES>::identifyFace(Node* nodes[3]) const
 {
 	for (unsigned i=0; i<4; i++)
 	{
@@ -134,18 +126,32 @@ unsigned Tet::identifyFace(Node* nodes[3]) const
 	}
 	return std::numeric_limits<unsigned>::max();
 }
-Element* Tet::reviseElement() const
+
+template <unsigned ORDER, unsigned NNODES>
+Element* TemplateTet<ORDER, NNODES>::reviseElement() const
 {
 	if (_nodes[0] == _nodes[1] || _nodes[1] == _nodes[2]) {
-		return new Tri(_nodes[0], _nodes[2], _nodes[3], _value);
+		MeshLib::Node** tri_nodes(new MeshLib::Node*[3]);
+		tri_nodes[0] = _nodes[0];
+		tri_nodes[1] = _nodes[2];
+		tri_nodes[2] = _nodes[3];
+		return new Tri(tri_nodes, _value);
 	}
 
 	if (_nodes[2] == _nodes[0]) {
-		return new Tri(_nodes[0], _nodes[1], _nodes[3], _value);
+		MeshLib::Node** tri_nodes(new MeshLib::Node*[3]);
+		tri_nodes[0] = _nodes[0];
+		tri_nodes[1] = _nodes[1];
+		tri_nodes[2] = _nodes[3];
+		return new Tri(tri_nodes, _value);
 	}
 
 	if (_nodes[0] == _nodes[3] || _nodes[1] == _nodes[3] || _nodes[2] == _nodes[3]) {
-		return new Tri(_nodes[0], _nodes[1], _nodes[2], _value);
+		MeshLib::Node** tri_nodes(new MeshLib::Node*[3]);
+		tri_nodes[0] = _nodes[0];
+		tri_nodes[1] = _nodes[1];
+		tri_nodes[2] = _nodes[2];
+		return new Tri(tri_nodes, _value);
 	}
 
 	// this should not happen
diff --git a/MeshLib/Elements/TemplateTri.h b/MeshLib/Elements/TemplateTri.h
new file mode 100644
index 0000000000000000000000000000000000000000..5fad5041b3e6e1ab81060d5490659ee92ecb8c5d
--- /dev/null
+++ b/MeshLib/Elements/TemplateTri.h
@@ -0,0 +1,128 @@
+/**
+ * Copyright (c) 2012, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ *
+ * \file TemplateTri.h
+ *
+ * Created on 2012-05-02 by Karsten Rink
+ */
+
+#ifndef TEMPLATETRI_H_
+#define TEMPLATETRI_H_
+
+#include "Edge.h"
+#include "Node.h"
+#include "Face.h"
+
+#include "MathTools.h"
+
+
+namespace MeshLib {
+
+/**
+ * This class represents a 2d triangle element. The following sketch shows the node and edge numbering.
+ * @anchor TriNodeAndEdgeNumbering
+ * @code
+ *
+ *          2
+ *          o
+ *         / \
+ *        /   \
+ *      2/     \1
+ *      /       \
+ *     /         \
+ *    0-----------1
+ *          0
+ *
+ * @endcode
+ */
+template <unsigned ORDER, unsigned NNODES>
+class TemplateTri : public Face
+{
+public:
+	/// Constructor with an array of mesh nodes.
+	TemplateTri(Node* nodes[NNODES], unsigned value = 0);
+
+	/// Copy constructor
+	TemplateTri(const TemplateTri<ORDER, NNODES> &tri);
+
+	/// Destructor
+	virtual ~TemplateTri();
+
+	/// Get the number of edges for this element.
+	unsigned getNEdges() const { return 3; };
+
+	/// Get the number of neighbors for this element.
+	unsigned getNNeighbors() const { return 3; };
+
+	/// Get the number of nodes for this element.
+	virtual unsigned getNNodes(unsigned order = 1) const
+	{
+		return order == ORDER ? NNODES : 3;
+	}
+
+	/**
+	 * Method returns the type of the element. In this case TRIANGLE will be returned.
+	 * @return MshElemType::TRIANGLE
+	 */
+	virtual MshElemType::type getType() const { return MshElemType::TRIANGLE; }
+
+	/// Returns true if these two indices form an edge and false otherwise
+	bool isEdge(unsigned idx1, unsigned idx2) const;
+
+	/**
+	 * Method clone is inherited from class Element. It makes a deep copy of the TemplateTri instance.
+	 * @return an exact copy of the object
+	 */
+	virtual Element* clone() const
+	{
+		return new TemplateTri<ORDER, NNODES>(*this);
+	}
+
+
+	/**
+	 * This method should be called after at least two nodes of the triangle
+	 * element are collapsed. As a consequence of the node collapsing an edge
+	 * of the triangle will be collapsed. If one of the edges is collapsed we
+	 * obtain an edge. In this case the method will create the appropriate
+	 * object of class Edge.
+	 * @return an Edge object or NULL
+	 */
+	virtual Element* reviseElement() const;
+
+protected:
+	/// Calculates the area of the triangle by returning half of the area of the corresponding parallelogram.
+	double computeVolume()
+	{
+		return MathLib::calcTriangleArea(_nodes[0]->getCoords(), _nodes[1]->getCoords(), _nodes[2]->getCoords());
+	}
+
+protected:
+	/// Return a specific edge node.
+	inline Node* getEdgeNode(unsigned edge_id, unsigned node_id) const
+	{
+		return _nodes[_edge_nodes[edge_id][node_id]];
+	}
+
+	/// Returns the ID of a face given an array of nodes.
+	unsigned identifyFace(Node* nodes[3]) const;
+
+	static const unsigned _edge_nodes[3][2];
+}; /* class */
+
+template <unsigned ORDER, unsigned NNODES>
+const unsigned TemplateTri<ORDER, NNODES>::_edge_nodes[3][2] = {
+		{0, 1}, // Edge 0
+		{1, 2}, // Edge 1
+		{0, 2}  // Edge 2
+	};
+
+} /* namespace */
+
+#include "TemplateTri.hpp"
+
+#endif /* TEMPLATETRI_H_ */
+
diff --git a/MeshLib/Elements/TemplateTri.hpp b/MeshLib/Elements/TemplateTri.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..20bd6ee5380b48ce219d5bd6c96fb4a8450711a3
--- /dev/null
+++ b/MeshLib/Elements/TemplateTri.hpp
@@ -0,0 +1,100 @@
+/**
+ * Copyright (c) 2012, 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 TemplateTri.hpp
+ *
+ *  Created on  Sep 27, 2012 by Thomas Fischer
+ */
+
+#ifndef TEMPLATETRI_HPP_
+#define TEMPLATETRI_HPP_
+
+namespace MeshLib {
+
+template <unsigned ORDER, unsigned NNODES>
+TemplateTri<ORDER,NNODES>::TemplateTri(Node* nodes[NNODES], unsigned value) :
+	Face(value)
+{
+	_nodes = nodes;
+	_neighbors = new Element*[3];
+	for (unsigned i=0; i<3; i++)
+		_neighbors[i] = NULL;
+	this->_area = this->computeVolume();
+}
+
+template <unsigned ORDER, unsigned NNODES>
+TemplateTri<ORDER,NNODES>::TemplateTri(const TemplateTri<ORDER, NNODES> &tri) :
+	Face(tri.getValue())
+{
+	_nodes = new Node*[NNODES];
+	for (unsigned i=0; i<NNODES; i++)
+	{
+		_nodes[i] = tri._nodes[i];
+	}
+
+	_neighbors = new Element*[3];
+	for (unsigned i=0; i<3; i++) {
+		_neighbors[i] = tri._neighbors[i];
+	}
+
+	_area = tri.getArea();
+}
+
+template <unsigned ORDER, unsigned NNODES>
+TemplateTri<ORDER,NNODES>::~TemplateTri()
+{}
+
+template <unsigned ORDER, unsigned NNODES>
+bool TemplateTri<ORDER,NNODES>::isEdge(unsigned idx1, unsigned idx2) const
+{
+	for (unsigned i(0); i<3; i++)
+	{
+		if (_edge_nodes[i][0]==idx1 && _edge_nodes[i][1]==idx2) return true;
+		if (_edge_nodes[i][1]==idx1 && _edge_nodes[i][0]==idx2) return true;
+	}
+	return false;
+}
+
+template <unsigned ORDER, unsigned NNODES>
+Element* TemplateTri<ORDER,NNODES>::reviseElement() const
+{
+	// try to create an edge
+	if (_nodes[0] == _nodes[1] || _nodes[1] == _nodes[2]) {
+		Node** nodes (new Node*[2]);
+		nodes[0] = _nodes[0];
+		nodes[1] = _nodes[2];
+		return new Edge(nodes, _value);
+	}
+
+	if (_nodes[0] == _nodes[2]) {
+		Node** nodes (new Node*[2]);
+		nodes[0] = _nodes[0];
+		nodes[1] = _nodes[1];
+		return new Edge(nodes, _value);
+	}
+
+	return NULL;
+}
+
+template <unsigned ORDER, unsigned NNODES>
+unsigned TemplateTri<ORDER,NNODES>::identifyFace(Node* nodes[3]) const
+{
+	for (unsigned i=0; i<3; i++)
+	{
+		unsigned flag(0);
+		for (unsigned j=0; j<2; j++)
+			for (unsigned k=0; k<2; k++)
+				if (_nodes[_edge_nodes[i][j]] == nodes[k])
+					flag++;
+		if (flag==2)
+			return i;
+	}
+	return std::numeric_limits<unsigned>::max();
+}
+
+} // namespace MeshLib
+
+#endif /* TEMPLATETRI_HPP_ */
diff --git a/MeshLib/Elements/Tet.h b/MeshLib/Elements/Tet.h
index 4958fcea391183309a0e27388592f11e5ef27694..73e920023b211448f50222be89b657954d8baff4 100644
--- a/MeshLib/Elements/Tet.h
+++ b/MeshLib/Elements/Tet.h
@@ -1,124 +1,23 @@
 /**
- * Copyright (c) 2012, OpenGeoSys Community (http://www.opengeosys.org)
+ * Copyright (c) 2012, OpenGeoSys Community (http://www.opengeosys.net)
  *            Distributed under a Modified BSD License.
  *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
+ *              http://www.opengeosys.org/LICENSE.txt
  *
  * \file Tet.h
  *
- * Created on 2012-05-02 by Karsten Rink
+ *  Created on  Sep 27, 2012 by Thomas Fischer
  */
 
 #ifndef TET_H_
 #define TET_H_
 
-#include "Cell.h"
+#include "TemplateTet.h"
 
 namespace MeshLib {
 
-/**
- * This class represents a 3d tetrahedron element. The following sketch shows the node and edge numbering.
- * @anchor TetrahedronNodeAndEdgeNumbering
- * @code
- *          3
- *         /|\
- *        / | \
- *      3/  |  \5
- *      /   |4  \
- *     /    |    \
- *    0.....|.....2
- *     \    |  2 /
- *      \   |   /
- *      0\  |  /1
- *        \ | /
- *         \|/
- *          1
- *
- * @endcode
- */
-class Tet : public Cell
-{
-public:
-	/// Constructor with an array of mesh nodes.
-	Tet(Node* nodes[4], unsigned value = 0);
-
-	/// Constructor using single mesh nodes.
-	Tet(Node* n0, Node* n1, Node* n2, Node* n3, unsigned value = 0);
-
-	/// Copy constructor
-	Tet(const Tet &tet);
-
-	/// Destructor
-	virtual ~Tet();
-
-	/// Returns the face i of the element.
-	const Element* getFace(unsigned i) const;
-
-	/// Get the number of edges for this element.
-	unsigned getNEdges() const { return 6; };
-
-	/// Get the number of nodes for face i.
-	unsigned getNFaceNodes(unsigned i) const { (void)i; return 3; };
-
-	/// Get the number of faces for this element.
-	unsigned getNFaces() const { return 4; };
-
-	/// Get the number of neighbors for this element.
-	unsigned getNNeighbors() const { return 4; };
-
-	/// Get the number of nodes for this element.
-	virtual unsigned getNNodes() const { return 4; };
+typedef TemplateTet<1,4> Tet;
 
-	/**
-	 * Method returns the type of the element. In this case TETRAHEDRON will be returned.
-	 * @return MshElemType::TETRAHEDRON
-	 */
-	virtual MshElemType::type getType() const { return MshElemType::TETRAHEDRON; }
-
-	/// Returns true if these two indeces form an edge and false otherwise
-	bool isEdge(unsigned i, unsigned j) const;
-
-	/**
-	 * Method clone is inherited from class Element. It makes a deep copy of the Tet instance.
-	 * @return an exact copy of the object
-	 */
-	virtual Element* clone() const;
-
-	/**
-	 * This method should be called after at least two nodes of the tetrahedron
-	 * element are collapsed. As a consequence of the node collapsing an edge
-	 * of the tetrahedron will be collapsed. If one edge is collapsed we obtain
-	 * a triangle.
-	 * @return a Triangle object or NULL
-	 */
-	virtual Element* reviseElement() const;
-
-protected:
-	/// Constructor without nodes (for use of derived classes)
-	Tet(unsigned value = 0);
-
-	/// Calculates the volume of a tetrahedron via the determinant of the matrix given by its four points.
-	double computeVolume();
-
-	/**
-	 * Return a specific edge node, see @ref TetrahedronNodeAndEdgeNumbering for numbering
-	 * @param edge_id the id/number of the edge, have to be an integer value in the interval [0,6]
-	 * @param node_id the id of the node within the edge (either 0 or 1)
-	 * @return a pointer to the internal Node
-	 */
-	inline Node* getEdgeNode(unsigned edge_id, unsigned node_id) const { return _nodes[_edge_nodes[edge_id][node_id]]; };
-
-	/// Returns the ID of a face given an array of nodes.
-	unsigned identifyFace(Node* nodes[3]) const;
-
-
-	static const unsigned _face_nodes[4][3];
-	static const unsigned _edge_nodes[6][2];
-
-}; /* class */
-
-} /* namespace */
+}
 
 #endif /* TET_H_ */
-
diff --git a/MeshLib/Elements/Tri.cpp b/MeshLib/Elements/Tri.cpp
deleted file mode 100644
index 4ff9fb8561ebff2f4700ce63046338cbee4d3083..0000000000000000000000000000000000000000
--- a/MeshLib/Elements/Tri.cpp
+++ /dev/null
@@ -1,120 +0,0 @@
-/**
- * Copyright (c) 2012, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- *
- * \file Tri.cpp
- *
- * Created on 2012-05-02 by Karsten Rink
- */
-
-#include "Tri.h"
-#include "Edge.h"
-#include "Node.h"
-
-#include "MathTools.h"
-
-namespace MeshLib {
-
-
-const unsigned Tri::_edge_nodes[3][2] =
-{
-	{0, 1}, // Edge 0
-	{1, 2}, // Edge 1
-	{0, 2}  // Edge 2
-};
-
-
-Tri::Tri(Node* nodes[3], unsigned value)
-	: Face(value)
-{
-	_nodes = nodes;
-	_neighbors = new Element*[3];
-	for (unsigned i=0; i<3; i++)
-		_neighbors[i] = NULL;
-	this->_area = this->computeVolume();
-}
-
-Tri::Tri(Node* n0, Node* n1, Node* n2, unsigned value)
-	: Face(value)
-{
-	_nodes = new Node*[3];
-	_nodes[0] = n0;
-	_nodes[1] = n1;
-	_nodes[2] = n2;
-	_neighbors = new Element*[3];
-	for (unsigned i=0; i<3; i++)
-		_neighbors[i] = NULL;
-	this->_area = this->computeVolume();
-}
-
-Tri::Tri(const Tri &tri)
-	: Face(tri.getValue())
-{
-	_nodes = new Node*[3];
-	_neighbors = new Element*[3];
-	for (unsigned i=0; i<3; i++)
-	{
-		_nodes[i] = tri._nodes[i];
-		_neighbors[i] = tri._neighbors[i];
-	}
-	_area = tri.getArea();
-}
-
-Tri::~Tri()
-{
-}
-
-bool Tri::isEdge(unsigned idx1, unsigned idx2) const
-{
-	for (unsigned i(0); i<3; i++)
-	{
-		if (_edge_nodes[i][0]==idx1 && _edge_nodes[i][1]==idx2) return true;
-		if (_edge_nodes[i][1]==idx1 && _edge_nodes[i][0]==idx2) return true;
-	}
-	return false;
-}
-
-Element* Tri::clone() const
-{
-	return new Tri(*this);
-}
-
-double Tri::computeVolume()
-{
-	return MathLib::calcTriangleArea(_nodes[0]->getCoords(), _nodes[1]->getCoords(), _nodes[2]->getCoords());
-}
-
-unsigned Tri::identifyFace(Node* nodes[3]) const
-{
-	for (unsigned i=0; i<3; i++)
-	{
-		unsigned flag(0);
-		for (unsigned j=0; j<2; j++)
-			for (unsigned k=0; k<2; k++)
-				if (_nodes[_edge_nodes[i][j]] == nodes[k])
-					flag++;
-		if (flag==2)
-			return i;
-	}
-	return std::numeric_limits<unsigned>::max();
-}
-
-Element* Tri::reviseElement() const
-{
-	// try to create an edge
-	if (_nodes[0] == _nodes[1] || _nodes[1] == _nodes[2]) {
-		return new Edge(_nodes[0], _nodes[2], _value);
-	}
-
-	if (_nodes[0] == _nodes[2]) {
-		return new Edge(_nodes[0], _nodes[1], _value);
-	}
-
-	return NULL;
-}
-
-}
-
diff --git a/MeshLib/Elements/Tri.h b/MeshLib/Elements/Tri.h
index b4cc0da6d9afa249bf02345f55586e73a94c8936..d0b4be2b047dd18d1d76f983030f563934b5eb07 100644
--- a/MeshLib/Elements/Tri.h
+++ b/MeshLib/Elements/Tri.h
@@ -1,104 +1,23 @@
 /**
- * Copyright (c) 2012, OpenGeoSys Community (http://www.opengeosys.org)
+ * Copyright (c) 2012, OpenGeoSys Community (http://www.opengeosys.net)
  *            Distributed under a Modified BSD License.
  *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
+ *              http://www.opengeosys.org/LICENSE.txt
  *
  * \file Tri.h
  *
- * Created on 2012-05-02 by Karsten Rink
+ *  Created on  Sep 27, 2012 by Thomas Fischer
  */
 
 #ifndef TRI_H_
 #define TRI_H_
 
-#include "Face.h"
+#include "TemplateTri.h"
 
 namespace MeshLib {
 
-/**
- * This class represents a 2d triangle element. The following sketch shows the node and edge numbering.
- * @anchor TriNodeAndEdgeNumbering
- * @code
- *
- *          2
- *          o
- *         / \
- *        /   \
- *      2/     \1
- *      /       \
- *     /         \
- *    0-----------1
- *          0
- *
- * @endcode
- */
-class Tri : public Face
-{
-public:
-	/// Constructor with an array of mesh nodes.
-	Tri(Node* nodes[3], unsigned value = 0);
-
-	/// Constructor using single mesh nodes.
-	Tri(Node* n0, Node* n1, Node* n2, unsigned value = 0);
-
-	/// Copy constructor
-	Tri(const Tri &tri);
-
-	/// Destructor
-	virtual ~Tri();
-
-	/// Get the number of edges for this element.
-	unsigned getNEdges() const { return 3; };
-
-	/// Get the number of neighbors for this element.
-	unsigned getNNeighbors() const { return 3; };
-
-	/// Get the number of nodes for this element.
-	virtual unsigned getNNodes() const { return 3; };
+typedef TemplateTri<1,3> Tri;
 
-	/**
-	 * Method returns the type of the element. In this case TRIANGLE will be returned.
-	 * @return MshElemType::TRIANGLE
-	 */
-	virtual MshElemType::type getType() const { return MshElemType::TRIANGLE; }
-
-	/// Returns true if these two indeces form an edge and false otherwise
-	bool isEdge(unsigned i, unsigned j) const;
-
-	/**
-	 * Method clone is inherited from class Element. It makes a deep copy of the Tri instance.
-	 * @return an exact copy of the object
-	 */
-	virtual Element* clone() const;
-
-	/**
-	 * This method should be called after at least two nodes of the triangle
-	 * element are collapsed. As a consequence of the node collapsing an edge
-	 * of the triangle will be collapsed. If one of the edges is collapsed we
-	 * obtain an edge. In this case the method will create the appropriate
-	 * object of class Edge.
-	 * @return an Edge object or NULL
-	 */
-	virtual Element* reviseElement() const;
-
-protected:
-	/// Calculates the area of the triangle by returning half of the area of the corresponding parallelogram.
-	double computeVolume();
-
-protected:
-	/// Return a specific edge node.
-	inline Node* getEdgeNode(unsigned edge_id, unsigned node_id) const { return _nodes[_edge_nodes[edge_id][node_id]]; };
-
-	/// Returns the ID of a face given an array of nodes.
-	unsigned identifyFace(Node* nodes[3]) const;
-
-	static const unsigned _edge_nodes[3][2];
-
-}; /* class */
-
-} /* namespace */
+}
 
 #endif /* TRI_H_ */
-
diff --git a/MeshLib/MshEditor.cpp b/MeshLib/MshEditor.cpp
index 9e91048d7d5cc4750231e3bc75545757c427558e..f83467b645d27c062eb918e87989dbabbc937be8 100644
--- a/MeshLib/MshEditor.cpp
+++ b/MeshLib/MshEditor.cpp
@@ -134,15 +134,17 @@ MeshLib::Mesh* MshEditor::getMeshSurface(const MeshLib::Mesh &mesh, const double
 	for (unsigned i=0; i<nNewElements; ++i)
 	{
 		MeshLib::Element* elem (sfc_elements[i]);
-		if (elem->getType() == MshElemType::TRIANGLE)
-			new_elements[i] = new MeshLib::Tri(sfc_nodes[node_id_map[elem->getNode(0)->getID()]],
-				                               sfc_nodes[node_id_map[elem->getNode(1)->getID()]],
-											   sfc_nodes[node_id_map[elem->getNode(2)->getID()]]);
-		else
-			new_elements[i] = new MeshLib::Quad(sfc_nodes[node_id_map[elem->getNode(0)->getID()]],
-				                                sfc_nodes[node_id_map[elem->getNode(1)->getID()]],
-												sfc_nodes[node_id_map[elem->getNode(2)->getID()]],
-												sfc_nodes[node_id_map[elem->getNode(3)->getID()]]);
+		if (elem->getType() == MshElemType::TRIANGLE) {
+			MeshLib::Node** tri_nodes (new MeshLib::Node*[3]);
+			for (unsigned k(0); k<3; k++)
+				tri_nodes[k] = sfc_nodes[node_id_map[elem->getNode(k)->getID()]];
+			new_elements[i] = new MeshLib::Tri(tri_nodes);
+		} else {
+			MeshLib::Node** quad_nodes (new MeshLib::Node*[4]);
+			for (unsigned k(0); k<3; k++)
+				quad_nodes[k] = sfc_nodes[node_id_map[elem->getNode(k)->getID()]];
+			new_elements[i] = new MeshLib::Quad(quad_nodes);
+		}
 		delete sfc_elements[i];
 	}
 
@@ -161,7 +163,7 @@ void MshEditor::get2DSurfaceElements(const std::vector<MeshLib::Element*> &all_e
 	{
 		for (unsigned i=0; i<nElements; ++i)
 		{
-			if (complete_surface) 
+			if (complete_surface)
 				sfc_elements.push_back(all_elements[i]);
 			else
 			{
@@ -244,7 +246,7 @@ std::vector<GeoLib::PointWithID*> MshEditor::getSurfaceNodes(const MeshLib::Mesh
 	std::vector<MeshLib::Node*> sfc_nodes;
 	std::vector<unsigned> node_id_map(mesh.getNNodes());
 	get2DSurfaceNodes(all_nodes, sfc_nodes, sfc_elements, node_id_map);
-	
+
 	const unsigned nElements (sfc_elements.size());
 	for (unsigned i=0; i<nElements; ++i)
 		delete sfc_elements[i];
diff --git a/Utils/SimpleMeshCreation/generateStructuredQuadMesh.cpp b/Utils/SimpleMeshCreation/generateStructuredQuadMesh.cpp
index 0c52a87f0b03248c4269fff3d00588f907ce4e18..9760e3df6d55c292f0ec9b0345ab3e7edf08c02c 100644
--- a/Utils/SimpleMeshCreation/generateStructuredQuadMesh.cpp
+++ b/Utils/SimpleMeshCreation/generateStructuredQuadMesh.cpp
@@ -2,7 +2,7 @@
  * Copyright (c) 2012, OpenGeoSys Community (http://www.opengeosys.net)
  *            Distributed under a Modified BSD License.
  *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.net/LICENSE.txt
+ *              http://www.opengeosys.org/LICENSE.txt
  *
  * \file generateStructuredQuadMesh.cpp
  *
diff --git a/scripts/cmake/Coverage.cmake b/scripts/cmake/Coverage.cmake
new file mode 100644
index 0000000000000000000000000000000000000000..4c5e17062ef06d76e53d1eb9132fb50617ac24f5
--- /dev/null
+++ b/scripts/cmake/Coverage.cmake
@@ -0,0 +1,17 @@
+INCLUDE(CodeCoverage)
+
+SET(COVERAGE_EXCLUDES
+	'/usr/*'
+	'${CMAKE_BINARY_DIR}/*'
+	'${CMAKE_SOURCE_DIR}/tests/*'
+	'${CMAKE_SOURCE_DIR}/BaseLib/zlib/*'
+	'${CMAKE_SOURCE_DIR}/BaseLib/logog/*'
+	'${CMAKE_SOURCE_DIR}/BaseLib/RapidXML/*'
+)
+
+IF(JENKINS_URL)
+	SETUP_TARGET_FOR_COVERAGE_COBERTURA(ctest_coverage ctest "ctest_coverage_results" "-j;${PROCESSOR_COUNT}")
+ELSE()
+	SETUP_TARGET_FOR_COVERAGE(ctest_coverage ctest "ctest_coverage_results" "-j;${PROCESSOR_COUNT}")
+	SETUP_TARGET_FOR_COVERAGE(ogs-gui_coverage ${EXECUTABLE_OUTPUT_PATH}/ogs-gui "ogs-gui_coverage_results")
+ENDIF()
\ No newline at end of file
diff --git a/scripts/cmake/Find.cmake b/scripts/cmake/Find.cmake
index 2277bcc4356ed75e0cc2fb57a25b88d85e9b2432..75ae2794957a4360c67df8acd04ff4464e062c73 100644
--- a/scripts/cmake/Find.cmake
+++ b/scripts/cmake/Find.cmake
@@ -25,6 +25,8 @@ FIND_PROGRAM(GPROF_PATH gprof DOC "GNU profiler gprof")
 
 FIND_PACKAGE(cppcheck)
 
+FIND_PACKAGE(PythonInterp)
+
 ######################
 ### Find libraries ###
 ######################
diff --git a/scripts/cmake/Test.cmake b/scripts/cmake/Test.cmake
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/scripts/cmake/cmake/CodeCoverage.cmake b/scripts/cmake/cmake/CodeCoverage.cmake
index 3529a098d06a981d1b16ec74be85c089b5f13676..0e7837fb94b7c27e91e0eaf30f5568c7c40988c8 100644
--- a/scripts/cmake/cmake/CodeCoverage.cmake
+++ b/scripts/cmake/cmake/CodeCoverage.cmake
@@ -6,8 +6,9 @@
 # 1. Copy this file into your cmake modules path
 # 2. Add the following line to your CMakeLists.txt:
 #      INCLUDE(CodeCoverage)
-# 
-# 3. Use the function SETUP_TARGET_FOR_COVERAGE to create a custom make target
+#
+# 3. SET(COVERAGE_EXCLUDES 'dir1/*' 'dir2/*')
+# 4. Use the function SETUP_TARGET_FOR_COVERAGE to create a custom make target
 #    which runs your test executable and produces a lcov code coverage report.
 #
 
@@ -15,7 +16,7 @@
 FIND_PROGRAM( GCOV_PATH gcov )
 FIND_PROGRAM( LCOV_PATH lcov )
 FIND_PROGRAM( GENHTML_PATH genhtml )
-FIND_PROGRAM( GCOVR_PATH gcovr PATHS ${CMAKE_SOURCE_DIR}/tests)
+FIND_PROGRAM( GCOVR_PATH gcovr PATHS ${CMAKE_SOURCE_DIR}/scripts/test)
 
 IF(NOT GCOV_PATH)
 	MESSAGE(FATAL_ERROR "gcov not found! Aborting...")
@@ -53,23 +54,23 @@ FUNCTION(SETUP_TARGET_FOR_COVERAGE _targetname _testrunner _outputname)
 
 	# Setup target
 	ADD_CUSTOM_TARGET(${_targetname}
-		
+
 		# Cleanup lcov
 		${LCOV_PATH} --directory . --zerocounters
-		
+
 		# Run tests
 		COMMAND ${_testrunner} ${ARGV3}
-		
+
 		# Capturing lcov counters and generating report
 		COMMAND ${LCOV_PATH} --directory . --capture --output-file ${_outputname}.info
-		COMMAND ${LCOV_PATH} --remove ${_outputname}.info 'tests/*' '/usr/*' --output-file ${_outputname}.info.cleaned
+		COMMAND ${LCOV_PATH} --remove ${_outputname}.info ${COVERAGE_EXCLUDES} --output-file ${_outputname}.info.cleaned
 		COMMAND ${GENHTML_PATH} -o ${_outputname} ${_outputname}.info.cleaned
 		COMMAND ${CMAKE_COMMAND} -E remove ${_outputname}.info ${_outputname}.info.cleaned
-		
+
 		WORKING_DIRECTORY ${CMAKE_BINARY_DIR}
 		COMMENT "Resetting code coverage counters to zero.\nProcessing code coverage counters and generating report."
 	)
-	
+
 	# Show info where to find the report
 	ADD_CUSTOM_COMMAND(TARGET ${_targetname} POST_BUILD
 		COMMAND ;
@@ -99,7 +100,7 @@ FUNCTION(SETUP_TARGET_FOR_COVERAGE_COBERTURA _targetname _testrunner _outputname
 		${_testrunner} ${ARGV3}
 
 		# Running gcovr
-		COMMAND ${GCOVR_PATH} -x -r ${CMAKE_SOURCE_DIR} -e '${CMAKE_SOURCE_DIR}/tests/'  -o ${_outputname}.xml
+		COMMAND ${GCOVR_PATH} -x -r ${CMAKE_SOURCE_DIR} -e ${COVERAGE_EXCLUDES} -o ${_outputname}.xml
 		WORKING_DIRECTORY ${CMAKE_BINARY_DIR}
 		COMMENT "Running gcovr to produce Cobertura code coverage report."
 	)
diff --git a/scripts/test/gcovr b/scripts/test/gcovr
new file mode 100755
index 0000000000000000000000000000000000000000..2c1f8800efd0a0f06db2ea150b5360e5f4ef311c
--- /dev/null
+++ b/scripts/test/gcovr
@@ -0,0 +1,986 @@
+#! /usr/bin/env python
+#
+# A report generator for gcov 3.4
+#
+# This routine generates a format that is similar to the format generated
+# by the Python coverage.py module.  This code is similar to the
+# data processing performed by lcov's geninfo command.  However, we
+# don't worry about parsing the *.gcna files, and backwards compatibility for
+# older versions of gcov is not supported.
+#
+# Outstanding issues
+#   - verify that gcov 3.4 or newer is being used
+#   - verify support for symbolic links
+#
+# gcovr is a FAST project.  For documentation, bug reporting, and
+# updates, see https://software.sandia.gov/trac/fast/wiki/gcovr
+#
+# _________________________________________________________________________
+#
+# FAST: Utilities for Agile Software Development
+# Copyright (c) 2008 Sandia Corporation.
+# This software is distributed under the BSD License.
+# Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
+# the U.S. Government retains certain rights in this software.
+# For more information, see the FAST README.txt file.
+#
+# $Revision: 2774 $
+# $Date: 2012-04-13 17:26:55 -0700 (Fri, 13 Apr 2012) $
+# _________________________________________________________________________
+#
+
+import copy
+import glob
+import os
+import re
+import subprocess
+import sys
+import time
+import xml.dom.minidom
+
+from optparse import OptionParser
+from string import Template
+from os.path import normpath
+
+__version__ = "2.5-prerelease"
+src_revision = "$Revision: 2774 $"
+gcov_cmd = "gcov"
+
+output_re = re.compile("[Cc]reating [`'](.*)'$")
+source_re = re.compile("cannot open (source|graph) file")
+
+def version_str():
+    ans = __version__
+    m = re.match('\$Revision:\s*(\S+)\s*\$', src_revision)
+    if m:
+        ans = ans + " (r%s)" % (m.group(1))
+    return ans
+
+#
+# Container object for coverage statistics
+#
+class CoverageData(object):
+
+    def __init__(self, fname, uncovered, covered, branches, noncode):
+        self.fname=fname
+        # Shallow copies are cheap & "safe" because the caller will
+        # throw away their copies of covered & uncovered after calling
+        # us exactly *once*
+        self.uncovered = copy.copy(uncovered)
+        self.covered   = copy.copy(covered)
+        self.noncode   = copy.copy(noncode)
+        # But, a deep copy is required here
+        self.all_lines = copy.deepcopy(uncovered)
+        self.all_lines.update(covered.keys())
+        self.branches = copy.deepcopy(branches)
+
+    def update(self, uncovered, covered, branches, noncode):
+        self.all_lines.update(uncovered)
+        self.all_lines.update(covered.keys())
+        self.uncovered.update(uncovered)
+        self.noncode.intersection_update(noncode)
+        for k in covered.keys():
+            self.covered[k] = self.covered.get(k,0) + covered[k]
+        for k in branches.keys():
+            for b in branches[k]:
+                d = self.branches.setdefault(k, {})
+                d[b] = d.get(b, 0) + branches[k][b]
+        self.uncovered.difference_update(self.covered.keys())
+
+    def uncovered_str(self):
+        if options.show_branch:
+            # Don't do any aggregation on branch results
+            tmp = []
+            for line in self.branches.keys():
+                for branch in self.branches[line]:
+                    if self.branches[line][branch] == 0:
+                        tmp.append(line)
+                        break
+
+            tmp.sort()
+            return ",".join([str(x) for x in tmp]) or ""
+
+        tmp = list(self.uncovered)
+        if len(tmp) == 0:
+            return ""
+
+        tmp.sort()
+        first = None
+        last = None
+        ranges=[]
+        for item in tmp:
+            if last is None:
+                first=item
+                last=item
+            elif item == (last+1):
+                last=item
+            else:
+                if len(self.noncode.intersection(range(last+1,item))) \
+                       == item - last - 1:
+                    last = item
+                    continue
+
+                if first==last:
+                    ranges.append(str(first))
+                else:
+                    ranges.append(str(first)+"-"+str(last))
+                first=item
+                last=item
+        if first==last:
+            ranges.append(str(first))
+        else:
+            ranges.append(str(first)+"-"+str(last))
+        return ",".join(ranges)
+
+    def coverage(self):
+        if ( options.show_branch ):
+            total = 0
+            cover = 0
+            for line in self.branches.keys():
+                for branch in self.branches[line].keys():
+                    total += 1
+                    cover += self.branches[line][branch] > 0 and 1 or 0
+        else:
+            total = len(self.all_lines)
+            cover = len(self.covered)
+
+        percent = total and str(int(100.0*cover/total)) or "--"
+        return (total, cover, percent)
+
+    def summary(self):
+        tmp = options.filter.sub('',self.fname)
+        if not self.fname.endswith(tmp):
+            # Do no truncation if the filter does not start matching at
+            # the beginning of the string
+            tmp = self.fname
+        tmp = tmp.ljust(40)
+        if len(tmp) > 40:
+            tmp=tmp+"\n"+" "*40
+
+        (total, cover, percent) = self.coverage()
+        return ( total, cover,
+                 tmp + str(total).rjust(8) + str(cover).rjust(8) + \
+                 percent.rjust(6) + "%   " + self.uncovered_str() )
+
+
+def resolve_symlinks(orig_path):
+    """
+    Return the normalized absolute path name with all symbolic links resolved
+    """
+    drive,tmp = os.path.splitdrive(os.path.abspath(orig_path))
+    if not drive:
+        drive = os.path.sep
+    parts = tmp.split(os.path.sep)
+    actual_path = [drive]
+    while parts:
+        actual_path.append(parts.pop(0))
+        if not os.path.islink(os.path.join(*actual_path)):
+            continue
+        actual_path[-1] = os.readlink(os.path.join(*actual_path))
+        tmp_drive, tmp_path = os.path.splitdrive(
+            resolve_symlinks(os.path.join(*actual_path)) )
+        if tmp_drive:
+            drive = tmp_drive
+        actual_path = [drive] + tmp_path.split(os.path.sep)
+    return os.path.join(*actual_path)
+
+
+def path_startswith(path, base):
+    return path.startswith(base) and (
+        len(base) == len(path) or path[len(base)] == os.path.sep )
+
+
+class PathAliaser(object):
+    def __init__(self):
+        self.aliases = {}
+        self.master_targets = set()
+        self.preferred_name = {}
+
+    def master_path(self, path):
+        match_found = False
+        while True:
+            for base, alias in self.aliases.items():
+                if path_startswith(path, base):
+                    path = alias + path[len(base):]
+                    match_found = True
+                    break
+            for master_base in self.master_targets:
+                if path_startswith(path, master_base):
+                    return path, master_base, True
+            if match_found:
+                sys.stderr.write(
+                    "(ERROR) violating fundamental assumption while walking "
+                    "directory tree.\n\tPlease report this to the gcovr "
+                    "developers.\n" )
+            return path, None, match_found
+
+    def unalias_path(self, path):
+        path = resolve_symlinks(path)
+        path, master_base, known_path = self.master_path(path)
+        if not known_path:
+            return path
+        # Try and resolve the preferred name for this location
+        if master_base in self.preferred_name:
+            return self.preferred_name[master_base] + path[len(master_base):]
+        return path
+
+    def add_master_target(self, master):
+        self.master_targets.add(master)
+
+    def add_alias(self, target, master):
+        self.aliases[target] = master
+
+    def set_preferred(self, master, preferred):
+        self.preferred_name[master] = preferred
+
+aliases = PathAliaser()
+
+# This is UGLY.  Here's why: UNIX resolves symbolic links by walking the
+# entire directory structure.  What that means is that relative links
+# are always relative to the actual directory inode, and not the
+# "virtual" path that the user might have traversed (over symlinks) on
+# the way to that directory.  Here's the canonical example:
+#
+#   a / b / c / testfile
+#   a / d / e --> ../../a/b
+#   m / n --> /a
+#   x / y / z --> /m/n/d
+#
+# If we start in "y", we will see the following directory structure:
+#   y
+#   |-- z
+#       |-- e
+#           |-- c
+#               |-- testfile
+#
+# The problem is that using a simple traversal based on the Python
+# documentation:
+#
+#    (os.path.join(os.path.dirname(path), os.readlink(result)))
+#
+# will not work: we will see a link to /m/n/d from /x/y, but completely
+# miss the fact that n is itself a link.  If we then naively attempt to
+# apply the "c" relative link, we get an intermediate path that looks
+# like "/m/n/d/e/../../a/b", which would get normalized to "/m/n/a/b"; a
+# nonexistant path.  The solution is that we need to walk the original
+# path, along with the full path of all links 1 directory at a time and
+# check for embedded symlinks.
+#
+def link_walker(path):
+    targets = [os.path.abspath(path)]
+    while targets:
+        target_dir = targets.pop(0)
+        actual_dir = resolve_symlinks(target_dir)
+        #print "target dir: %s  (%s)" % (target_dir, actual_dir)
+        master_name, master_base, visited = aliases.master_path(actual_dir)
+        if visited:
+            #print "  ...root already visited as %s" % master_name
+            aliases.add_alias(target_dir, master_name)
+            continue
+        if master_name != target_dir:
+            aliases.set_preferred(master_name, target_dir)
+            aliases.add_alias(target_dir, master_name)
+        aliases.add_master_target(master_name)
+        #print "  ...master name = %s" % master_name
+        #print "  ...walking %s" % target_dir
+        for root, dirs, files in os.walk(target_dir, topdown=True):
+            #print "    ...reading %s" % root
+            for d in dirs:
+                tmp = os.path.abspath(os.path.join(root, d))
+                #print "    ...checking %s" % tmp
+                if os.path.islink(tmp):
+                    #print "      ...buffering link %s" % tmp
+                    targets.append(tmp)
+            yield root, dirs, files
+
+
+def search_file(expr, path):
+    """
+    Given a search path, recursively descend to find files that match a
+    regular expression.
+    """
+    ans = []
+    pattern = re.compile(expr)
+    if path is None or path == ".":
+        path = os.getcwd()
+    elif not os.path.exists(path):
+        raise IOError("Unknown directory '"+path+"'")
+    for root, dirs, files in link_walker(path):
+        for name in files:
+            if pattern.match(name):
+                name = os.path.join(root,name)
+                if os.path.islink(name):
+                    ans.append( os.path.abspath(os.readlink(name)) )
+                else:
+                    ans.append( os.path.abspath(name) )
+    return ans
+
+
+#
+# Get the list of datafiles in the directories specified by the user
+#
+def get_datafiles(flist, options):
+    allfiles=[]
+    for dir in flist:
+        if options.verbose:
+            sys.stdout.write( "Scanning directory %s for gcda/gcno files...\n"
+                              % (dir, ) )
+        files = search_file(".*\.gc(da|no)$", dir)
+        # gcno files will *only* produce uncovered results; however,
+        # that is useful information for the case where a compilation
+        # unit is never actually exercised by the test code.  So, we
+        # will process gcno files, but ONLY if there is no corresponding
+        # gcda file.
+        gcda_files = [file for file in files if file.endswith('gcda')]
+        tmp = set(gcda_files)
+        gcno_files = [ file for file in files if
+                       file.endswith('gcno') and file[:-2]+'da' not in tmp ]
+        if options.verbose:
+            sys.stdout.write(
+                "Found %d files (and will process %d)\n" %
+                ( len(files), len(gcda_files) + len(gcno_files) ) )
+        allfiles.extend(gcda_files)
+        allfiles.extend(gcno_files)
+    return allfiles
+
+
+def process_gcov_data(file, covdata, options):
+    INPUT = open(file,"r")
+    #
+    # Get the filename
+    #
+    line = INPUT.readline()
+    segments=line.split(':',3)
+    if len(segments) != 4 or not segments[2].lower().strip().endswith('source'):
+        raise RuntimeError('Fatal error parsing gcov file, line 1: \n\t"%s"' % line.rstrip())
+    fname = aliases.unalias_path(os.path.abspath((segments[-1]).strip()))
+    if options.verbose:
+        sys.stdout.write("Parsing coverage data for file %s\n" % fname)
+    #
+    # Return if the filename does not match the filter
+    #
+    if not options.filter.match(fname):
+        if options.verbose:
+            sys.stdout.write("  Filtering coverage data for file %s\n" % fname)
+        return
+    #
+    # Return if the filename matches the exclude pattern
+    #
+    for i in range(0,len(options.exclude)):
+        if options.exclude[i].match(options.filter.sub('',fname)) or \
+               options.exclude[i].match(fname) or \
+               options.exclude[i].match(os.path.abspath(fname)):
+            if options.verbose:
+                sys.stdout.write("  Excluding coverage data for file %s\n" % fname)
+            return
+    #
+    # Parse each line, and record the lines
+    # that are uncovered
+    #
+    noncode   = set()
+    uncovered = set()
+    covered   = {}
+    branches  = {}
+    #first_record=True
+    lineno = 0
+    for line in INPUT:
+        segments=line.split(":",2)
+        tmp = segments[0].strip()
+        if len(segments) > 1:
+            try:
+                lineno = int(segments[1].strip())
+            except:
+                pass # keep previous line number!
+
+        if tmp[0] == '#':
+            uncovered.add( lineno )
+        elif tmp[0] in "0123456789":
+            covered[lineno] = int(segments[0].strip())
+        elif tmp[0] == '-':
+            # remember certain non-executed lines
+            code = segments[2].strip()
+            if len(code) == 0 or code == "{" or code == "}" or \
+               code.startswith("//") or code == 'else':
+                noncode.add( lineno )
+        elif tmp.startswith('branch'):
+            fields = line.split()
+            try:
+                count = int(fields[3])
+                branches.setdefault(lineno, {})[int(fields[1])] = count
+            except:
+                # We ignore branches that were "never executed"
+                pass
+        elif tmp.startswith('call'):
+            pass
+        elif tmp.startswith('function'):
+            pass
+        elif tmp[0] == 'f':
+            pass
+            #if first_record:
+                #first_record=False
+                #uncovered.add(prev)
+            #if prev in uncovered:
+                #tokens=re.split('[ \t]+',tmp)
+                #if tokens[3] != "0":
+                    #uncovered.remove(prev)
+            #prev = int(segments[1].strip())
+            #first_record=True
+        else:
+            sys.stderr.write(
+                "(WARNING) Unrecognized GCOV output: '%s'\n"
+                "\tThis is indicitive of a gcov output parse error.\n"
+                "\tPlease report this to the gcovr developers." % tmp )
+    #
+    # If the file is already in covdata, then we
+    # remove lines that are covered here.  Otherwise,
+    # initialize covdata
+    #
+    if not fname in covdata:
+        covdata[fname] = CoverageData(fname,uncovered,covered,branches,noncode)
+    else:
+        covdata[fname].update(uncovered,covered,branches,noncode)
+    INPUT.close()
+
+#
+# Process a datafile (generated by running the instrumented application)
+# and run gcov with the corresponding arguments
+#
+# This is trickier than it sounds: The gcda/gcno files are stored in the
+# same directory as the object files; however, gcov must be run from the
+# same directory where gcc/g++ was run.  Normally, the user would know
+# where gcc/g++ was invoked from and could tell gcov the path to the
+# object (and gcda) files with the --object-directory command.
+# Unfortunately, we do everything backwards: gcovr looks for the gcda
+# files and then has to infer the original gcc working directory.
+#
+# In general, (but not always) we can assume that the gcda file is in a
+# subdirectory of the original gcc working directory, so we will first
+# try ".", and on error, move up the directory tree looking for the
+# correct working directory (letting gcov's own error codes dictate when
+# we hit the right directory).  This covers 90+% of the "normal" cases.
+# The exception to this is if gcc was invoked with "-o ../[...]" (i.e.,
+# the object directory was a peer (not a parent/child) of the cwd.  In
+# this case, things are really tough.  We accept an argument
+# (--object-directory) that SHOULD BE THE SAME as the one povided to
+# gcc.  We will then walk that path (backwards) in the hopes of
+# identifying the original gcc working directory (there is a bit of
+# trial-and-error here)
+#
+def process_datafile(filename, covdata, options):
+    #
+    # Launch gcov
+    #
+    abs_filename = os.path.abspath(filename)
+    (dirname,fname) = os.path.split(abs_filename)
+    #(name,ext) = os.path.splitext(base)
+
+    potential_wd = []
+    starting_dir = os.getcwd()
+    errors=[]
+    Done = False
+
+    if options.objdir:
+        src_components = abs_filename.split(os.sep)
+        components = normpath(options.objdir).split(os.sep)
+        idx = 1
+        while idx <= len(components):
+            if idx > len(src_components):
+                break
+            if components[-1*idx] != src_components[-1*idx]:
+                break
+            idx += 1
+        if idx > len(components):
+            pass # a parent dir; the normal process will find it
+        elif components[-1*idx] == '..':
+            dirs = [ os.path.join(src_components[:len(src_components)-idx+1]) ]
+            while idx <= len(components) and components[-1*idx] == '..':
+                tmp = []
+                for d in dirs:
+                    for f in os.listdir(d):
+                        x = os.path.join(d,f)
+                        if os.path.isdir(x):
+                            tmp.append(x)
+                dirs = tmp
+                idx += 1
+            potential_wd = dirs
+        else:
+            if components[0] == '':
+                # absolute path
+                tmp = [ options.objdir ]
+            else:
+                # relative path: check relative to both the cwd and the
+                # gcda file
+                tmp = [ os.path.join(x, options.objdir) for x in
+                        [os.path.dirname(abs_filename), os.getcwd()] ]
+            potential_wd = [ testdir for testdir in tmp
+                             if os.path.isdir(testdir) ]
+            if len(potential_wd) == 0:
+                errors.append("ERROR: cannot identify the location where GCC "
+                              "was run using --object-directory=%s\n" %
+                              options.objdir)
+            # Revert to the normal
+            #sys.exit(1)
+
+    # no objdir was specified (or it was a parent dir); walk up the dir tree
+    if len(potential_wd) == 0:
+        wd = os.path.split(abs_filename)[0]
+        while True:
+            potential_wd.append(wd)
+            wd = os.path.split(wd)[0]
+            if wd == potential_wd[-1]:
+                break
+
+    cmd = [ gcov_cmd, abs_filename,
+            "--branch-counts", "--branch-probabilities", "--preserve-paths",
+            '--object-directory', dirname ]
+
+    # NB: We are lazy English speakers, so we will only parse English output
+    env = dict(os.environ)
+    env['LC_ALL'] = 'en_US'
+
+
+    while len(potential_wd) > 0 and not Done:
+        # NB: either len(potential_wd) == 1, or all entires are absolute
+        # paths, so we don't have to chdir(starting_dir) at every
+        # iteration.
+        os.chdir(potential_wd.pop(0))
+
+
+        #if options.objdir:
+        #    cmd.extend(["--object-directory", Template(options.objdir).substitute(filename=filename, head=dirname, tail=base, root=name, ext=ext)])
+
+        if options.verbose:
+            sys.stdout.write("Running gcov: '%s' in '%s'\n" % ( ' '.join(cmd), os.getcwd() ))
+        (out, err) = subprocess.Popen( cmd, env=env,
+                                       stdout=subprocess.PIPE,
+                                       stderr=subprocess.PIPE ).communicate()
+        out=out.decode('utf-8')
+        err=err.decode('utf-8')
+
+        # find the files that gcov created
+        gcov_files = {'active':[], 'filter':[], 'exclude':[]}
+        for line in out.splitlines():
+            found = output_re.search(line.strip())
+            if found is not None:
+                fname = found.group(1)
+                if not options.gcov_filter.match(fname):
+                    if options.verbose:
+                        sys.stdout.write("Filtering gcov file %s\n" % fname)
+                    gcov_files['filter'].append(fname)
+                    continue
+                exclude=False
+                for i in range(0,len(options.gcov_exclude)):
+                    if options.gcov_exclude[i].match(options.gcov_filter.sub('',fname)) or \
+                           options.gcov_exclude[i].match(fname) or \
+                           options.gcov_exclude[i].match(os.path.abspath(fname)):
+                        exclude=True
+                        break
+                if not exclude:
+                    gcov_files['active'].append(fname)
+                elif options.verbose:
+                    sys.stdout.write("Excluding gcov file %s\n" % fname)
+                    gcov_files['exclude'].append(fname)
+
+        if source_re.search(err):
+            # gcov tossed errors: try the next potential_wd
+            errors.append(err)
+        else:
+            # Process *.gcov files
+            for fname in gcov_files['active']:
+                process_gcov_data(fname, covdata, options)
+            Done = True
+
+        if not options.keep:
+            for group in gcov_files.values():
+                for fname in group:
+                    os.remove(fname)
+
+    os.chdir(starting_dir)
+    if options.delete:
+        if not abs_filename.endswith('gcno'):
+            os.remove(abs_filename)
+
+    if not Done:
+        sys.stderr.write(
+            "(WARNING) GCOV produced the following errors processing %s:\n"
+            "\t   %s"
+            "\t(gcovr could not infer a working directory that resolved it.)\n"
+            % ( filename, "\t   ".join(errors) ) )
+
+#
+# Produce the classic gcovr text report
+#
+def print_text_report(covdata):
+    def _num_uncovered(key):
+        (total, covered, percent) = covdata[key].coverage()
+        return total - covered
+    def _percent_uncovered(key):
+        (total, covered, percent) = covdata[key].coverage()
+        if covered:
+            return -1.0*covered/total
+        else:
+            return total or 1e6
+    def _alpha(key):
+        return key
+
+    if options.output:
+        OUTPUT = open(options.output,'w')
+    else:
+        OUTPUT = sys.stdout
+    total_lines=0
+    total_covered=0
+    # Header
+    OUTPUT.write("-"*78 + '\n')
+    a = options.show_branch and "Branch" or "Lines"
+    b = options.show_branch and "Taken" or "Exec"
+    OUTPUT.write("File".ljust(40) + a.rjust(8) + b.rjust(8)+ "  Cover   Missing\n")
+    OUTPUT.write("-"*78 + '\n')
+
+    # Data
+    keys = list(covdata.keys())
+    keys.sort(key=options.sort_uncovered and _num_uncovered or \
+              options.sort_percent and _percent_uncovered or _alpha)
+    for key in keys:
+        (t, n, txt) = covdata[key].summary()
+        total_lines += t
+        total_covered += n
+        OUTPUT.write(txt + '\n')
+
+    # Footer & summary
+    OUTPUT.write("-"*78 + '\n')
+    percent = total_lines and str(int(100.0*total_covered/total_lines)) or "--"
+    OUTPUT.write("TOTAL".ljust(40) + str(total_lines).rjust(8) + \
+          str(total_covered).rjust(8) + str(percent).rjust(6)+"%" + '\n')
+    OUTPUT.write("-"*78 + '\n')
+
+    # Close logfile
+    if options.output:
+        OUTPUT.close()
+
+#
+# Produce an XML report in the Cobertura format
+#
+def print_xml_report(covdata):
+    branchTotal = 0
+    branchCovered = 0
+    lineTotal = 0
+    lineCovered = 0
+
+    options.show_branch = True
+    for key in covdata.keys():
+        (total, covered, percent) = covdata[key].coverage()
+        branchTotal += total
+        branchCovered += covered
+
+    options.show_branch = False
+    for key in covdata.keys():
+        (total, covered, percent) = covdata[key].coverage()
+        lineTotal += total
+        lineCovered += covered
+
+    impl = xml.dom.minidom.getDOMImplementation()
+    docType = impl.createDocumentType(
+        "coverage", None,
+        "http://cobertura.sourceforge.net/xml/coverage-03.dtd" )
+    doc = impl.createDocument(None, "coverage", docType)
+    root = doc.documentElement
+    root.setAttribute( "line-rate", lineTotal == 0 and '0.0' or
+                       str(float(lineCovered) / lineTotal) )
+    root.setAttribute( "branch-rate", branchTotal == 0 and '0.0' or
+                       str(float(branchCovered) / branchTotal) )
+    root.setAttribute( "timestamp", str(int(time.time())) )
+    root.setAttribute( "version", "gcovr %s" % (version_str(),) )
+
+    # Generate the <sources> element: this is either the root directory
+    # (specified by --root), or the CWD.
+    sources = doc.createElement("sources")
+    root.appendChild(sources)
+
+    # Generate the coverage output (on a per-package basis)
+    packageXml = doc.createElement("packages")
+    root.appendChild(packageXml)
+    packages = {}
+    source_dirs = set()
+
+    keys = list(covdata.keys())
+    keys.sort()
+    for f in keys:
+        data = covdata[f]
+        dir = options.filter.sub('',f)
+        if f.endswith(dir):
+            src_path = f[:-1*len(dir)]
+            if len(src_path) > 0:
+                while dir.startswith(os.path.sep):
+                    src_path += os.path.sep
+                    dir = dir[len(os.path.sep):]
+                source_dirs.add(src_path)
+        else:
+            # Do no truncation if the filter does not start matching at
+            # the beginning of the string
+            dir = f
+        (dir, fname) = os.path.split(dir)
+
+        package = packages.setdefault(
+            dir, [ doc.createElement("package"), {},
+                   0, 0, 0, 0 ] )
+        c = doc.createElement("class")
+        lines = doc.createElement("lines")
+        c.appendChild(lines)
+
+        class_lines = 0
+        class_hits = 0
+        class_branches = 0
+        class_branch_hits = 0
+        for line in data.all_lines:
+            hits = data.covered.get(line, 0)
+            class_lines += 1
+            if hits > 0:
+                class_hits += 1
+            l = doc.createElement("line")
+            l.setAttribute("number", str(line))
+            l.setAttribute("hits", str(hits))
+            branches = data.branches.get(line)
+            if branches is None:
+                l.setAttribute("branch", "false")
+            else:
+                b_hits = 0
+                for v in branches.values():
+                    if v > 0:
+                        b_hits += 1
+                coverage = 100*b_hits/len(branches)
+                l.setAttribute("branch", "true")
+                l.setAttribute( "condition-coverage",
+                                "%i%% (%i/%i)" %
+                                (coverage, b_hits, len(branches)) )
+                cond = doc.createElement('condition')
+                cond.setAttribute("number", "0")
+                cond.setAttribute("type", "jump")
+                cond.setAttribute("coverage", "%i%%" % ( coverage ) )
+                class_branch_hits += b_hits
+                class_branches += float(len(branches))
+                conditions = doc.createElement("conditions")
+                conditions.appendChild(cond)
+                l.appendChild(conditions)
+
+            lines.appendChild(l)
+
+        className = fname.replace('.', '_')
+        c.setAttribute("name", className)
+        c.setAttribute("filename", os.path.join(dir, fname))
+        c.setAttribute("line-rate", str(class_hits / (1.0*class_lines or 1.0)))
+        c.setAttribute( "branch-rate",
+                        str(class_branch_hits / (1.0*class_branches or 1.0)) )
+        c.setAttribute("complexity", "0.0")
+
+        package[1][className] = c
+        package[2] += class_hits
+        package[3] += class_lines
+        package[4] += class_branch_hits
+        package[5] += class_branches
+
+    for packageName, packageData in packages.items():
+        package = packageData[0];
+        packageXml.appendChild(package)
+        classes = doc.createElement("classes")
+        package.appendChild(classes)
+        classNames = list(packageData[1].keys())
+        classNames.sort()
+        for className in classNames:
+            classes.appendChild(packageData[1][className])
+        package.setAttribute("name", packageName.replace(os.sep, '.'))
+        package.setAttribute("line-rate", str(packageData[2]/(1.0*packageData[3] or 1.0)))
+        package.setAttribute( "branch-rate", str(packageData[4] / (1.0*packageData[5] or 1.0) ))
+        package.setAttribute("complexity", "0.0")
+
+
+    # Populate the <sources> element: this is either the root directory
+    # (specified by --root), or relative directories based
+    # on the filter, or the CWD
+    if options.root is not None:
+        source = doc.createElement("source")
+        source.appendChild(doc.createTextNode(options.root))
+        sources.appendChild(source)
+    elif len(source_dirs) > 0:
+        cwd = os.getcwd()
+        for d in source_dirs:
+            source = doc.createElement("source")
+            if d.startswith(cwd):
+                reldir = d[len(cwd):].lstrip(os.path.sep)
+            elif cwd.startswith(d):
+                i = 1
+                while normpath(d) != \
+                          normpath(os.path.join(*tuple([cwd]+['..']*i))):
+                    i += 1
+                reldir = os.path.join(*tuple(['..']*i))
+            else:
+                reldir = d
+            source.appendChild(doc.createTextNode(reldir))
+            sources.appendChild(source)
+    else:
+        source = doc.createElement("source")
+        source.appendChild(doc.createTextNode('.'))
+        sources.appendChild(source)
+
+    xmlString = doc.toprettyxml()
+    #xml.dom.ext.PrettyPrint(doc)
+    if options.output is None:
+        sys.stdout.write(xmlString+'\n')
+    else:
+        OUTPUT = open(options.output, 'w')
+        OUTPUT.write(xmlString +'\n')
+        OUTPUT.close()
+
+
+##
+## MAIN
+##
+
+#
+# Create option parser
+#
+parser = OptionParser()
+parser.add_option("--version",
+        help="Print the version number, then exit",
+        action="store_true",
+        dest="version",
+        default=False)
+parser.add_option("-v","--verbose",
+        help="Print progress messages",
+        action="store_true",
+        dest="verbose",
+        default=False)
+parser.add_option('--object-directory',
+        help="Specify the directory that contains the gcov data files.  gcovr must be able to identify the path between the *.gcda files and the directory where gcc was originally run.  Normally, gcovr can guess correctly.  This option overrides gcovr's normal path detection and can specify either the path from gcc to the gcda file (i.e. what was passed to gcc's '-o' option), or the path from the gcda file to gcc's original working directory.",
+        action="store",
+        dest="objdir",
+        default=None)
+parser.add_option("-o","--output",
+        help="Print output to this filename",
+        action="store",
+        dest="output",
+        default=None)
+parser.add_option("-k","--keep",
+        help="Keep temporary gcov files",
+        action="store_true",
+        dest="keep",
+        default=False)
+parser.add_option("-d","--delete",
+        help="Delete the coverage files after they are processed",
+        action="store_true",
+        dest="delete",
+        default=False)
+parser.add_option("-f","--filter",
+        help="Keep only the data files that match this regular expression",
+        action="store",
+        dest="filter",
+        default=None)
+parser.add_option("-e","--exclude",
+        help="Exclude data files that match this regular expression",
+        action="append",
+        dest="exclude",
+        default=[])
+parser.add_option("--gcov-filter",
+        help="Keep only gcov data files that match this regular expression",
+        action="store",
+        dest="gcov_filter",
+        default=None)
+parser.add_option("--gcov-exclude",
+        help="Exclude gcov data files that match this regular expression",
+        action="append",
+        dest="gcov_exclude",
+        default=[])
+parser.add_option("-r","--root",
+        help="Defines the root directory.  This is used to filter the files, and to standardize the output.",
+        action="store",
+        dest="root",
+        default=None)
+parser.add_option("-x","--xml",
+        help="Generate XML instead of the normal tabular output.",
+        action="store_true",
+        dest="xml",
+        default=None)
+parser.add_option("-b","--branches",
+        help="Tabulate the branch coverage instead of the line coverage.",
+        action="store_true",
+        dest="show_branch",
+        default=None)
+parser.add_option("-u","--sort-uncovered",
+        help="Sort entries by increasing number of uncovered lines.",
+        action="store_true",
+        dest="sort_uncovered",
+        default=None)
+parser.add_option("-p","--sort-percentage",
+        help="Sort entries by decreasing percentage of covered lines.",
+        action="store_true",
+        dest="sort_percent",
+        default=None)
+parser.usage="gcovr [options]"
+parser.description="A utility to run gcov and generate a simple report that summarizes the coverage"
+#
+# Process options
+#
+(options, args) = parser.parse_args(args=sys.argv)
+if options.version:
+    sys.stdout.write(
+        "gcovr %s\n"
+        "\n"
+        "Copyright (2008) Sandia Corporation. Under the terms of Contract\n"
+        "DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government\n"
+        "retains certain rights in this software.\n"
+        % (version_str(),) )
+    sys.exit(0)
+if options.objdir:
+    tmp = options.objdir.replace('/',os.sep).replace('\\',os.sep)
+    while os.sep+os.sep in tmp:
+        tmp = tmp.replace(os.sep+os.sep, os.sep)
+    if normpath(options.objdir) != tmp:
+        sys.stderr.write(
+            "(WARNING) relative referencing in --object-directory.\n"
+            "\tthis could cause strange errors when gcovr attempts to\n"
+            "\tidentify the original gcc working directory.\n")
+#
+# Setup filters
+#
+for i in range(0,len(options.exclude)):
+    options.exclude[i] = re.compile(options.exclude[i])
+if options.filter is not None:
+    options.filter = re.compile(options.filter)
+elif options.root is not None:
+    if not options.root:
+        sys.stderr.write(
+            "(ERROR) empty --root option.\n"
+            "\tRoot specifies the path to the root directory of your project\n"
+            "\tand cannot be an empty string.\n")
+        sys.exit(1)
+    options.filter = re.compile(re.escape(os.path.abspath(options.root)+os.sep))
+if options.filter is None:
+    options.filter = re.compile('')
+#
+for i in range(0,len(options.gcov_exclude)):
+    options.gcov_exclude[i] = re.compile(options.gcov_exclude[i])
+if options.gcov_filter is not None:
+    options.gcov_filter = re.compile(options.gcov_filter)
+else:
+    options.gcov_filter = re.compile('')
+#
+# Get data files
+#
+if len(args) == 1:
+    datafiles = get_datafiles(["."], options)
+else:
+    datafiles = get_datafiles(args[1:], options)
+#
+# Get coverage data
+#
+covdata = {}
+for file in datafiles:
+    process_datafile(file,covdata,options)
+if options.verbose:
+    sys.stdout.write("Gathered coveraged data for "+str(len(covdata))+" files\n")
+#
+# Print report
+#
+if options.xml:
+    print_xml_report(covdata)
+else:
+    print_text_report(covdata)